[R-sig-dyn-mod] ] Bayesian model fitting and deSolve

Thank you Dr Petzoldt, I did review the examples and I think I am
very close to solving the problem! I decided to go with using the
modCost function to generate my f for modMCMC, and am getting
reasonable results. I will spend some more time on it, still having
some issues with defining the function for  a normal prior...
Besides, that I would like to ask if the group could help me
identifying why the differential equation system does not pick up my
second input of rate 800 for 0.5 h starting at 12 h and ending at
12.5? I have gone through it several times, but can not find the
error. I appreciate the help.

[...]

I cannot see an "error". As you write yourself, it is an input *rate* so 
an input of 800 means 800 per time unit, i.e. 400 in the interval from 
12 and 12.5.

The second factor is your relatively high value of "k, so parts of y1 
are "decomposed" before the input is ready.

Therefore, you cannot see the "800" in the figure. Please do not mix the 
concepts of events (immediate change of states) with forcing (or input), 
that is a rate.


Fort a better understanding, reduce the time step (for better resolution 
of graphics) and the k parameter.

Hope it helps

ThPe




require(deSolve)

pkmod <- function(t, y, p) {
   inp <- forc(t)
   dy1 <- - p["k"]* y[1] + inp
   list(c(dy1), inp=inp)
}
intimes <- c(0,0.5,12,12.5,24)
input   <- c(800,0,800,0,0)
forc <- approxfun(intimes, input, method="constant")
times <- seq(0,24,by=0.01)
yini <- c(Central = 0)
#p    <- c(k=0.108)
p    <- c(k=0.001)
result <- ode(yini, times, pkmod, p, rtol = 1e-8, atol = 1e-8)

plot(result)