[R-sig-dyn-mod] Turbidostat

[Uta Berger]

Hello Thomas,
thank you very much for your quick help.

The second solution (the definition of korr as a second state variable) works also well. However, I do not have a chance to fit korr then using the modFit from the FME package. This problem remains with the approxfun solution, I guess. I probably need a separate model (simple exponential) for the first time steps :-(

I paste my model solution below - just in case.

Best regards,
uta.

library(deSolve)
library(FME)

maxTime <- 15
korr_s <- 0.3  #not fitted but badly assumed - just as an example
signal <- approxfun(1:maxTime, c(rep(korr_s,3),rep(1,10),rep(korr_s,2)))  

turbidostat <- function(time, E0, parms){
    korr <- signal(time)
    with(as.list(c(E0,parms)), {
    dE <- r*korr*E
    list(c(dE))
  })
}

E0 <- c(E=0.1)
parms <- c(r=0.5)
times <- seq(1,15,0.1)

eventdat <- data.frame(var = c("E","E","E","E","E"), time = c(7,9,11,13,15),value=c(rep(0.5,5)),method=c("rep","rep","rep","rep","rep"))


#####observed data
Obs <- (data.frame(time = c(0,1,2,3,7,9,11,13,15),                 
                   E = c(0.1,0.11,0.2,0.3,1.8,1.4,1.45,1.46,1.0)))         

out <- ode(E0,times,turbidostat,parms, events=list(data=eventdat))

plot(Obs$time,Obs$E,pch=16,cex=2.0)
lines(out, type = "l", lwd = 2)

Cost <- function(parms){
                    model <- ode(E0,times,turbidostat,parms, events=list(data=eventdat))
                    modCost(model=model, obs = Obs, x = "time")
                   }
Cost(parms)

Fit <- modFit(f=Cost,p=parms)
deviance(Fit)
Fit$par