ODE-based nonlinear least square in R

214 Views Asked by At

I have to set a nonlinear least square regressor on a function that comes from the integration of an ode. Something like

#rm(list=ls())
library(deSolve)

set.seed(1) # for reproducibility
t = seq(0,2*pi,by=2*pi/50) # sampling times
alpha = 1; beta = 1; sigma = 0.05 # test parameters
x = beta/alpha*(1-exp(-alpha*t))+rnorm(length(t),sd=sigma) # signal

D = data.frame(t=t,x=x)
plot(x~t,D,xlim=c(0,2*pi),ylim=c(0,1.2))

dxdt = function(t,x,p){ # differential equation
  alpha = p[1]
  beta = p[2]
  list(-alpha*x+beta)
}
f = function(t,p){ # de integrated (times=c(0,t) to ensure first time=0)
    ode(y=c(x=0),func=dxdt,times=c(0,t),parms=p)[-1,]
}

sol=f(t,c(1,1)) # solution
lines(sol) # plot the solution

argminf = function(data){ # nonlinear regressor
   nls(x~f(t,p),data=data,start=list(p=c(1,1)))
}

hatf = argminf(D)

in calling argminf() i get back

Error in qr.qty(QR, resid) : 
'qr' and 'y' must have the same number of rows

what is my misunderstanding in setting this code?

1

There are 1 best solutions below

0
On

You should apply the following corrections

f = function(t, p) { # de integrated (times = c(0, t) to ensure first time = 0)
  ode(y = c(x = 0), func = dxdt, times = c(0, t), parms = p
  )[-1, 2]
}
sol = f(t, c(1, 1)) # solution
lines(t, sol) # plot the solution

argminf = function(data) { # nonlinear regressor
  nls(x ~ f(t, p), data = data, start = list(p = c(1, 1)))
}

hatf = argminf(D)
lines(t, fitted(hatf))