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

Thomas Petzoldt Thomas.Petzoldt at TU-Dresden.de
Thu Apr 25 23:13:53 CEST 2013


On 4/25/2013 10:29 PM, Andras Farkas wrote:
> 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)



More information about the R-sig-dynamic-models mailing list