euler {pomp}R Documentation

Dynamical models based on stochastic Euler algorithms

Description

Facilities for simulating stochastic Markov processes based on the Euler algorithm.

Usage

euler.simulate(xstart, times, params, euler.step.fun, delta.t,
               statenames, paramnames, zeronames, tcovar, covar,
               PACKAGE)
euler.density(x, times, params, euler.dens.fun,
              statenames = character(0), paramnames = character(0),
              tcovar, covar, log = FALSE, PACKAGE)

Arguments

xstart Matrix (dimensions nvar x nsim) of states at initial time times[1].
x Matrix (dimensions nvar x nsim x ntimes) of states at times times.
times Vector of times (length ntimes) at which states are required or given.
params Matrix containing parameters of the model. The nsim columns of params correspond to those of xstart.
euler.step.fun Name of a compiled, loaded native function containing the model simulator. This function will be called to take the actual Euler steps. It must be of type "euler_step_sim" as defined in the header "euler.h", which is included with the package.
euler.dens.fun Name of a compiled, loaded native function containing the model probability density function. This function will be called to compute the probability of the actual Euler steps. It must be of type "euler_step_pdf" as defined in the header "euler.h", which is included with the package.
delta.t Time interval of Euler steps.
statenames Names of state variables, in the order they will be expected by the routine named in euler.step.fun and euler.dens.fun.
paramnames Names of parameters, in the order they will be expected by the routine named in euler.step.fun and euler.dens.fun.
zeronames Names of additional variables which will be zeroed before each time in times. These are useful, e.g., for storing accumulations of state variables.
tcovar Times at which covariates are measured.
covar Matrix of covariates. This should have dimensions length(tcovar) x ncovar, where ncovar is the number of covariates.
log logical; if TRUE, probabilities p are given as log(p).
PACKAGE an optional argument that specifies to which dynamically loaded library we restrict the search for euler.step.fun. If this is '"base"', we search in the R executable itself.

Value

euler.simulate returns a nvar x nsim x ntimes array, where nvar is the number of state variables, nsim is the number of replicate simulations (= number of columns of xstart and params), and ntimes is the length of times. If x is this array, x[,,1] will be identical to xstart; the rownames of x and xstart will also coincide.
euler.density returns a nsim x ntimes-1 array. If f is this array, f[i,j] is the likelihood of a transition from x[,i,j] to x[,i,j+1] in exactly one Euler step of duration times[j+1]-times[j].

Author(s)

Aaron A. King (kingaa at umich dot edu)

See Also

eulermultinom, pomp

Examples

## set up a lookup table for basis functions for the seasonality
tbasis <- seq(0,20,by=1/52)
basis <- periodic.bspline.basis(tbasis,nbasis=3)

## some parameters
params <- c(
            gamma=26,mu=0.2,iota=0.01,
            beta1=1200,beta2=2100,beta3=300,
            beta.sd=0.1,
            pop=2.1e5,
            rho=0.6,
            S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200
            )

## set up the pomp object
## the C codes "sir_euler_simulator" and "sir_euler_density" are included in the "examples" directory
po <- pomp(
           time=seq(1/52,20,by=1/52),
           data=rbind(measles=numeric(52*20)),
           t0=0,
           tbasis=tbasis,
           basis=basis,
           dt=1/52/20,
           statenames=c("S","I","R","cases","W","trans1"),
           paramnames=c("gamma","mu","iota","beta1","beta.sd","pop"),
           zeronames=c("cases"),
           rprocess=function(xstart,times,params,dt,tbasis,basis,
             statenames,paramnames,zeronames,...){
             euler.simulate(
                            xstart=xstart,
                            times=times,
                            params=params,
                            euler.step.fun="sir_euler_simulator",
                            delta.t=dt,
                            statenames=statenames,
                            paramnames=paramnames,
                            zeronames=zeronames,
                            tcovar=tbasis,
                            covar=basis
                            )
           },
           dprocess=function(x,times,params,log,tbasis,basis,
             statenames,paramnames,...){
             euler.density(
                           x=x,
                           times=times,
                           params=params,
                           euler.dens.fun="sir_euler_density",
                           statenames=statenames,
                           paramnames=paramnames,
                           tcovar=tbasis,
                           covar=basis,
                           log=log
                           )
           },
           rmeasure=function(x,times,params,...){
             nsims <- ncol(params)
             ntimes <- length(times)
             array(
                   data=rbinom(
                     n=nsims*ntimes,
                     size=x["cases",,],
                     prob=exp(params["rho",])
                     ),
                   dim=c(1,nsims,ntimes),
                   dimnames=list("measles",NULL,NULL)
                   )
           },
           dmeasure=function(y,x,times,params,log=FALSE,...){
             nreps <- ncol(params)
             ntimes <- length(times)
             f <- array(dim=c(nreps,ntimes))
             for (k in 1:nreps)
               f[k,] <- dbinom(
                               x=y,
                               size=x["cases",k,],
                               prob=exp(params["rho",k]),
                               log=log
                               )
             f
           },
           initializer=function(params,t0,...){
             p <- exp(params)
             pop <- p["pop"]
             fracs <- p[c("S.0","I.0","R.0")]
             x0 <- c(
                     round(pop*fracs/sum(fracs)), # make sure the three compartments sum to 'pop' initially
                     rep(0,9)   # zeros for 'cases', 'W', and the transition numbers
                     )
             names(x0) <- c(
                            "S","I","R","cases","W",
                            paste("trans",1:7,sep="")
                            )
             x0
           }
           )

## simulate from the model
tic <- Sys.time()
x <- simulate(po,params=log(params),nsim=3)
toc <- Sys.time()
print(toc-tic)
plot(x[[1]])

t <- seq(0,4/52,1/52/20)
X <- simulate(po,params=log(params),nsim=10,states=TRUE,obs=TRUE,times=t)

f <- dprocess(
              po,
              x=X$states[,,31:40],
              times=t[31:40],
              params=matrix(
                log(params),
                nrow=length(params),
                ncol=10,
                dimnames=list(names(params),NULL)
                ),
              log=TRUE
              )
apply(f,1,sum)

g <- dmeasure(
              po,
              y=t(as.matrix(X$obs[,7,])),
              x=X$states,
              times=t,
              params=matrix(
                log(params),
                nrow=length(params),
                ncol=10,
                dimnames=list(names(params),NULL)
                ),
              log=TRUE
              )
apply(g,1,sum)

[Package pomp version 0.19-1 Index]