[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