[R-sig-dyn-mod] Failed attempt to use modFit- part 2
Thomas Petzoldt
Thomas.Petzoldt at TU-Dresden.de
Sat Dec 29 18:28:32 CET 2012
Hi,
here another ODE model that fits your data rather well. The approach
uses deSolve without simecol. It applies a two-equation growth model
that splits the population in a non-reproducing and a an actively
reproducing fraction.
Note also that the parameters are again highly correlated and that
versions with fixed parameter "a" would also result in a plausible
curves, but with very different r and K.
Thomas P.
##===================================================================
library(FME)
library(deSolve)
## model function
model <- function(t, y, p) {
with(as.list(c(y, p)), {
dy1 <- -a * y1 # "inactive" part of population
dy2 <- a * y1 + r * y2 * (1 - y2 / K) # growing population
list(c(dy1, dy2), pob = y1 + y2) # total pop = y1 + y2
})
}
## cost function
fCost <- function(p) {
parms[names(p)] <- p # "names" allow to fit parameter subsets
out <- ode(init, times = censo$time, model, parms = parms)
cost <- modCost(out, censo)
cat(parms, " : ", cost$model, "\n")
cost
}
## data set
censo <- read.table("censos.csv", header = TRUE)
init <- c(y1 = 1732411, y2 = 0)
parms <- c(r = 0.05, K = 3.0e7, a = 0.001)
times <- 1873:2100
out <- ode(init, times, model, parms)
plot(censo, xlim = c(1850, 2100), ylim = c(0, 5e7))
lines(out[, 1], out[, 4]) ## initial parameters
##--------------------------------------------------------------------
## fit model (all parameters)
result <- modFit(fCost, p = parms[c("r", "K", "a")])
## alternative: leave "a" untouched
#result <- modFit(fCost, p = parms[c("r", "K")])
## retrieve and display parameters
(parms[names(coef(result))] <- coef(result))
deviance(result)
## simulate fitted model and plot results
out <- ode(init, times = 1873:2100, model, parms)
lines(out[, 1], out[, 4], col = "red")
legend("topright", lty = 1, col = c("black", "red"),
legend = c("initial parms", "fitted"))
-------------- next part --------------
A non-text attachment was scrubbed...
Name: censos.csv
Type: application/vnd.ms-excel
Size: 259 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-dynamic-models/attachments/20121229/4b50b2f8/attachment.xlb>
More information about the R-sig-dynamic-models
mailing list