Hi Himanshu,
the use of optim is described on its help page.
In addition to this package FME provides additional functionality for
fitting ODE models. See FME help files, package vignettes and the
example below for details.
Hope it helps
Thomas Petzoldt
################################################################
library(deSolve)
library(FME)
## function derivs
derivs <- function(time, y, pars) {
with (as.list(c(y, pars)), {
dy = a1 * Y * (1 - Y/Ymax) - a2 * (Y + 0.001)
dz = a3 * Y - a4 * Z;
return(list(c(dy, dz)))
})
}
## parameters
pars <- c(a1 = 0.9, a2 = 0.7, a3 = 0.06, a4 = 0.02)
Ymax <- c(0.8)
## initial values
y <- c(Y = 0.2, Z = 0.1)
## time steps
times <- c(seq(0, 10, 0.1))
## numerical solution of ODE
out <- ode(y = y, parms = pars, times = times, func = derivs)
## example observation data
yobs <- data.frame(
time = 0:7,
Y = c(0.2, 0.195, 0.19, 0.187, 0.185, 0.183, 0.182, 0.18),
Z = c(0.1, 0.11, 0.119, 0.128, 0.136, 0.145, 0.153, 0.16)
)
## 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(a1 = 1, a2 = 0.6, a3 = 0.1, a4 = 0.1)
## fit the model; nprint = 1 shows intermediate results
fit <- modFit(f = modelCost, p = startpars, control = list(nprint = 1))
## Note the high correlation between a1 and a2, resp a3 and a4
## that can make parameter identification difficult
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)
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
