[R-sig-dyn-mod] FME modMCMC and forcing model parameters

[Andras Farkas]

Dear All,
 
please provide insights to the following: we have:
 
parKM <-structure(list(KM1 = 3.33, KM2 = 3.33), .Names = c("KM1", "KM2"
))
parVC<-structure(list(VC1 = 27.1, VC2 = 27.1), .Names = c("VC1", "VC2"
))
parVP <-structure(list(VP1 = 127, VP2 = 127), .Names = c("VP1", "VP2"
))
parVMAX <-structure(list(VMAX1 = 43.9, VMAX2 = 43.9), .Names = c("VMAX1", 
                                                                 "VMAX2"))
parQ <-structure(list(Q1 = 35.1, Q2 = 35.1), .Names = c("Q1", "Q2"))
time=c(0,2.42, 11.97, 14.17, 23.97)
cp.a=c(0,3.051335071, 0.218260571, 2.99395674, 0.321142542)
observed <-matrix (nc=2,byrow=2,data=rbind(time,cp.a))
colnames(observed) <- c("time", "cp.a")

pars <-pars <-c(parKM,parVC,parVP,parVMAX,parQ)
t <-seq(0,24+24,by=0.01)
intimesout <-c(0, 2, 12, 14, 24, 25, 48)
inputout <-c(225, 0, 225, 0, 300, 0, 0)
forc <- approxfun(intimesout, inputout, method="constant", rule=2)
derivs <- function(t, state, pars) {
  parintimes <- c(0, 12)
  KMinput   <- as.list(parKM)
  KMforc <- approxfun(parintimes, KMinput , method="linear", rule=2)
  VCinput   <- as.list(parVC)
  VCforc <- approxfun(parintimes, VCinput , method="linear", rule=2)
  VPinput   <- as.list(parVP)
  VPforc <- approxfun(parintimes, VPinput , method="linear", rule=2)
  VMAXinput   <- as.list(parVMAX)
  VMAXforc <- approxfun(parintimes, VMAXinput , method="linear", rule=2)
  Qinput   <- as.list(parQ)
  Qforc <- approxfun(parintimes, Qinput , method="linear", rule=2)  
  KM <- KMforc(t)
  VC <- VCforc(t)
  VP <- VPforc(t)
  VMAX <- VMAXforc(t)
  Q <- Qforc(t)
  inp <- forc(t)
  dy1 <- inp - (Q/VC)*state[1]+(Q/VP)*state[2]-((VMAX*state[1])/VC)/(KM+(state[1]/VC))
  dy2 <- (Q/VC) * state[1] - (Q/VP) *state[2]
  dc <- state[1]/VC
  return(list(c(dy1,dy2),cp=dc))
  
}
model <- function(pars,times=seq(0,24+24,by=0.01)) {
  state <- c(a=0,b=0)
  return(ode(y = state, times = times, func = derivs, parms = pars,method="lsoda",rtol = 1e-8,atol =1e-8,hmax=0.01,maxsteps=50000))
}
tout <-c(0, 2, 2.42, 11.97, 12, 14, 14.17, 23.97, 24, 25, 48)
 
obj <- function(x, parset = names(x)) {
  pars[parset] <- x
  out <- model(pars, tout)
  return(modCost(obs = observed, model = out))
}

final <- modMCMC(p = c(unlist(parKM),unlist(parVC),unlist(parVP),unlist(parVMAX),unlist(parQ)), f = obj, prior=NULL,lower = rep(0.0000001, length(pars)),jump = NULL, niter =50,updatecov=10, wvar0 = 1,var0=NULL,burninlength = 5)
 
 
my goal here is to estimate best fit parameters for KM1-2,VC1-2,VP1-2,VMAX1-2, and Q1-2 usong the model and observed data... #1 for all parameters is forced in at time 0 versus #2 at time 12. The transition from #1 to #2 should be linear over time as per the forcing functions above... Once the code runs (takes a while) the final$bestpar values remain the same as the starting point. Any thoughts as to why it may be that the fitting of the model to the observed data is not working in this case?
 
appreciate your input,
 
Andras 

	[[alternative HTML version deleted]]