1.2 What’s the proper way to cite pomp?
@@ -214,166 +297,166 @@3 Implementing POMP mo
3.1 How can I include a vector of variables in a C snippet?
A C snippet can make use of any feature of the C language. In particular, we can use pointers to give access to arrays. For example, consider the following, which implement a simple chain of random variables.
-rSim <- Csnippet("
- double *x = &X1;
- int i, n = (int) N;
- for (i = 0; i < n-1; i++) x[i] = x[i+1];
- x[n-1] = rnorm(0,1);
-")
-
-rInit <- Csnippet("
- double *x = &X1;
- int i, n = (int) N;
- for (i=0; i < n; i++) x[i] = 0.0;
-")
-
-rMeas <- Csnippet("
- Y1 = rnorm(X2,2);
-")
+rSim <- Csnippet("
+ double *x = &X1;
+ int i, n = (int) N;
+ for (i = 0; i < n-1; i++) x[i] = x[i+1];
+ x[n-1] = rnorm(0,1);
+")
+
+rInit <- Csnippet("
+ double *x = &X1;
+ int i, n = (int) N;
+ for (i=0; i < n; i++) x[i] = 0.0;
+")
+
+rMeas <- Csnippet("
+ Y1 = rnorm(X2,2);
+")
The following simulates and plots using the above.
-library(ggplot2)
-library(tidyr)
-
-simulate(
- times=1:20,t0=0,
- rprocess=euler(rSim,delta.t=1),
- rmeasure=rMeas,
- rinit=rInit,
- obsnames="Y1",
- statenames=sprintf("X%d",1:5),
- paramnames="N",
- params=c(N=5),
- format="data.frame"
-) |>
- pivot_longer(-c(.id,time)) |>
- ggplot(aes(x=time,y=value,group=name,color=name))+
- geom_line()+
- theme_bw()
+library(ggplot2)
+library(tidyr)
+
+simulate(
+ times=1:20,t0=0,
+ rprocess=euler(rSim,delta.t=1),
+ rmeasure=rMeas,
+ rinit=rInit,
+ obsnames="Y1",
+ statenames=sprintf("X%d",1:5),
+ paramnames="N",
+ params=c(N=5),
+ format="data.frame"
+) |>
+ pivot_longer(-c(.id,time)) |>
+ ggplot(aes(x=time,y=value,group=name,color=name))+
+ geom_line()+
+ theme_bw()
3.2 Can I write a pomp for multivariate data?
Yes, this is no problem. The data you supply to pomp can contain multiple variables. You simply refer to these variables by name in writing the rmeasure and dmeasure C snippets. The following example illustrates this for a multivariate Ornstein-Uhlenbeck process (as in pompExample("ou2")) with a slightly peculiar measurement model.
-simulate(
- times=1:100, t0=0,
- rmeasure=Csnippet("
- y1 = rnorm(x1+x2,sigma1);
- y2 = rnorm(x2-x1,sigma2);
- "),
- dmeasure=Csnippet("
- lik = dnorm(y1,x1+x2,sigma1,1)+dnorm(y2,x2-x1,sigma2,1);
- lik = (give_log) ? lik : exp(lik);
- "),
- rprocess=discrete_time(
- step.fun=Csnippet("
- double tx1, tx2;
- tx1 = rnorm(a11*x1 + a12*x2,nu1);
- tx2 = rnorm(a21*x1 + a22*x2,nu2);
- x1 = tx1; x2 = tx2;
- "),
- delta.t=1),
- rinit=Csnippet("
- x1 = 0;
- x2 = 0;
- "),
- obsnames=c("y1","y2"),
- statenames=c("x1","x2"),
- paramnames=c("a11","a12","a21","a22","sigma1","sigma2","nu1","nu2"),
- params=c(a11=0.5,a12=-0.1,a21=0.2,a22=-1,nu1=0.3,nu2=0.1,sigma1=0.1,sigma2=0.3)
-) -> obj
+simulate(
+ times=1:100, t0=0,
+ rmeasure=Csnippet("
+ y1 = rnorm(x1+x2,sigma1);
+ y2 = rnorm(x2-x1,sigma2);
+ "),
+ dmeasure=Csnippet("
+ lik = dnorm(y1,x1+x2,sigma1,1)+dnorm(y2,x2-x1,sigma2,1);
+ lik = (give_log) ? lik : exp(lik);
+ "),
+ rprocess=discrete_time(
+ step.fun=Csnippet("
+ double tx1, tx2;
+ tx1 = rnorm(a11*x1 + a12*x2,nu1);
+ tx2 = rnorm(a21*x1 + a22*x2,nu2);
+ x1 = tx1; x2 = tx2;
+ "),
+ delta.t=1),
+ rinit=Csnippet("
+ x1 = 0;
+ x2 = 0;
+ "),
+ obsnames=c("y1","y2"),
+ statenames=c("x1","x2"),
+ paramnames=c("a11","a12","a21","a22","sigma1","sigma2","nu1","nu2"),
+ params=c(a11=0.5,a12=-0.1,a21=0.2,a22=-1,nu1=0.3,nu2=0.1,sigma1=0.1,sigma2=0.3)
+) -> obj
3.3 How do I deal with missing data?
Missing data are a common complication. pomp can handle NAs in the data, but if it is needed, the measurement model probability density function, dmeasure must be written so as to deal with NAs appropriately. For example, look at the following variant of the SIR model describing the influenza outbreak in a boarding school:
-library(ggplot2)
-
-read.csv(text="
-B,day
-NA,0
-1,1
-6,2
-26,3
-73,4
-NA,5
-293,6
-258,7
-236,8
-191,9
-124,10
-69,11
-26,12
-11,13
-4,14
-") -> dat
-
-dat |>
- ggplot(aes(x=day, y=B)) +
- geom_line() + geom_point()
+library(ggplot2)
+
+read.csv(text="
+B,day
+NA,0
+1,1
+6,2
+26,3
+73,4
+NA,5
+293,6
+258,7
+236,8
+191,9
+124,10
+69,11
+26,12
+11,13
+4,14
+") -> dat
+
+dat |>
+ ggplot(aes(x=day, y=B)) +
+ geom_line() + geom_point()
Data are missing at days 0 and 5.
We create a pomp object implementing a simple SIR process model and a binomial measurement model, as in the original example. The only difference is in the dmeasure:
-library(pomp)
-
-sir_step <- Csnippet("
- double dN_SI = rbinom(S,1-exp(-Beta*I/N*dt));
- double dN_IR = rbinom(I,1-exp(-gamma*dt));
- S -= dN_SI;
- I += dN_SI - dN_IR;
- R += dN_IR;
- H += dN_IR;
-")
-
-sir_init <- Csnippet("
- S = N-1;
- I = 1;
- R = 0;
- H = 0;
-")
-
-rmeas <- Csnippet("B = rbinom(H,rho);")
-
-dmeas <- Csnippet("
- if (ISNA(B)) {
- lik = (give_log) ? 0 : 1;
- } else {
- lik = dbinom(B,H,rho,give_log);
- }
-")
-
-dat |>
- pomp(time="day",t0=0,
- rprocess=euler(sir_step,delta.t=1/12),
- rinit=sir_init,
- rmeasure=rmeas,
- dmeasure=dmeas,
- accumvars="H",
- paramnames=c("Beta","gamma","N", "rho"),
- statenames=c("S","I","R","H")
- ) -> sir_na
+library(pomp)
+
+sir_step <- Csnippet("
+ double dN_SI = rbinom(S,1-exp(-Beta*I/N*dt));
+ double dN_IR = rbinom(I,1-exp(-gamma*dt));
+ S -= dN_SI;
+ I += dN_SI - dN_IR;
+ R += dN_IR;
+ H += dN_IR;
+")
+
+sir_init <- Csnippet("
+ S = N-1;
+ I = 1;
+ R = 0;
+ H = 0;
+")
+
+rmeas <- Csnippet("B = rbinom(H,rho);")
+
+dmeas <- Csnippet("
+ if (ISNA(B)) {
+ lik = (give_log) ? 0 : 1;
+ } else {
+ lik = dbinom(B,H,rho,give_log);
+ }
+")
+
+dat |>
+ pomp(time="day",t0=0,
+ rprocess=euler(sir_step,delta.t=1/12),
+ rinit=sir_init,
+ rmeasure=rmeas,
+ dmeasure=dmeas,
+ accumvars="H",
+ paramnames=c("Beta","gamma","N", "rho"),
+ statenames=c("S","I","R","H")
+ ) -> sir_na
In the above, to check for missing data, we use the ISNA macro, which is part of R’s C API. Note that the dmeasure returns a likelihood of 1 (log likelihood 0) for the missing data. [What’s the probability of seeing something if you don’t look?] When there is an observation, it returns a binomial likelihood as usual.
Our simulations now fill in simulated values for the missing data,
-sir_na |>
- simulate(
- params=c(Beta=3,gamma=2,rho=0.9,N=2600),
- nsim=10,
- format="data.frame",
- include.data=TRUE
- ) |>
- ggplot(aes(x=day,y=B,group=.id,color=.id=="data"))+
- geom_line()+
- guides(color="none")+
- theme_bw()
+sir_na |>
+ simulate(
+ params=c(Beta=3,gamma=2,rho=0.9,N=2600),
+ nsim=10,
+ format="data.frame",
+ include.data=TRUE
+ ) |>
+ ggplot(aes(x=day,y=B,group=.id,color=.id=="data"))+
+ geom_line()+
+ guides(color="none")+
+ theme_bw()
and the particle filter handles the missing data correctly:
-sir_na |>
- pfilter(Np=1000,params=c(Beta=3,gamma=2,rho=0.9,N=2600)) |>
- as.data.frame() |>
- pivot_longer(-day) |>
- ggplot(aes(x=day,y=value,color=name))+
- guides(color="none")+
- geom_line()+
- facet_wrap(~name,ncol=1,scales='free_y')+
- theme_bw()
+sir_na |>
+ pfilter(Np=1000,params=c(Beta=3,gamma=2,rho=0.9,N=2600)) |>
+ as.data.frame() |>
+ pivot_longer(-day) |>
+ ggplot(aes(x=day,y=value,color=name))+
+ guides(color="none")+
+ geom_line()+
+ facet_wrap(~name,ncol=1,scales='free_y')+
+ theme_bw()
In the above particle filter computation, notice that the effective sample size (ESS) is full, as it should be, when the missing data contribute no information.
@@ -382,103 +465,103 @@ 3.4 I have t0 much l
3.1 How can I include a vector of variables in a C snippet?
A C snippet can make use of any feature of the C language. In particular, we can use pointers to give access to arrays. For example, consider the following, which implement a simple chain of random variables.
-rSim <- Csnippet("
- double *x = &X1;
- int i, n = (int) N;
- for (i = 0; i < n-1; i++) x[i] = x[i+1];
- x[n-1] = rnorm(0,1);
-")
-
-rInit <- Csnippet("
- double *x = &X1;
- int i, n = (int) N;
- for (i=0; i < n; i++) x[i] = 0.0;
-")
-
-rMeas <- Csnippet("
- Y1 = rnorm(X2,2);
-")rSim <- Csnippet("
+ double *x = &X1;
+ int i, n = (int) N;
+ for (i = 0; i < n-1; i++) x[i] = x[i+1];
+ x[n-1] = rnorm(0,1);
+")
+
+rInit <- Csnippet("
+ double *x = &X1;
+ int i, n = (int) N;
+ for (i=0; i < n; i++) x[i] = 0.0;
+")
+
+rMeas <- Csnippet("
+ Y1 = rnorm(X2,2);
+")The following simulates and plots using the above.
-library(ggplot2)
-library(tidyr)
-
-simulate(
- times=1:20,t0=0,
- rprocess=euler(rSim,delta.t=1),
- rmeasure=rMeas,
- rinit=rInit,
- obsnames="Y1",
- statenames=sprintf("X%d",1:5),
- paramnames="N",
- params=c(N=5),
- format="data.frame"
-) |>
- pivot_longer(-c(.id,time)) |>
- ggplot(aes(x=time,y=value,group=name,color=name))+
- geom_line()+
- theme_bw()library(ggplot2)
+library(tidyr)
+
+simulate(
+ times=1:20,t0=0,
+ rprocess=euler(rSim,delta.t=1),
+ rmeasure=rMeas,
+ rinit=rInit,
+ obsnames="Y1",
+ statenames=sprintf("X%d",1:5),
+ paramnames="N",
+ params=c(N=5),
+ format="data.frame"
+) |>
+ pivot_longer(-c(.id,time)) |>
+ ggplot(aes(x=time,y=value,group=name,color=name))+
+ geom_line()+
+ theme_bw()
3.2 Can I write a pomp for multivariate data?
Yes, this is no problem. The data you supply to pomp can contain multiple variables. You simply refer to these variables by name in writing the rmeasure and dmeasure C snippets. The following example illustrates this for a multivariate Ornstein-Uhlenbeck process (as in pompExample("ou2")) with a slightly peculiar measurement model.
simulate(
- times=1:100, t0=0,
- rmeasure=Csnippet("
- y1 = rnorm(x1+x2,sigma1);
- y2 = rnorm(x2-x1,sigma2);
- "),
- dmeasure=Csnippet("
- lik = dnorm(y1,x1+x2,sigma1,1)+dnorm(y2,x2-x1,sigma2,1);
- lik = (give_log) ? lik : exp(lik);
- "),
- rprocess=discrete_time(
- step.fun=Csnippet("
- double tx1, tx2;
- tx1 = rnorm(a11*x1 + a12*x2,nu1);
- tx2 = rnorm(a21*x1 + a22*x2,nu2);
- x1 = tx1; x2 = tx2;
- "),
- delta.t=1),
- rinit=Csnippet("
- x1 = 0;
- x2 = 0;
- "),
- obsnames=c("y1","y2"),
- statenames=c("x1","x2"),
- paramnames=c("a11","a12","a21","a22","sigma1","sigma2","nu1","nu2"),
- params=c(a11=0.5,a12=-0.1,a21=0.2,a22=-1,nu1=0.3,nu2=0.1,sigma1=0.1,sigma2=0.3)
-) -> objsimulate(
+ times=1:100, t0=0,
+ rmeasure=Csnippet("
+ y1 = rnorm(x1+x2,sigma1);
+ y2 = rnorm(x2-x1,sigma2);
+ "),
+ dmeasure=Csnippet("
+ lik = dnorm(y1,x1+x2,sigma1,1)+dnorm(y2,x2-x1,sigma2,1);
+ lik = (give_log) ? lik : exp(lik);
+ "),
+ rprocess=discrete_time(
+ step.fun=Csnippet("
+ double tx1, tx2;
+ tx1 = rnorm(a11*x1 + a12*x2,nu1);
+ tx2 = rnorm(a21*x1 + a22*x2,nu2);
+ x1 = tx1; x2 = tx2;
+ "),
+ delta.t=1),
+ rinit=Csnippet("
+ x1 = 0;
+ x2 = 0;
+ "),
+ obsnames=c("y1","y2"),
+ statenames=c("x1","x2"),
+ paramnames=c("a11","a12","a21","a22","sigma1","sigma2","nu1","nu2"),
+ params=c(a11=0.5,a12=-0.1,a21=0.2,a22=-1,nu1=0.3,nu2=0.1,sigma1=0.1,sigma2=0.3)
+) -> obj3.3 How do I deal with missing data?
Missing data are a common complication. pomp can handle NAs in the data, but if it is needed, the measurement model probability density function, dmeasure must be written so as to deal with NAs appropriately. For example, look at the following variant of the SIR model describing the influenza outbreak in a boarding school:
library(ggplot2)
-
-read.csv(text="
-B,day
-NA,0
-1,1
-6,2
-26,3
-73,4
-NA,5
-293,6
-258,7
-236,8
-191,9
-124,10
-69,11
-26,12
-11,13
-4,14
-") -> dat
-
-dat |>
- ggplot(aes(x=day, y=B)) +
- geom_line() + geom_point()library(ggplot2)
+
+read.csv(text="
+B,day
+NA,0
+1,1
+6,2
+26,3
+73,4
+NA,5
+293,6
+258,7
+236,8
+191,9
+124,10
+69,11
+26,12
+11,13
+4,14
+") -> dat
+
+dat |>
+ ggplot(aes(x=day, y=B)) +
+ geom_line() + geom_point()
Data are missing at days 0 and 5.
We create a pomp object implementing a simple SIR process model and a binomial measurement model, as in the original example. The only difference is in the dmeasure:
library(pomp)
-
-sir_step <- Csnippet("
- double dN_SI = rbinom(S,1-exp(-Beta*I/N*dt));
- double dN_IR = rbinom(I,1-exp(-gamma*dt));
- S -= dN_SI;
- I += dN_SI - dN_IR;
- R += dN_IR;
- H += dN_IR;
-")
-
-sir_init <- Csnippet("
- S = N-1;
- I = 1;
- R = 0;
- H = 0;
-")
-
-rmeas <- Csnippet("B = rbinom(H,rho);")
-
-dmeas <- Csnippet("
- if (ISNA(B)) {
- lik = (give_log) ? 0 : 1;
- } else {
- lik = dbinom(B,H,rho,give_log);
- }
-")
-
-dat |>
- pomp(time="day",t0=0,
- rprocess=euler(sir_step,delta.t=1/12),
- rinit=sir_init,
- rmeasure=rmeas,
- dmeasure=dmeas,
- accumvars="H",
- paramnames=c("Beta","gamma","N", "rho"),
- statenames=c("S","I","R","H")
- ) -> sir_nalibrary(pomp)
+
+sir_step <- Csnippet("
+ double dN_SI = rbinom(S,1-exp(-Beta*I/N*dt));
+ double dN_IR = rbinom(I,1-exp(-gamma*dt));
+ S -= dN_SI;
+ I += dN_SI - dN_IR;
+ R += dN_IR;
+ H += dN_IR;
+")
+
+sir_init <- Csnippet("
+ S = N-1;
+ I = 1;
+ R = 0;
+ H = 0;
+")
+
+rmeas <- Csnippet("B = rbinom(H,rho);")
+
+dmeas <- Csnippet("
+ if (ISNA(B)) {
+ lik = (give_log) ? 0 : 1;
+ } else {
+ lik = dbinom(B,H,rho,give_log);
+ }
+")
+
+dat |>
+ pomp(time="day",t0=0,
+ rprocess=euler(sir_step,delta.t=1/12),
+ rinit=sir_init,
+ rmeasure=rmeas,
+ dmeasure=dmeas,
+ accumvars="H",
+ paramnames=c("Beta","gamma","N", "rho"),
+ statenames=c("S","I","R","H")
+ ) -> sir_naIn the above, to check for missing data, we use the ISNA macro, which is part of R’s C API. Note that the dmeasure returns a likelihood of 1 (log likelihood 0) for the missing data. [What’s the probability of seeing something if you don’t look?] When there is an observation, it returns a binomial likelihood as usual.
Our simulations now fill in simulated values for the missing data,
-sir_na |>
- simulate(
- params=c(Beta=3,gamma=2,rho=0.9,N=2600),
- nsim=10,
- format="data.frame",
- include.data=TRUE
- ) |>
- ggplot(aes(x=day,y=B,group=.id,color=.id=="data"))+
- geom_line()+
- guides(color="none")+
- theme_bw()sir_na |>
+ simulate(
+ params=c(Beta=3,gamma=2,rho=0.9,N=2600),
+ nsim=10,
+ format="data.frame",
+ include.data=TRUE
+ ) |>
+ ggplot(aes(x=day,y=B,group=.id,color=.id=="data"))+
+ geom_line()+
+ guides(color="none")+
+ theme_bw()
and the particle filter handles the missing data correctly:
-sir_na |>
- pfilter(Np=1000,params=c(Beta=3,gamma=2,rho=0.9,N=2600)) |>
- as.data.frame() |>
- pivot_longer(-day) |>
- ggplot(aes(x=day,y=value,color=name))+
- guides(color="none")+
- geom_line()+
- facet_wrap(~name,ncol=1,scales='free_y')+
- theme_bw()sir_na |>
+ pfilter(Np=1000,params=c(Beta=3,gamma=2,rho=0.9,N=2600)) |>
+ as.data.frame() |>
+ pivot_longer(-day) |>
+ ggplot(aes(x=day,y=value,color=name))+
+ guides(color="none")+
+ geom_line()+
+ facet_wrap(~name,ncol=1,scales='free_y')+
+ theme_bw()
In the above particle filter computation, notice that the effective sample size (ESS) is full, as it should be, when the missing data contribute no information.
As described in the documentation (?pomp), pomp allows you to define “accumulator variables” that collect events occurring between observations. The example above illustrates this.
In that example, we took t0 equal to the first observation time. Sometimes we want to take t0 to be much earlier. For example, if we wish to posit that the initial state of the unobserved Markov process is distributed at t0 according to its stationary distribution, one way to achieve this is to put t0 so early that simulations will equilibrate before the first observation is made.
This is straightforward, but the presence of accumulator variables leads to a small twist. Let’s look at the boarding-school flu example to illustrate.
-library(pomp)
-
-sir_step <- Csnippet("
- double dN_SI = rbinom(S,1-exp(-Beta*I/N*dt));
- double dN_IR = rbinom(I,1-exp(-gamma*dt));
- S -= dN_SI;
- I += dN_SI - dN_IR;
- R += dN_IR;
- H += dN_IR;
- ")
-
-sir_init <- Csnippet("
- S = N-1;
- I = 1;
- R = 0;
- H = 0;
-")
-
-rmeas <- Csnippet("B = rbinom(H,rho);")
-
-dmeas <- Csnippet("if (ISNA(B)) {
- lik = (give_log) ? 0 : 1;
-} else {
- lik = dbinom(B,H,rho,give_log);
-}")
-
-read.csv(text="
-day,B
-1,3
-2,8
-3,28
-4,76
-5,222
-6,293
-7,257
-8,237
-9,192
-10,126
-11,70
-12,28
-13,12
-14,5
-") -> dat
-
-dat |>
- pomp(time="day",t0=-8,
- rprocess=euler(sir_step,delta.t=1/12),
- rinit=sir_init,
- rmeasure=rmeas,
- dmeasure=dmeas,
- accumvars="H",
- paramnames=c("Beta","gamma","N", "rho"),
- statenames=c("S","I","R","H")
- ) -> bsflulibrary(pomp)
+
+sir_step <- Csnippet("
+ double dN_SI = rbinom(S,1-exp(-Beta*I/N*dt));
+ double dN_IR = rbinom(I,1-exp(-gamma*dt));
+ S -= dN_SI;
+ I += dN_SI - dN_IR;
+ R += dN_IR;
+ H += dN_IR;
+ ")
+
+sir_init <- Csnippet("
+ S = N-1;
+ I = 1;
+ R = 0;
+ H = 0;
+")
+
+rmeas <- Csnippet("B = rbinom(H,rho);")
+
+dmeas <- Csnippet("if (ISNA(B)) {
+ lik = (give_log) ? 0 : 1;
+} else {
+ lik = dbinom(B,H,rho,give_log);
+}")
+
+read.csv(text="
+day,B
+1,3
+2,8
+3,28
+4,76
+5,222
+6,293
+7,257
+8,237
+9,192
+10,126
+11,70
+12,28
+13,12
+14,5
+") -> dat
+
+dat |>
+ pomp(time="day",t0=-8,
+ rprocess=euler(sir_step,delta.t=1/12),
+ rinit=sir_init,
+ rmeasure=rmeas,
+ dmeasure=dmeas,
+ accumvars="H",
+ paramnames=c("Beta","gamma","N", "rho"),
+ statenames=c("S","I","R","H")
+ ) -> bsfluNote that, as above, we’ve allowed for the possibility of NAs in the data.
Now let’s look at the data and some simulations from the model.
-library(ggplot2)
-
-bsflu |>
- as.data.frame() |>
- ggplot(aes(x=day,y=B))+
- geom_point()+geom_line()+
- theme_bw()library(ggplot2)
+
+bsflu |>
+ as.data.frame() |>
+ ggplot(aes(x=day,y=B))+
+ geom_point()+geom_line()+
+ theme_bw()
bsflu |>
- simulate(params=c(Beta=1.5,gamma=1,rho=0.9,N=2600),
- nsim=10,format="d",include.data=TRUE) |>
- ggplot(aes(x=day,y=B,group=.id,color=.id=="data"))+
- geom_line()+
- guides(color="none")+
- theme_bw()bsflu |>
+ simulate(params=c(Beta=1.5,gamma=1,rho=0.9,N=2600),
+ nsim=10,format="d",include.data=TRUE) |>
+ ggplot(aes(x=day,y=B,group=.id,color=.id=="data"))+
+ geom_line()+
+ guides(color="none")+
+ theme_bw()
The data plot looks fine, but what’s going on with the simulations? It’s easy to understand: we are assuming that there is one infectious host at t = t0 = -8 days. In most of the simulations, this leads to a number of new infections, which the H variable accumulates from t=-8 to t=1, the time of our first observation. But our first observation is not the number of new cases to have occurred in the last 9 days, but the number that occurred (and were reported) in the last 1 day.
We can fix this by introducing a fictitious “observation” at t=0, with missing data, i.e., by assuming we observed nothing at t=0:
library(dplyr)
-
-dat |>
- bind_rows(c(day=0,B=NA)) |>
- arrange(day) |>
- pomp(time="day",t0=-8,
- rprocess=euler(sir_step,delta.t=1/12),
- rinit=sir_init,
- rmeasure=rmeas,
- dmeasure=dmeas,
- accumvars="H",
- paramnames=c("Beta","gamma","N", "rho"),
- statenames=c("S","I","R","H")
- ) -> bsflu
-
-bsflu |>
- simulate(params=c(Beta=1.5,gamma=1,rho=0.9,N=2600),
- nsim=10,format="d",include.data=TRUE) |>
- filter(day>0) |>
- ggplot(aes(x=day,y=B,group=.id,color=.id=="data"))+
- geom_line()+
- guides(color="none")+
- theme_bw()library(dplyr)
+
+dat |>
+ bind_rows(c(day=0,B=NA)) |>
+ arrange(day) |>
+ pomp(time="day",t0=-8,
+ rprocess=euler(sir_step,delta.t=1/12),
+ rinit=sir_init,
+ rmeasure=rmeas,
+ dmeasure=dmeas,
+ accumvars="H",
+ paramnames=c("Beta","gamma","N", "rho"),
+ statenames=c("S","I","R","H")
+ ) -> bsflu
+
+bsflu |>
+ simulate(params=c(Beta=1.5,gamma=1,rho=0.9,N=2600),
+ nsim=10,format="d",include.data=TRUE) |>
+ filter(day>0) |>
+ ggplot(aes(x=day,y=B,group=.id,color=.id=="data"))+
+ geom_line()+
+ guides(color="none")+
+ theme_bw()