Estimating parameters of a function for curve fitting

Specific Domains Optimization (Mathematical)
1

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.

######
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?

2

You might try Optim.jl.

3

there is LsqFit.jl for just fitting of non linear functions to data.

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)
4

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

# 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

TypeError: non-boolean (Missing) used in boolean context
5

add this:

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
6

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.

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
7

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.

8

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.

9

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

10

What’s wrong with the examples here?

11

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).

12

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)