Dear All,
please help with some thoughts on overcoming the following issues, if possible:
#R Code
require(deSolve)
require(FME)
pars <- list(k = 0.06,v=18)
intimes <- c(0,0.5,12,12.5,50)
input <- c(800,0,800,0,0)
forc <- approxfun(intimes, input, method="constant")
model <- function(pars, times=seq(0, 50, by = 1)) {
derivs <- function(t, state, pars) {
with(as.list(c(state, pars)), {
inp <- forc(t)
dy1 <- - k* y[1] + inp/v
return(list(c(dy1)))
})
}
state <- c(y = 0)
return(lsoda(y = state, times = times, func = derivs, parms = pars))
}
observed <- matrix
(nc=2,byrow=2,data=c(0,0,0.5,19.7,3,17.4,5.5,15.3,8,13.5,11.5,11.3,13,29.8,15.5,26.3,16.5,25,18,23.2,24,17.2,30,12.7,40,7.7,49.5,4.8))
colnames(observed) <- c("time", "y")
obj <- function(x, parset = names(x)) {
pars[parset] <- x
tout <- seq(0, 50, by = 1)
out <- model(pars, tout)
return(modCost(obs = observed, model = out))
}
prior <- function(p)
return( sum(((p-c(0.06,18))/c(0.04,6))^2 ))
final <- modMCMC(p = c(k=0.06,v=18), f = obj,
prior=prior,lower=c(0.0001,0.1),jump = NULL, niter = 1000,updatecov=100, wvar0
= 1,var0=NULL,burninlength = 50)
summary(final)
plot(final)
the values for "observed" were generated with k=0.05, v=20, which are the
retunrs I expect from the modMCMC run. Please comment on:
1. is there any way to get modMCMC to run multiple times faster with this data?
2. my accepted number of runs is very low. Is there a way to increase that?
3. any thoughts on improving convergence?
4. at times I only have one or two "observed" data points to work with, which
often results in "suboptimal" outputs from the modMCMC runs. any thoughts on
how to improve identifiability of parameter values with less "observed" data
points being available?
thank you for the help,
Andras
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
