kingaa
commented
4 years ago You raise several interesting issues. Among them are:
- What is the relationship between the Euler and the Gillespie simulators?
- What is the relationship between stochastic simulations and a deterministic skeleton?
- What is the role of exact simulations in data analysis?
Perhaps there are others that you would like to discuss as well.
On the first issue: the Gillespie algorithm is an exact simulator for a continuous-time, homogeneous generalized birth-death process. In such a system, stochasticity is purely demographic and per-capita event rates are constant. The Euler method, by contrast, is an approximation scheme, and a first-order one at that. In particular, the Euler process converges to the Gillespie one as the Euler time-step tends to zero. Below, I have modified your code to use a much smaller timestep. You see that the distributions agree to a much higher degree of precision.
On the second issue, there are several reasons to expect that there may be large differences between sample paths of a stochastic process and trajectories of its deterministic skeleton. For one, the deterministic skeleton is not uniquely defined. Second, if the noise is not small, stochastic trajectories may exhibit large deviations from the deterministic trajectories. Third, even if the noise is small, dynamics that are purely transient in the deterministic system will be manifest in the stochastic system due to re-excitation of transients. Given all this, it is reasonable to ask why anyone takes interest in the skeleton at all, and indeed some of the interest is purely historical. However, it is sometimes the case that a deterministic analysis can shed light on the behavior of a stochastic system.
On the third issue, some philosophical reflections are in order. One might expect that it is very important that one have a high-fidelity stochastic simulator—an exact one if possible—in any data analysis. Certainly (it would seem) development and deployment of exact simulation methods would be a useful thing for a mathematically-inclined person to devote time to. But consider the following: How much effort is it worth expending to obtain exact solutions for a model that is itself a poor approximation to the reality one is attempting to describe? How much effort should go into improving the fidelity of the simulations to the model, and how much into making the model/reality connection more faithful?
For example, inasmuch as the SIR model we have been discussing is meant to explain a series of observations from a real system, it is clear that many, potentially important, effects are missing: Are contact rates really constant, or do they perhaps vary over the course of a day? If so, wouldn't it be better to include diurnal variation in the contact rates? Doing so would take Gillespie off the table, unfortunately, although one could develop a custom exact stochastic simulation algorithm that would only necessitate integrating various hazards over short time intervals. One could go on. The point is not that exact simulation methods are without value in a data analysis, but that they are costly, and that one must consider what value one obtains for the price. There can be no single, general answer to the questions I posed in the preceding paragraph, but the issue is present in every data analysis.
To put it another way, in simulation-based analysis, all commerce between model and data passes through the simulator. Thus, in a perfectly exact sense, the simulator is the model that attempts to explain the data. From this point of view, the model equations are merely approximations of the simulator, more or less worthy of study according to how well they represent the simulations.
library(pomp)
library(tidyverse)
## define the steppers
seir_step <-Csnippet("
double rate[3];
double dN[3];
rate[0]=theta * I; // Infection rate
rate[1]=omega; // rate of E to I transitions
rate[2]=mu; // recovery rate
reulermultinom(1,S,&rate[0],dt,&dN[0]); // new infections
reulermultinom(1,E,&rate[1],dt,&dN[1]); // new progressions
reulermultinom(1,I,&rate[2],dt,&dN[2]); // new recoveries
S+=-dN[0]; // update the number of Susceptibles
E+=dN[0]-dN[1]; // update the number of Infections
I+=dN[1]-dN[2]; // update the number of Infections
R+=dN[2]; // update the number of Recoveries
E2I+=dN[1];"
)
seir_gillespie <-
gillespie_hl(
exposure = list("rate = theta * I * S;",
c(S=-1,E=1,I=0,R=0,E2I=0)),
infection = list("rate = omega * E;",
c(S=0,E=-1,I=1 ,R=0,E2I=1)),
recovery = list("rate = mu * I;",
c(S=0,E=0,I=-1,R=1,E2I=0))
)
seir_odes <- Csnippet("
DS = -theta * S * I;
DE = theta * S * I - omega * E;
DI = omega * E - mu * I;
DR = mu * I;
DE2I= omega * E;
")
## parameters
popsize <- 5e4
pars <- c(
theta = 1.3 / popsize / (2.5/7), # basic reproduction number
omega = 1 / (2/7), # latent period duration = 2 days
mu = 1 / (2.5/7), # infectious period duration = 2.5 days
phi = 36, # negative binomial overdispersion
popsize = popsize
)
parnames <- names(pars)
statenames <- c("S", "E", "I", "R", "E2I") #state space variable name
seir <- pomp(
data=data.frame(time=1:30,cases=NA),
times = "time",
t0 = 0, # initial time point
skeleton = vectorfield(seir_odes), # simulates from the latent process
params = pars,
paramnames = parnames,
statenames = statenames,
accumvars = c("E2I"),
rinit = function(popsize, t0, ...) {
return(c(S = popsize-50, E = 25, I = 25, R = 0, E2I = 0))
}
)
## Generate simulations: NOTE WE HAVE TWO SETS OF EULER SIMULATIONS
seir %>%
simulate(
rprocess=euler(seir_step,delta.t=1/7),
statenames=statenames,
paramnames=parnames,
nsim = 1e3,
format="d"
) -> euler_sims1
seir %>%
simulate(
rprocess=euler(seir_step,delta.t=1/700),
statenames=statenames,
paramnames=parnames,
nsim=1e3,
format="d"
) -> euler_sims2
seir %>%
simulate(
rprocess=seir_gillespie,
statenames=statenames,
paramnames=parnames,
nsim=1e3,
format="d"
) -> gillespie_sims
seir %>%
trajectory(format="d") -> ode_path
## summarize and plot ----------------------------------------------------------------
bind_rows(
Euler1 = euler_sims1 %>% select(time,rep=.id,E2I),
Euler2 = euler_sims2 %>% select(time,rep=.id,E2I),
Gillespie = gillespie_sims %>% select(time,rep=.id,E2I),
.id = "simulator"
) %>%
group_by(simulator,time) %>%
summarize(
p=c(0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95),
q=quantile(E2I,prob=p),
label=sprintf("Q%d%%",100*p)
) %>%
select(-p) %>%
spread(label,q) %>%
ungroup() -> quants
## plot
quants %>%
ggplot(aes(x = time, y = `Q50%`)) +
geom_ribbon(aes(ymin = `Q5%`, ymax = `Q95%`, fill=simulator), color=NA, alpha = 0.3) +
geom_ribbon(aes(ymin = `Q10%`, ymax = `Q90%`, fill=simulator), color=NA, alpha = 0.3) +
geom_ribbon(aes(ymin = `Q25%`, ymax = `Q75%`, fill=simulator), color=NA, alpha = 0.3) +
geom_line(aes(color=simulator)) +
geom_line(data = ode_path, aes(x= time, y = E2I), colour = "black") +
guides(color=FALSE)+
labs(y="incidence")+
scale_fill_brewer(
type = "qual", palette = 1,
labels=c(Euler1="Euler, dt=1/7",Euler2="Euler, dt=1/700",Gillespie="Gillespie")) +
scale_color_brewer(type = "qual", palette = 1) +
theme_minimal()
fintzij
I am trying to code up an SEIR model in pomp and am seeing a rather large discrepancy in the distribution of sample paths that I generate when simulating with the Euler-multinomial stepper compared with both gillespie simulation and the deterministic ODE skeleton. I've attached some code below, which is largely based on examples in the pomp documentation. I'm fairly certain that the gillespie/ode algorithms are producing the correct paths as I have also benchmarked the results against another package, but I can't for the life of me figure out where I'm going wrong in the specification of the Euler-multinomial stepper. I'd really appreciate any advice you could offer as I've run out of rocks to look under. It's really a simple model so if it's an error on my end I'm sure it's a rather silly one, but I just can't seem to find it.
Thank you!!!