Thank you Dr Petzoldt, apreciate the help. Revising my parameters as you have suggested helped better understand the question at hand... Working with the modMCMC function have raised another question. I got the system to run and give me answers for the most part, but not always do I get one that I expect (see code below).... Also, the answeres for parameters that I get are not different enough when I change the independent variables to a degree that should definitelly result in a significant magnitude of change in the estimated parameters. I would apreciate the groups help to point out what I could do better. here is the code:
 
require(deSolve)
require(FME)
pkmod <- function(t, y, p) {
  inp <- forc(t)
  dy1 <- - p["k"]* y[1] + inp
  list(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 <- c(0,13,18)
yini <- c(Central = 0)
p    <- c(k=0.07)
result <- lsoda(yini, times, pkmod, p, rtol = 1e-08, atol = 1e-08)
result <- data.frame(result)
names(result) <- c("time","Concentration")
plot(result$Concentration ~ result$time, type="b", main="Central compartment", xlab="time", ylab="Amount")
modelout <-cbind(result[1],result[2])
time <-c(0,13,18)
Concentration <-c(0,628.3,514.4)
obs <-cbind(time,Concentration)
cost <-function(p){
  
  modCost(model=modelout,obs=obs,x="time",y=NULL,err=NULL,cost=NULL,weight="none",scaleVar = FALSE)
}
 
prior <- function(p)
  -2*sum(log(dnorm(p, mean = 0.07, sd = 0.03)))
  
final <-modMCMC(f=cost, p=p, jump = NULL, lower=c(0.00000001), prior = prior,niter = 5000,  burninlength = 200)
summary(final)
plot(final)

If we run this code, then the # of accepted runs is 0, which I found suprising and makes me believe that the code is off somewhere... The model output to be used in modCost was calculated with a value of 0.04 (the value I am expecting to be returned), which is 1 SD from the prior mean, so I would believe that the system can find it. I must say that I am not sure about the prior function being written accuratelly, although the returned mean and sd parameters although are not what I am looking for, they do make sense from the distributions perspective in summary(final). Appreciate your insights.
 
thanks,
 
Andras
--- On Thu, 4/25/13, Thomas Petzoldt <Thomas.Petzoldt@TU-Dresden.de> wrote:


From: Thomas Petzoldt <Thomas.Petzoldt@TU-Dresden.de>
Subject: Re: [R-sig-dyn-mod] ] Bayesian model fitting and deSolve
To: "Special Interest Group for Dynamic Simulation Models in R" <r-sig-dynamic-models@r-project.org>
Date: Thursday, April 25, 2013, 5:13 PM


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)

_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

	[[alternative HTML version deleted]]

