# Simple curve fitting with CurveFit (or alternative)

**URL:** <https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737>\
**Category:** General Usage\
**Created:** [April 7, 2021, 8:34am UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737 "2021-04-07T08:34:21Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![StatisticalMouse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/statisticalmouse/32/43370_2.png) [@StatisticalMouse](https://discourse.julialang.org/u/StatisticalMouse)\
**Post date:** [April 7, 2021, 8:34am UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/1 "2021-04-07T08:34:21Z")

</div>

I’d like to try different curves against data to see which of them would best fit. So in this example

```julia
x0 = collect(1.0:100.0)
y0 = x .^ 4.5
f0 = curve_fit(LinearFit, x0, y0)

```

I’d like to learn just how bad it is, but is there an easy way to get that out of it?

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [April 7, 2021, 10:22am UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/2 "2021-04-07T10:22:05Z")

</div>

Please check this [post](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/38).

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [April 7, 2021, 10:40am UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/3 "2021-04-07T10:40:36Z")

</div>

Just in case it might be useful, sharing one linear regression example below, where the correlation coefficient **r** is displayed and other regression statistics are computed too:

```julia
using Plots, Random, Distributions, GLM, DataFrames, Printf

# input data
n = 30
x = LinRange(0,1, n) .* 10
e = rand(Normal(0,mean(x)/4),n)
parms = Dict(:b0 => 6, :b1 => 3.2)
y = parms[:b0] .+ parms[:b1] .* x .+ 3*e
data = DataFrame(X=x, Y=y)

# linear regression
ols = lm(@formula(Y ~ 1 + X), data)
a, b = coef(ols)
r = r2(ols)

# plotting
yhat = parms[:b0] .+ parms[:b1] .* x
err = y - yhat; z0 = 0*err;
str = "Y = " * @sprintf("%.2f",a) * " + " * @sprintf("%.2f",b) * " X (r = " * @sprintf("%.2f",r) * ")" 
plot(x,yhat, label=str, lw=2, c=:green, yerror=(z0,err), ribbon=(z0,err), fc=:cyan, fa=0.1)
scatter!(x, y, mc=:red, xlabel="X", ylabel="Y", label="Input data: Y = Y(X)", legend=:bottomright)

```

Some stats output after command:

```julia
julia> ols = lm(@formula(Y ~ 1 + X), data)

Y ~ 1 + X

Coefficients:
───────────────────────────────────────────────────────────────────────
               Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
───────────────────────────────────────────────────────────────────────
(Intercept) 5.16706 1.31032 3.94 0.0005 2.48298 7.85113
X 3.172 0.225023 14.10 <1e-13 2.71106 3.63294

```

![linear_regression_with_corr_coeff](https://global.discourse-cdn.com/julialang/original/3X/f/e/fe0e03622d5b92054c2dddb98707f7f5d9d381e8.png)

---

<div class="post-metadata">

**Author:** ![StatisticalMouse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/statisticalmouse/32/43370_2.png) [@StatisticalMouse](https://discourse.julialang.org/u/StatisticalMouse)\
**Post date:** [April 7, 2021, 12:04pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/4 "2021-04-07T12:04:31Z")

</div>

Thanks, I’ll look into these examples. The problem itself should be relatively simple, the actual data is a straight line in a log-log plot, it’s just not quite that easy to post its code here. (The wikipedia page for log-log plots is also somewhat opaque.)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [April 7, 2021, 12:47pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/5 "2021-04-07T12:47:15Z")

</div>

Check out this document too on [lin-log and log-log regressions](https://kenbenoit.net/assets/courses/ME104/logmodels2.pdf).

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [April 7, 2021, 12:58pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/6 "2021-04-07T12:58:05Z")

</div>

I have a very simple package to make fits really, really simple, in common cases: [GitHub - m3g/EasyFit.jl: Easy interface for obtaining fits for 2D data](https://github.com/m3g/EasyFit.jl)

---

<div class="post-metadata">

**Author:** ![feanor12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/feanor12/32/8212_2.png) [@feanor12](https://discourse.julialang.org/u/feanor12)\
**Post date:** [April 7, 2021, 1:36pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/7 "2021-04-07T13:36:05Z")

</div>

Just for fun:  
To fit a polynom you can use

```julia
"""
Fit y = f(x) to nth order polynomial.
Returns list of polynomial coefficients.
"""
function polyfit(x,y,n)
	@assert length(x) > n
	A = [k^n for k in x, n in 0:n]
	A\y
end

```

To reason about the fit you would need to look for measures often found under [goodness of fit](https://en.wikipedia.org/wiki/Goodness_of_fit), but I am not sure which are implemented by julia packages. If you want to look into one of those methods I think the package [HypothesisTests](https://juliastats.org/HypothesisTests.jl/stable/) is a good resource.

Also if you try to fit a [log transformed exponential function](https://mathworld.wolfram.com/LeastSquaresFittingExponential.html) you should use a weighted least square fit to compensate for the transform.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 8, 2021, 11:54am UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/8 "2021-04-08T11:54:40Z")

</div>

> [@feanor12](#):
>
> `A = [k^n for k in x, n in 0:n]`

Don’t do this — the [Vandermonde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix) is generally [ill-conditioned](https://arxiv.org/abs/1504.02118).

---

<div class="post-metadata">

**Author:** ![StatisticalMouse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/statisticalmouse/32/43370_2.png) [@StatisticalMouse](https://discourse.julialang.org/u/StatisticalMouse)\
**Post date:** [April 8, 2021, 12:10pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/9 "2021-04-08T12:10:06Z")

</div>

I got this working with minimal effort (on the package side, I had some other unrelated Julia problems though). When I take log10 of the original x and y values, and do `fit = fitlinear(x_log,y_log)`, I get a nice match.

 ![Screenshot 2021-04-08 at 15.09.00](https://global.discourse-cdn.com/julialang/original/3X/d/5/d59b30084d5b89f556c7dbd5c852f05301be206d.png)  
I can also see that the average square residue is low.

```julia
 ------------------- Linear Fit ------------- 

 Equation: y = ax + b 

 With: a = -0.9999071242153111
       b = 5.319512210970783e-5

 Pearson correlation coefficient, R = 0.9999734052571524
 Average square residue = 9.731441862500589e-6

 Predicted Y: ypred = [2.999774567768043, 2.6987725305011243...
 residues = [-0.0008626824944601985, -0.007119649365998182...

 -------------------------------------------- 

```

---

<div class="post-metadata">

**Author:** ![StatisticalMouse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/statisticalmouse/32/43370_2.png) [@StatisticalMouse](https://discourse.julialang.org/u/StatisticalMouse)\
**Post date:** [April 8, 2021, 12:13pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/10 "2021-04-08T12:13:17Z")

</div>

When I tried to fit against the original (non-log) data, the only one that gave good results was `fit = fitspline(x,y)`.

 ![Screenshot 2021-04-08 at 15.11.29](https://global.discourse-cdn.com/julialang/original/3X/7/8/780bc2d97e5d956c148ff8554924e317dac1958a.png)  
It’s just that spline doesn’t really tell me anything useful.

```julia
 -------- Spline fit --------------------------- 

 x spline: x = [0.001, 0.011...
 y spline: y = [1005.2007999999998, 90.34990000000002...

 ----------------------------------------------- 

```

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [April 8, 2021, 1:46pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/11 "2021-04-08T13:46:45Z")

</div>

A spline will always look well, but you cannot use them to really get any information about your data.

If your `x` and `y` values are related by something like `y = A/x` it is expected that the log/log plot is linear (`log(y) = A - log(x)`) while the plot of the original data is not.

---

<div class="post-metadata">

**Author:** ![StatisticalMouse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/statisticalmouse/32/43370_2.png) [@StatisticalMouse](https://discourse.julialang.org/u/StatisticalMouse)\
**Post date:** [April 8, 2021, 2:23pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/12 "2021-04-08T14:23:13Z")

</div>

Yes. X ^ -1 is close. Can I fit that with EasyFit?

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [April 8, 2021, 2:40pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/13 "2021-04-08T14:40:10Z")

</div>

No, I didn’t implement that (you can fit anything using directly [LsqFit](https://github.com/JuliaNLSolvers/LsqFit.jl)). Although the approach of fitting a straight line in the log/log plot is the best approach. Linear fits are always better.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [April 8, 2021, 3:00pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/14 "2021-04-08T15:00:01Z")

</div>

LsqFit seems to be more general and conveniently provides the confidence intervals for the parameters estimated.  
A correlation coefficient might also be obtained from the approximation of the covariance matrix. An attempt below, `r` formula to be confirmed.

_PS: in your case it is possible and better to work with a log-log scale._

```julia
using LsqFit, Plots, Printf

# input noisy data
t = 0.01:0.01:1
y = 420*t.^(-0.13 .+ 0.02*rand(length(t)))

# LsqFit
m(t, p) = p[1] * t.^(-p[2])
p0 = [0.5, 0.5]
fit = curve_fit(m, t, y, p0)
# fit parameters and confidence interval
p = fit.param
confidence_interval(fit, 0.1) # 10% confidence intervals

# correlation coefficient = covariance of 2 variables / product of their standard deviations
cov = estimate_covar(fit)
r = cov[1,2]/sqrt(cov[1,1]*cov[2,2]) # TODO: confirm formula

# plot results
str = @sprintf("%.1f",p[1])*"*t^(-" * @sprintf("%.3f",p[2]) * ")" * @sprintf("; r ~ %.2f",r)
scatter(t,y, mc=:reds, ms=3, xlabel="t", label="input data")
plot!(t,m(t,fit.param),ly=0.25,lc=:blue, ylabel="y(t)", label="LsqFit: "*str)

```

```julia
julia> p = fit.param
 419.2404632519413
   0.12214587254371537

julia> confidence_interval(fit, 0.1) # 10% confidence intervals
 (418.14933388256867, 420.331592621314)
 (0.12055572549189375, 0.12373601959553698)

julia> r = cov[1,2]/sqrt(cov[1,1]*cov[2,2]) # TODO: confirm formula
-0.7431633567883187

```

![power_law_LsqFit](https://global.discourse-cdn.com/julialang/original/3X/e/6/e62c2665604a9e001e146d7cc641f743851fade4.png)

---

<div class="post-metadata">

**Author:** ![StatisticalMouse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/statisticalmouse/32/43370_2.png) [@StatisticalMouse](https://discourse.julialang.org/u/StatisticalMouse)\
**Post date:** [April 8, 2021, 3:29pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/15 "2021-04-08T15:29:36Z")

</div>

After thinking about this, I agree. The linear plot is really hard to read if you want any useful info out of it. The only reason against the log-log plot is that I’m not so used to such plots. Well, time to learn.

I’ll take a look at the solutions a bit later.

---

<div class="post-metadata">

**Author:** ![StatisticalMouse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/statisticalmouse/32/43370_2.png) [@StatisticalMouse](https://discourse.julialang.org/u/StatisticalMouse)\
**Post date:** [April 8, 2021, 3:34pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/16 "2021-04-08T15:34:51Z")

</div>

Also, if this was a part of a machine learning problem, trying gradient descent on the linear data would be a lost cause. I ’knew’ that already, but visualizing it makes is clearer.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [April 8, 2021, 3:45pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/17 "2021-04-08T15:45:53Z")

</div>

If you can linearize the plot, you probably should do the linear fit and afterwards recover the non-linear parameters from that fit.

Now, how you maybe using this for optimizing anything is another story. The function `1/x` is also minimum at `x -> infty`, thus minimizing it or the linearized version will take you to the same place.

---

<div class="post-metadata">

**Author:** ![StatisticalMouse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/statisticalmouse/32/43370_2.png) [@StatisticalMouse](https://discourse.julialang.org/u/StatisticalMouse)\
**Post date:** [April 8, 2021, 3:52pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/18 "2021-04-08T15:52:11Z")

</div>

How would I get the lin-lin equation from this log-log fit? I spent some reading of the wikipedia page, but not enough to fully get it.

---

<div class="post-metadata">

**Author:** ![StatisticalMouse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/statisticalmouse/32/43370_2.png) [@StatisticalMouse](https://discourse.julialang.org/u/StatisticalMouse)\
**Post date:** [April 8, 2021, 7:16pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/19 "2021-04-08T19:16:46Z")

</div>

Ok, got it. To change from log-log to lin-lin one can do:  
`y₂ = MathConstants.e .^ ( a .* log.(x) .+ b .* log.(10))`

I did use WolframAlpha, I fed it `log10(y) = a * log10(x) + b `  
It says the answer is approximate form. I then changed it to Julia form.

```julia
log₁₀y == a * log₁₀x + b 
10 ^ log₁₀y == 10 ^ ( a * log₁₀x + b )
y == 10 ^ ( a * log₁₀x + b )
y == 10 ^ ( a * log₁₀x ) * 10 ^ b
magic..
y == MathConstants.e .^ ( a .* log.(x) .+ b .* log.(10))

```

Above is pseudo-code.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [April 8, 2021, 7:29pm UTC](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737/20 "2021-04-08T19:29:27Z")

</div>

Indeed. That can be written in Julia as:

```julia
y2 = exp(fit.b) * x.^fit.a

```

(where `fit.b` and `fit.a` are the intercept and slope of the linear fit of the log-log plot - with base e. If you use base 10, use `log10`, and then y\_2 = 10^{fit.b} x^{fit.a}.

[Next page](https://discourse.julialang.org/t/simple-curve-fitting-with-curvefit-or-alternative/58737.md?page=2)