Dear List,
 
wonder if you could provide some input in the following (I am doing something wrong, but can not tell what it is...)
 
the equations here describe a one compartment intravenous infusion model with analytical solutions (ifelse to distinguish during and post infusion calculations). The problem I am facing is that for some reason the parameters (pars) are not getting estimated by modMCMC as I would like it to (the problem is with my code I am sure somewhere, NOT with modMCMC).... My calc function and model do produce what I need (see plot(model(pars),type="l")), I am estimating the cost, but then get stuck on modMCMC. While it will do the runs, the bestpar will be equal to my prior information. Your input is greatly appreciated,
 
thanks,
 
Andras
 
 
 
require(deSolve)
require(FME)
tin <-c(0.5,0.5)
D <-c(400,400)
dtimes <-c(0,12,24)
pars <-list(k = 0.2,v=30)
calc <- function(t) {
  dc1 <-ifelse(t<=tin[1],((D[1]/tin[1])*(1/(pars$k*pars$v)))*(1-(exp(-pars$k*t))),((D[1]/tin[1])*(1/(pars$k*pars$v)))*(1-(exp(-pars$k*tin[1])))*exp(-pars$k*(t-tin[1])) ) 
  return(c(dc1))
}
model <- function(pars,t=seq(0, dtimes[2], by = 0.1)) {
  return(data.frame(time=t,y=calc(t)))
}
plot(model(pars),type="l")
observed <- matrix (nc=2,byrow=2,data=c(0.8,15,8,8))
colnames(observed) <- c("time", "y")
tout <-  c(observed[,"time"])
obj <- function(x, parset = names(x)) {
  pars[parset] <- x
  out <- model(pars, tout)
  return(modCost(obs = observed, model = out))
}
prior <- function(p)
  return( sum(((p-c(unlist(unname(pars))))/c(unlist(unname(pars))*0.4))^2 ))
final <- modMCMC(p = unlist(pars), f = obj, prior=prior,lower=c(0.0001,0.1),jump = NULL, niter = 10000,updatecov=10, wvar0 = 1,var0=NULL,burninlength = 100)
 
best1 <-as.list(final$bestpar)
out1 <- model(pars)
plot(out1,type="l",ylim=c(0,30))
pars <-as.list(final$bestpar)
out2 <- model(pars)
lines(out2,col="red")
points(observed)
	[[alternative HTML version deleted]]

