Hello list:

I am setting up a textbook example of a Verlhust population growth model in class (see code).  I have actual census data as indicated in the code and i would like to fit the "alfa" and "L" parameters using fitOdeModel.   The problem is that the fit routine seems to be leaving the L parameter unchanged: (results bellow)

> resultado
$par
       alfa           L 
3.65782e-02 1.00000e+08 

$value
[1] 0.002420281

$counts
function gradient 
      53        4 

$convergence
[1] 0

$message
NULL

I suspect this has to do with the initial values of "alfa" and "L", which differ in several orders of magnitude (alfa is around 1 to 5 percent, L is probably in the order of tens of millions).  I know there is a "scalepar" argument to fitOdeModel, but I've not managed to get this working.  Any suggestions?

Thanks in advance and regards,
jose romero

The code:

library("simecol")
#esta es la data de la poblacion venezolana según censo
censo <- data.frame(ano=c(1936,1950,1961,1971,1981),
           poblacion=c(3194543,5043838,
           7523999,10721522,14602480))
#se define el modelo de verlhust
verlhust <- new("odeModel",
            main=function (time,init,parms) {
                x <- init
                alfa <- parms["alfa"]
                L <- parms["L"]
                dx <- alfa*x*(1-x/L)
                list(dx)
            },
            parms=c(alfa=0.04,L=100000000),
            times=seq(from=1936,to=2010,by=0.5),
            init=c(poblacion=3194543),
            solver=rk4)
vsim1 <- out(sim(verlhust))
sim1 <- subset(vsim1,time%in%censo$ano)
#ahora vamos a ajustar los parametros de acuerdo
#a la data empirica.
resultado <- fitOdeModel(verlhust,
                         obstime=censo$ano,
                         yobs=subset(censo,select=poblacion),
                         fn=ssqOdeModel,
                         initialize=FALSE,
                         method="BFGS",
                         lower=c(alfa=0,L=0))
#actualizamos el valor del parámetro alfa y volvemos a simular
resultado
parms(verlhust) <- resultado$par
vsim2 <- out(sim(verlhust))
sim2 <- subset(vsim2,time%in%censo$ano)
#ahora comparamos
sim1
cat("Error sim1: ",sum(abs(sim1$poblacion-censo$poblacion)),"\n")
cat("Proyeccion 2010 : ",vsim1$poblacion[vsim1$time==2010],"\n")
sim2
cat("Error sim2: ",sum(abs(sim2$poblacion-censo$poblacion)),"\n")
cat("Proyeccion 2010 : ",vsim2$poblacion[vsim2$time==2010],"\n")


	[[alternative HTML version deleted]]

