# Estimating parameters of a function for curve fitting

**URL:** <https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404>\
**Category:** Optimization (Mathematical)\
**Tags:** curve-fitting\
**Created:** [December 17, 2019, 7:48pm UTC](https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404 "2019-12-17T19:48:57Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![affans](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/affans/32/11911_2.png) [@affans](https://discourse.julialang.org/u/affans)\
**Post date:** [December 17, 2019, 7:48pm UTC](https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404/1 "2019-12-17T19:48:57Z")

</div>

I have some data that has missing values. I am trying to fit a curve on this data, using Richard’s Model which is

\begin{aligned} I(t)=& \frac{K}{\left[1+e^{-r a(t-\tau)}\right] \frac{1}{a}} \end{aligned} 

where K is fixed and I am trying to estimate a, \tau, and r. So for any time t, I want to be able to compute I(t) and it should be close to my data. I do have missing points in my time series.

I know how to do this in `R`, with the following code (here \tau is replaced with b.

```julia
######
richards.model <- function(prm, K, x) {
  a = prm[['a']]
  b = prm[['b']]
  r = prm[['r']]
  denom <- (1+exp(-r*a*(x-b)))^(1/a)
  res <- K / denom
  return(res)
}

errfct <- function(prm,K,dat) {
  xx <- dat$t
  yy <- richards.model(prm , K=K, x=xx)
  return( sum( (yy-dat$cum.cases)^2, na.rm = T))
}

## compute on all data points
K <- 5500
pfit <- optim(par = list(a=2,b=80,r=0.1), 
              fn = errfct, dat = dat, K=K)

```

So the `optim` function calls the `errfct` with an estimate of parameters. The `errfct` function simply calculates I(t) with the given parameters, and calculates the sum of square errors. I am assuming the `optim` function is trying to minimize the least squares error.

Is there an equivalent package in Julia?

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [December 17, 2019, 8:06pm UTC](https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404/2 "2019-12-17T20:06:08Z")

</div>

You might try [`Optim.jl`](https://julianlsolvers.github.io/Optim.jl/stable/).

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [December 17, 2019, 8:07pm UTC](https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404/3 "2019-12-17T20:07:52Z")

</div>

there is [LsqFit.jl](https://github.com/JuliaNLSolvers/LsqFit.jl) for just fitting of non linear functions to data.

```julia
using Pkg #package manager
Pkg.add("LsqFit") #adding LsqFit to your Packages
using LsqFit #including LsqFit in your enviroment
#p = [a,tau,r] #vector, this is the order of parameters im gonna use
@. model(t,p,K) = p[1]*K/(1+exp(-p[1]*p[3]*(t-p[2])) #variable K to modify later
model0(K) = (t,p) -> model(t,p,K) #this is to create a function with a fixed K
p0 = [0.5,0.5,0.5] #initial parameters
xdata = ..... #your data, in this case t
ydata = ..... # experimental I(t) data
fit = curve_fit(model0, xdata, ydata, p0)

```

---

<div class="post-metadata">

**Author:** ![affans](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/affans/32/11911_2.png) [@affans](https://discourse.julialang.org/u/affans)\
**Post date:** [December 17, 2019, 8:33pm UTC](https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404/4 "2019-12-17T20:33:51Z")

</div>

Thanks, this is exactly what I needed except it dosn’t work with missing data.

```julia
# fit 
@. richards_model(x, p) = 5500 / (1+exp(-p[3]*p[1]*(x-p[2])))^(1/p[1])
p0 = [2, 80, 0.1]
xdata = Vector(df.time)
ydata = Vector(df.cumcases)
fit = curve_fit(richards_model, xdata, ydata, p0)

```

because `typeof(ydata) = Array{Union{Missing, Int64},1}`

The error thrown is

```julia
TypeError: non-boolean (Missing) used in boolean context

```

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [December 17, 2019, 8:40pm UTC](https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404/5 "2019-12-17T20:40:59Z")

</div>

add this:

```julia
idx = findall(!ismissing,ydata) #find the indices of all non missing data in ydata
xdata2 = xdata[idx] #xdata but only in the places where ydata is not missing
ydata2 = ydata[idx] #all not missing values in ydata

```

---

<div class="post-metadata">

**Author:** ![affans](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/affans/32/11911_2.png) [@affans](https://discourse.julialang.org/u/affans)\
**Post date:** [December 17, 2019, 10:10pm UTC](https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404/6 "2019-12-17T22:10:54Z")

</div>

Thanks, that did it. There is a function in the source code that checks the health of data. I may submit a PR to check for missing data as well.

```julia
function check_data_health(xdata, ydata)

    if any(isinf, xdata) || any(isinf, ydata) || any(isnan, xdata) || any(isnan, ydata)

        error("Data contains `Inf` or `NaN` values and a fit cannot be performed")

    end

end

```

---

<div class="post-metadata">

**Author:** ![pkofod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pkofod/32/2179_2.png) [@pkofod](https://discourse.julialang.org/u/pkofod)\
**Post date:** [December 19, 2019, 3:32pm UTC](https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404/7 "2019-12-19T15:32:05Z")

</div>

> [@affans](#):
>
> Thanks, that did it. There is a function in the source code that checks the health of data. I may submit a PR to check for missing data as well.

Please do that! I’ll be sure to review it. Ideally we should just support missing values by ignoring the relevant rows, but as a first step I’d be happy to merge an updated data health check.

---

<div class="post-metadata">

**Author:** ![pkofod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pkofod/32/2179_2.png) [@pkofod](https://discourse.julialang.org/u/pkofod)\
**Post date:** [December 19, 2019, 3:33pm UTC](https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404/8 "2019-12-19T15:33:06Z")

</div>

But let me also tell you the following: I know the LM is very popular for NLS, but in my own experience, BFGS directly on the ssr can be just as good an approach, if not better.

---

<div class="post-metadata">

**Author:** ![Humphrey\_Lee](https://avatars.discourse-cdn.com/v4/letter/h/a88e4f/32.png) [@Humphrey\_Lee](https://discourse.julialang.org/u/Humphrey_Lee)\
**Post date:** [April 7, 2020, 9:39am UTC](https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404/9 "2020-04-07T09:39:55Z")

</div>

I am looking for similar curve fitting (minimization of SSE) Julia codes. Do you manage to get this working? Thanks.

---

<div class="post-metadata">

**Author:** ![pkofod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pkofod/32/2179_2.png) [@pkofod](https://discourse.julialang.org/u/pkofod)\
**Post date:** [April 7, 2020, 2:45pm UTC](https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404/10 "2020-04-07T14:45:38Z")

</div>

What’s wrong with the examples here?

---

<div class="post-metadata">

**Author:** ![Humphrey\_Lee](https://avatars.discourse-cdn.com/v4/letter/h/a88e4f/32.png) [@Humphrey\_Lee](https://discourse.julialang.org/u/Humphrey_Lee)\
**Post date:** [April 7, 2020, 11:58pm UTC](https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404/11 "2020-04-07T23:58:14Z")

</div>

This example is perfectly fine. It also showed method to tackle missing data. What I’m missing is the original question which is to translate the R code into Julia equivalent (e.g. using Optim.jl, GLM.jl or JuMP.jl). Pardon me, I am Julia novice.

Much of my experiences is from Excel solver. I also use similar approach (min SSE) in R. I did use curve\_fitting\_solver\_function, e.g. curve\_fit() function to solve equations - did that in Python using SciPy. What I do like about minimization of SSE is that I can modify the minimization function, e.g. to MLE or have flow control (if-else).

---

<div class="post-metadata">

**Author:** ![Humphrey\_Lee](https://avatars.discourse-cdn.com/v4/letter/h/a88e4f/32.png) [@Humphrey\_Lee](https://discourse.julialang.org/u/Humphrey_Lee)\
**Post date:** [April 9, 2020, 8:37am UTC](https://discourse.julialang.org/t/estimating-parameters-of-a-function-for-curve-fitting/32404/12 "2020-04-09T08:37:36Z")

</div>

For those interested, below is an example of SSE minimization using solver in Julia. It can be easily modified for the posted question.

```
using Optim, Gadfly, Cairo

# Julia ver. 1.4, Gadfly ver. 1.2.1, Cairo ver. 1.0.2

# y = f(x1, x2) = β_1*x1 + β_2*x2 + β_3 # simply a function
# β_1, β_2 & β_3 are parameters to be solved by the Optim solver
# x1 and x2 are the variables - indepedent variables = inputs
# y is the dependent variable = output

# input data
x1 = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
x2 = [2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0]
#y = 1.0 + 2.0 * x + [-0.3, 0.3, -0.1]
# output data
y = [3.0, 6.0, 11.0, 18.0, 27.0, 38.0, 51.0, 66.0, 83.0, 102.0]

# define function to fit
function fun_y(betas, x1, x2)
	z = betas[1] * x1^2 + betas[2] * x2^2 + betas[3]
	return z
end

# define the sum squared errors
function sqerror(betas)
	err = 0.0
	for i in 1:length(x)
		#pred_i = betas[1] + betas[2] * x[i]^2
		pred_i = fun_y(betas, x1[i], x1[i])
		err += (y[i] - pred_i)^2
	end
	return err
end

# run the solver
res = optimize(sqerror, [0.0, 0.0, 0.0], LBFGS())
println(res)
par = Optim.minimizer(res)

println(par)

# compute the fitted parameters and get the predicted result
pred_i = []
for i in 1:length(y)
	new_i = fun_y(par, x1[i], x2[i]) #par[1] + par[2] * x[i] ^2
	push!(pred_i, new_i)
end

# put into Dataframe for plotting
df = DataFrame()
df[!, :x1] = x1
df[!, :x2] = x2
df[!, :y] = y
df[!, :yy] = pred_i

# plotting
p1 = plot(df, x="x1", y="y",
	Geom.point,
	layer(df, x="x1", y="yy", Geom.line, Theme(default_color=colorant"red")),
	#layer(df, x="x2", y="yy", Geom.line, Theme(default_color=colorant"red")),
	Coord.Cartesian(ymin=0, ymax=110, xmin=0, xmax=12, yflip=false),
)

curr_jl_dir = @ __DIR__
draw(PDF(curr_jl_dir * "\\y_x1.pdf", 5inch, 5inch), p1)

p2 = plot(df, x="x2", y="y",
	Geom.point,
	layer(df, x="x2", y="yy", Geom.line, Theme(default_color=colorant"red")),
	Coord.Cartesian(ymin=0, ymax=110, xmin=0, xmax=12, yflip=false),
)

draw(PDF(curr_jl_dir * "\\y_x2.pdf", 5inch, 5inch), p2)

p3 = plot(df, x="y", y="yy",
	Geom.point,
)

draw(PDF(curr_jl_dir * "\\y_yy.pdf", 5inch, 5inch), p3)

```