pomp (version 0.21-3)

euler: Dynamical models based on stochastic Euler algorithms

Description

Facilities for simulating discrete-time Markov processes and continuous-time Markov processes using the Euler algorithm.

Usage

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

Arguments

xstart

Matrix (dimensions nvar x nrep) of states at initial time times[1].

x

Matrix (dimensions nvar x nrep 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 nrep columns of params correspond to those of xstart.

step.fun

This can be either an R function or a compiled, dynamically loaded native function containing the model simulator. It should be written to take a single Euler step from a single point in state space. If it is a native function, it must be of type

dens.fun

This can be either an R function or a compiled, dynamically loaded native function containing the model transition probability density function. This function will be called to compute the probability of the actual Euler steps. It must be of type

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.

covarnames

Names of columns of the matrix of covariates covar, 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).

...

if step.fn (or dens.fn) is an R function, then additional arguments will be passed to it. If step.fn (or dens.fn) is a native routine, then additional arguments are ignored.

PACKAGE

an optional argument that specifies to which dynamically loaded library we restrict the search for the native routines. If this is '"base"', we search in the R executable itself.

Value

• euler.simulate returns a nvar x nrep x ntimes array, where nvar is the number of state variables, nrep 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 nrep 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].

See Also

eulermultinom, pomp

Examples

## set up a lookup table for basis functions for the seasonality
tbasis <- seq(0,25,by=1/52)
basis <- periodic.bspline.basis(tbasis,nbasis=3)
colnames(basis) <- paste("seas",1:3,sep='')

## some parameters
params <- c(
            gamma=26,mu=0.02,iota=0.01,
            beta1=1200,beta2=1800,beta3=600,
            beta.sd=1e-3,
            pop=2.1e6,
            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", "sir_euler_density", and "sir_ODE" are included
## in the "examples" directory (file "sir.c")
po <- pomp(
           time=seq(1/52,20,by=1/52),
           data=rbind(measles=numeric(52*20)),
           t0=0,
           tcovar=tbasis,
           covar=basis,
           delta.t=1/52/20,
           statenames=c("S","I","R","cases","W","B","dW"),
           paramnames=c("gamma","mu","iota","beta1","beta.sd","pop"),
           covarnames=c("seas1"),
           zeronames=c("cases"),
           step.fun="sir_euler_simulator",
           dens.fun="sir_euler_density",
           skeleton="sir_ODE",
           rprocess=euler.simulate,
           dprocess=euler.density,
           measurement.model=measles~binom(size=cases,prob=exp(rho)),
           initializer=function(params,t0,...){
             p <- exp(params)
             with(
                  as.list(p),
                  {
                    fracs <- 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","B","SI","SD","IR","ID","RD","dW")
                    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=rbind(measles=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)