CODEWITHSUNDEEP
roptimizationode
P.Galanis
P.Galanis

Reputation: 1

Optimise an ODE and recalculate the initial value in R

I was wondering if I could have some input in the following problem: I try to have an optimal value returned for my function, after using modCost and modFit to get results.

So what essentially I try to do is toy around with a basic forecast, optimizing to get a good fit (which I think works). But I would also like the 'yi' parameter to be included as a parameter in the Modfit function.

Is there a specific function I can use for that perhaps or a workaround? Any help/suggestion would be greatly appreciated.

Apologies for the poorly written code, first post etc :)

library(deSolve) 
library(FME) 

## function derivs 
derivs <- function(time, y, pars) { 
   with (as.list(c(y, pars)), { 
     dy = a * y   
     return(list(c(dy))) 
   }) 
} 

## parameters 
pars <- c(a = 0.1) 

yi=6

## initial values 
y <- c(Y = yi) 

## time steps 
times <- c(seq(2005, 2017, 1)) 

## numerical solution of ODE 
out <- ode(y = y, parms = pars, times = times, func = derivs) 

## example observation data 
yobs <- data.frame( 
   time = 2005:2017, 
   Y = c(3,6,9,10,12,8,7,9,14,16,18,15,19) 
) 

## model cost function, see help file and vignette for details 
modelCost <- function(p) { 
   out <- ode(y = y, func = derivs, parms = p, times = yobs$time) 
   return(modCost(out, yobs)) 
} 

## start values for the parameters 
startpars <- c(a = 0.1) 

## fit the model; nprint = 1 shows intermediate results 
fit <- modFit(f = modelCost, p = startpars, control = list(nprint = 1)) 


summary(fit) 

## graphical result 
out2 <- ode(y = y, parms = startpars, times = times, func = derivs) 
out3 <- ode(y = y, parms = fit$par, times = times, func = derivs) 
plot(out, out2, out3, obs = yobs) 

legend("topleft", legend=c("original", "startpars", "fitted"), 
   col = 1:3, lty = 1:3)

Upvotes: 0

Views: 432

Answers (1)

Dan
Dan

Reputation: 12074

Adapting only the following lines of code seems to work:

## model cost function, see help file and vignette for details 
modelCost <- function(p) { 
  out <- ode(y = p["Y"], func = derivs, parms = p["a"], times = yobs$time) 
  return(modCost(out, yobs)) 
} 

## start values for the parameters 
startpars <- c(Y = 6, a = 0.1)

You'll notice in the solution that both the initial condition and parameter have been optimised for:

#$par
#         Y          a 
#5.81276026 0.09872004

Note that the differential equation has an analytical solution: y = yi * exp(a * (time - time_0)). Consequently, a simpler approach is to fit this directly to the data using optim:

# My cost function
ls_cost <- function(p){
  sum((p["Y"] * exp(p["a"] * (yobs$time - min(yobs$time))) - yobs$Y)^2)
}

# Optimise for initial condition and 'a'
optim(par = startpars, fn = ls_cost)

which gives,

#$par
#         Y          a 
#5.81256899 0.09872287

Pretty close to the other approach.

Now, comparing the two using microbenchmark we get the following:

Unit: milliseconds
                                 expr      min       lq     mean  median       uq      max neval
 modFit(f = modelCost, p = startpars) 36.22868 37.52432 40.68432 38.3405 40.02509 85.91644   100

Unit: milliseconds
                                 expr      min       lq     mean   median       uq      max neval
 optim(par = startpars, fn = ls_cost) 1.927786 1.980567 2.082507 2.010147 2.094988 5.522633   100

So, the latter approach is an order of magnitude faster.

Upvotes: 1

Related Questions