Posterior predictive simulation

[fkrauer]

Hi all

I have an ODE model for infectious disease dynamics and I fit the model to a time series of incident deaths using the DEzs sampler. Some parameters of my parameter set theta are estimated, and some are fixed.

I want to simulate model trajectories by sampling n-times from the posterior distribution, simulate a trajectory (i.e. the number of deaths) with the theta from each draw, and calculating the percentiles of all trajectories for the posterior predictive intervals.

I have written a function for the simulation of the trajectories (below). Given that the DEzs runs 3 chains in parallel by default, what is the best strategy for sampling from the posterior? Is it appropriate to sample only the thetas from the first list of the MCMC object, or should I choose any of the three lists randomly?

Thanks for any help.

Cheers from Norway Fabienne

trajsim <- function(model, theta, trace, times, ndraw) {

#Draw n fitted parameter vectors theta from the MCMC object
sample <- getSample(trace, parametersOnly = TRUE, numSamples=ndraw)

traj.rep <- c()
pb <- txtProgressBar(min = 0, max = ndraw, style = 3)

for (i in 1:ndraw) {

#Define the theta as the estimated parameters from the current draw and the fixed parameters
theta.sample <- c(sample[i,], theta[!is.na(theta)])

#Simulate trajectory for selected theta
traj <- suppressWarnings(model(theta.sample, times))
modelpoint <- traj[["Dhinc"]]
sim <- data.frame(cbind(time=traj[["time"]], simobs=NA))

#loop over all time points of the trajectory to generate simulated observations resulting from poisson-error process
for (j in 1:nrow(traj)) {
sim[j,2] <- rpois(n = 1, lambda = modelpoint[j])
}

traj.sim <- merge(traj, sim, by = "time")

# Combine all trajectories in one dataframe
traj.rep <- rbind(traj.rep, cbind(replicate=i, as.data.frame(t(theta.sample)), traj.sim))

#print progress
setTxtProgressBar(pb, i)
}
return(traj.rep)
} 

[fkrauer]

Sorry, I just realised that getSample() returns only one set of parameters per draw anyway, so I don't have to select the chain to sample from.

Issue closed :)

[florianhartig]

yes, exactly!