[Pomp-commits] r790 - branches/mif2/inst/doc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jan 7 17:18:09 CET 2013
Author: kingaa
Date: 2013-01-07 17:18:08 +0100 (Mon, 07 Jan 2013)
New Revision: 790
Removed:
branches/mif2/inst/doc/bsmc-ricker-flat-prior.rda
branches/mif2/inst/doc/fullnat.bst
branches/mif2/inst/doc/gompertz-multi-mif.rda
branches/mif2/inst/doc/gompertz-trajmatch.rda
branches/mif2/inst/doc/nlf-block-boot.rda
branches/mif2/inst/doc/nlf-boot.rda
branches/mif2/inst/doc/nlf-fit-from-truth.rda
branches/mif2/inst/doc/nlf-fits.rda
branches/mif2/inst/doc/nlf-lag-tests.rda
branches/mif2/inst/doc/nlf-multi-short.rda
branches/mif2/inst/doc/pomp.bib
branches/mif2/inst/doc/ricker-mif.rda
branches/mif2/inst/doc/ricker-probe-match.rda
Modified:
branches/mif2/inst/doc/advanced_topics_in_pomp.Rnw
branches/mif2/inst/doc/advanced_topics_in_pomp.pdf
branches/mif2/inst/doc/intro_to_pomp.Rnw
branches/mif2/inst/doc/intro_to_pomp.pdf
Log:
- all vignettes moved to the web
Modified: branches/mif2/inst/doc/advanced_topics_in_pomp.Rnw
===================================================================
--- branches/mif2/inst/doc/advanced_topics_in_pomp.Rnw 2013-01-07 16:07:24 UTC (rev 789)
+++ branches/mif2/inst/doc/advanced_topics_in_pomp.Rnw 2013-01-07 16:18:08 UTC (rev 790)
@@ -1,697 +1,4 @@
-\documentclass[10pt,reqno,final]{amsart}
%\VignetteIndexEntry{Advanced topics in pomp}
-\usepackage{times}
-\usepackage[utf8]{inputenc}
-\usepackage[pdftex]{graphicx}
-\usepackage[round]{natbib}
-\usepackage{paralist}
-\usepackage{float}
-
-\setlength{\textwidth}{6.25in}
-\setlength{\textheight}{8.75in}
-\setlength{\evensidemargin}{0in}
-\setlength{\oddsidemargin}{0in}
-\setlength{\topmargin}{-.35in}
-\setlength{\parskip}{.1in}
-\setlength{\parindent}{0.0in}
-\setcounter{secnumdepth}{1}
-\setcounter{tocdepth}{1}
-
-\floatstyle{ruled}
-\newfloat{textbox}{t}{tbx}
-\floatname{textbox}{Box}
-
-\newcommand{\code}[1]{\texttt{#1}}
-\newcommand{\pkg}[1]{\textsf{#1}}
-\newcommand{\R}{\textsf{R}}
-\newcommand{\expect}[1]{\mathbb{E}\left[#1\right]}
-
-\title[Advanced topics in pomp]{Advanced topics in pomp}
-
-\author[A. A. King]{Aaron A. King}
-
-\address{
- A. A. King,
- Departments of Ecology \& Evolutionary Biology and Mathematics,
- University of Michigan,
- Ann Arbor, Michigan
- 48109-1048 USA
-}
-
-\email{kingaa at umich dot edu}
-
-\urladdr{http:pomp.r-forge.r-project.org}
-
-\date{\today, \pkg{pomp}~version~\Sexpr{packageDescription("pomp",fields="Version")}}
-
-\SweaveOpts{echo=T,results=verbatim,print=F,eps=F,pdf=T,keep.source=T,resolution=100}
-
+\documentclass{article}
\begin{document}
-
-\thispagestyle{empty}
-
-\maketitle
-
-\tableofcontents
-
-<<set-opts,echo=F,results=hide>>=
- glop <- options(keep.source=TRUE,width=60,continue=" ",prompt=" ")
- library(pomp)
- pdf.options(useDingbats=FALSE)
- set.seed(5384959)
-@
-
-This document discusses accelerating \pkg{pomp} by vectorizing your codes and/or using native (C or FORTRAN) codes.
-It also introduces \pkg{pomp}'s low-level interface for code developers.
-
-\section{Accelerating your codes: vectorizing \code{rprocess} and using native codes}
-
-In the ``Introduction to pomp'' vignette, we used \emph{plug-ins} provided by the package to specify the \code{rprocess} component of partially-observed Markov process models.
-The \code{rprocess} plug-ins require you to write a simulator for a single realization of the process, for a single set of parameters, from one time to another.
-\pkg{pomp} then calls this code many times---using potentially many different parameter values, states, and times---whenever it simulates the process, computes likelihood via Monte Carlo integration, etc.
-The inference methods implemented in \pkg{pomp} are quite computationally intensive, which puts a premium on the speed of your codes.
-Sometimes, you can realize substantial speed-up of your code by vectorizing it.
-This necessitates foregoing the relative simplicity of the plug-in-based implementation and writing \code{rprocess} ``from scratch''.
-Here, we'll develop a vectorized version of \code{rprocess} in \R\ code, then we'll see what the same thing looks like coded in C.
-We'll compare these different versions in terms of their speed at simulation.
-
-We'll use a discrete-time bivariate AR(1) process with normal measurement error as our example.
-In this model, the state process $X_{t}\in\mathbb{R}^2$ satisfies
-\begin{equation}\label{eq:ou-proc}
- X_{t} = \alpha\,X_{t-1}+\sigma\,\varepsilon_{t}.
-\end{equation}
-The measurement process is
-\begin{equation}\label{eq:ou-obs}
- Y_{t} = \beta\,X_{t}+\tau\,\xi_{t}.
-\end{equation}
-In these equations, $\alpha$ and and $\beta$ are $2\times 2$ constant matrices.
-$\xi_{t}$ and $\varepsilon_{t}$ are mutually-independent families of i.i.d.\ bivariate standard normal random variables.
-$\sigma$ is a lower-triangular matrix such that $\sigma\sigma^T$ is the variance-covariance matrix of $X_{t+1}\vert X_{t}$.
-We'll assume that each component of $X$ is measured independently and with the same error, $\tau$, so that the variance-covariance matrix of $Y_{t}\vert X_{t}$ has $\tau^2$ on the diagonal and zeros elsewhere.
-
-An implementation of this model is included in the package as a \code{pomp} object;
-load it by executing \code{data(ou2)}.
-
-\subsection{An unvectorized implementation using \R\ code only}
-Before we set about vectorizing the codes, let's have a look at what a plug-in based implementation written entirely in \R\ might look like.
-<<plugin-R-code,echo=T>>=
-data(ou2)
-ou2.dat <- as.data.frame(ou2)
-
-pomp(
- data=ou2.dat[c("time","y1","y2")],
- times="time",
- t0=0,
- rprocess=discrete.time.sim(
- step.fun=function (x, t, params, ...) {
- eps <- rnorm(n=2,mean=0,sd=1) # noise terms
- xnew <- c(
- x1=params["alpha.1"]*x["x1"]+params["alpha.3"]*x["x2"]+
- params["sigma.1"]*eps[1],
- x2=params["alpha.2"]*x["x1"]+params["alpha.4"]*x["x2"]+
- params["sigma.2"]*eps[1]+params["sigma.3"]*eps[2]
- )
- names(xnew) <- c("x1","x2")
- xnew
- }
- )
- ) -> ou2.Rplug
-@
-<<plugin-R-code-sim,echo=F>>=
-tic <- Sys.time()
-simdat.Rplug <- simulate(ou2.Rplug,params=coef(ou2),nsim=1000,states=T)
-toc <- Sys.time()
-etime.Rplug <- toc-tic
-@
-Notice how we specify the process model simulator using the \code{rprocess} plug-in \code{discrete.time.sim}.
-The latter function's \code{step.fun} argument is itself a function that simulates one realization of the process for one timestep and one set of parameters.
-When we vectorize the code, we'll do many realizations at once.
-
-\subsection{Vectorizing the process simulator using \R\ code only}
-Now, to write a vectorized \code{rprocess} in \R, we must write a function that simulates \code{nrep} realizations of the unobserved process.
-Each of these realizations may start at a different point in state space and each may have a different set of parameters.
-Moreover, this function must be capable of simulating the process over an arbitrary time interval and must be capable of reporting the unobserved states at arbitrary times in that interval.
-We'll accomplish this by writing an \R\ function with arguments \code{xstart}, \code{params}, and \code{times}.
-About these inputs, we must assume:
-\begin{enumerate}
-\item \code{xstart} will be a matrix, each column of which is a vector of initial values of the state process.
- Each state variable (matrix row) will be named.
-\item \code{params} will be a matrix, the columns of which are parameter vectors.
- The parameter names will be in the matrix column-names.
-\item \code{times} will be a vector of times at which realizations of the state process are required.
- We will have \code{times[k]} $\le$ \code{times[k+1]} for all indices \code{k}, but we cannot assume that the entries of \code{times} will be unique.
-\item The initial states \code{xstart} are assumed to obtain at time \code{times[1]}.
-\end{enumerate}
-This function must return a rank-3 array, which has the realized values of the state process at the requested times.
-This array must have rownames.
-Here is one implementation of such a simulator.
-<<vectorized-R-code>>=
-ou2.Rvect.rprocess <- function (xstart, times, params, ...) {
- nrep <- ncol(xstart) # number of realizations
- ntimes <- length(times) # number of timepoints
- ## unpack the parameters (for legibility only)
- alpha.1 <- params["alpha.1",]
- alpha.2 <- params["alpha.2",]
- alpha.3 <- params["alpha.3",]
- alpha.4 <- params["alpha.4",]
- sigma.1 <- params["sigma.1",]
- sigma.2 <- params["sigma.2",]
- sigma.3 <- params["sigma.3",]
- ## x is the array of states to be returned: it must have rownames
- x <- array(0,dim=c(2,nrep,ntimes))
- rownames(x) <- rownames(xstart)
- ## xnow holds the current state values
- x[,,1] <- xnow <- xstart
- tnow <- times[1]
- for (k in seq.int(from=2,to=ntimes,by=1)) {
- tgoal <- times[k]
- while (tnow < tgoal) { # take one step at a time
- eps <- array(rnorm(n=2*nrep,mean=0,sd=1),dim=c(2,nrep))
- tmp <- alpha.1*xnow['x1',]+alpha.3*xnow['x2',]+
- sigma.1*eps[1,]
- xnow['x2',] <- alpha.2*xnow['x1',]+alpha.4*xnow['x2',]+
- sigma.2*eps[1,]+sigma.3*eps[2,]
- xnow['x1',] <- tmp
- tnow <- tnow+1
- }
- x[,,k] <- xnow
- }
- x
-}
-@
-We can put this into a pomp object that is the same as \code{ou2.Rplug} in every way except in its \code{rprocess} slot by doing
-<<vectorized-R-pomp>>=
-ou2.Rvect <- pomp(ou2.Rplug,rprocess=ou2.Rvect.rprocess)
-@
-
-Let's pick some parameters and simulate some data to see how long it takes this code to run.
-<<theta>>=
-theta <- c(
- x1.0=-3, x2.0=4,
- tau=1,
- alpha.1=0.8, alpha.2=-0.5, alpha.3=0.3, alpha.4=0.9,
- sigma.1=3, sigma.2=-0.5, sigma.3=2
- )
-@
-<<vectorized-R-code-sim>>=
-tic <- Sys.time()
-simdat.Rvect <- simulate(ou2.Rvect,params=theta,states=T,nsim=1000)
-toc <- Sys.time()
-etime.Rvect <- toc-tic
-@
-<<echo=F>>=
-units(etime.Rvect) <- units(etime.Rplug)
-@
-Doing \Sexpr{dim(simdat.Rvect)[2]} simulations of \code{ou2.Rvect} took \Sexpr{round(etime.Rvect,2)}~\Sexpr{units(etime.Rvect)}.
-Compared to the \Sexpr{round(etime.Rplug,2)}~\Sexpr{units(etime.Rplug)} it took to run \Sexpr{dim(simdat.Rplug)[2]} simulations of \code{ou2.Rplug}, this is a \Sexpr{round(as.numeric(etime.Rplug)/as.numeric(etime.Rvect))}-fold speed-up.
-
-\subsection{Using \R's byte compiler}
-
-From version {2.13}, \R\ has provided byte-compilation facilities, in the base package \code{compiler}.
-Let's see to what extent we can speed up our codes by byte-compiling the components of our \code{pomp} object.
-<<plugin-byte-code,echo=T>>=
-require(compiler)
-pomp(
- ou2.Rplug,
- rprocess=discrete.time.sim(
- step.fun=cmpfun(
- function (x, t, params, ...) {
- eps <- rnorm(n=2,mean=0,sd=1) # noise terms
- xnew <- c(
- x1=params["alpha.1"]*x["x1"]+params["alpha.3"]*x["x2"]+
- params["sigma.1"]*eps[1],
- x2=params["alpha.2"]*x["x1"]+params["alpha.4"]*x["x2"]+
- params["sigma.2"]*eps[1]+params["sigma.3"]*eps[2]
- )
- names(xnew) <- c("x1","x2")
- xnew
- },
- options=list(optimize=3)
- )
- )
- ) -> ou2.Bplug
-@
-<<plugin-byte-code-sim,echo=F>>=
-tic <- Sys.time()
-simdat.Bplug <- simulate(ou2.Bplug,params=coef(ou2),nsim=1000,states=T)
-toc <- Sys.time()
-etime.Bplug <- toc-tic
-@
-Doing these \Sexpr{dim(simdat.Bplug)[2]} simulations of \code{ou2.Bplug} took \Sexpr{round(etime.Bplug,2)}~\Sexpr{units(etime.Bplug)}.
-This is a \Sexpr{round(as.numeric(etime.Rplug)/as.numeric(etime.Bplug),1)}-fold speed-up relative to the plug-in code written in \R.
-
-We can byte-compile the vectorized \R\ code, too, and compare its performance:
-<<vectorized-byte-code>>=
-ou2.Bvect <- pomp(ou2.Rplug,rprocess=cmpfun(ou2.Rvect.rprocess,options=list(optimize=3)))
-@
-<<vectorized-byte-code-sim>>=
-tic <- Sys.time()
-simdat.Bvect <- simulate(ou2.Bvect,params=theta,states=T,nsim=1000)
-toc <- Sys.time()
-etime.Bvect <- toc-tic
-@
-<<echo=F>>=
-units(etime.Bvect) <- units(etime.Rplug)
-@
-This code shows a \Sexpr{round(as.numeric(etime.Rplug)/as.numeric(etime.Bvect))}-fold speed-up relative to the plug-in code written in \R.
-
-
-\subsection{Accelerating the code using C: a plug-in based implementation}
-
-As we've seen, we can usually acheive big accelerations using compiled native code.
-A one-step simulator written in C for use with the \code{discrete.time.sim} plug-in is included with the package and can be viewed by doing
-<<view-ou2-source,eval=F>>=
-file.show(file=system.file("examples/ou2.c",package="pomp"))
-@
-The one-step simulator is in function \code{ou2\_step}.
-Prototypes for the one-step simulator and other functions are in the \code{pomp.h} header file;
-view it by doing
-<<view-pomp.h,eval=F,results=hide>>=
-file.show(file=system.file("include/pomp.h",package="pomp"))
-@
-We can put the one-step simulator into the \code{pomp} object and simulate as before by doing
-<<plugin-C-code>>=
-ou2.Cplug <- pomp(
- ou2.Rplug,
- rprocess=discrete.time.sim("ou2_step"),
- paramnames=c(
- "alpha.1","alpha.2","alpha.3","alpha.4",
- "sigma.1","sigma.2","sigma.3",
- "tau"
- ),
- statenames=c("x1","x2"),
- obsnames=c("y1","y2")
- )
-@
-<<plugin-C-sim>>=
-tic <- Sys.time()
-simdat.Cplug <- simulate(ou2.Cplug,params=theta,states=T,nsim=100000)
-toc <- Sys.time()
-etime.Cplug <- toc-tic
-@
-Note that \code{ou2\_step} is written in such a way that we must specify \code{paramnames}, \code{statenames}, and \code{obsnames}.
-These \Sexpr{dim(simdat.Cplug)[2]} simulations of \code{ou2.Cplug} took \Sexpr{round(etime.Cplug,2)}~\Sexpr{units(etime.Cplug)}.
-<<echo=F>>=
-units(etime.Cplug) <- units(etime.Rplug)
-speedup <- as.numeric(etime.Rplug/dim(simdat.Rplug)[2])/as.numeric(etime.Cplug/dim(simdat.Cplug)[2])
-@
-This is a \Sexpr{round(speedup)}-fold speed-up relative to \code{ou2.Rplug}.
-
-
-\subsection{A vectorized C implementation}
-
-The function \code{ou2\_adv} is a fully vectorized version of the simulator written in C.
-View this code by doing
-<<eval=F,results=hide>>=
-<<view-ou2-source>>
-@
-This function is called in the following \code{rprocess} function.
-Notice that the call to \code{ou2\_adv} uses the \verb+.C+ convention.
-<<vectorized-C-code>>=
-ou2.Cvect.rprocess <- function (xstart, times, params, ...) {
- nvar <- nrow(xstart)
- npar <- nrow(params)
- nrep <- ncol(xstart)
- ntimes <- length(times)
- array(
- .C("ou2_adv",
- X = double(nvar*nrep*ntimes),
- xstart = as.double(xstart),
- par = as.double(params),
- times = as.double(times),
- n = as.integer(c(nvar,npar,nrep,ntimes))
- )$X,
- dim=c(nvar,nrep,ntimes),
- dimnames=list(rownames(xstart),NULL,NULL)
- )
-}
-@
-
-The call that constructs the \code{pomp} object is:
-<<vectorized-C-code-pomp>>=
-ou2.Cvect <- pomp(
- ou2.Rplug,
- rprocess=ou2.Cvect.rprocess
- )
-@
-<<vectorized-C-code-sim,echo=T>>=
-tic <- Sys.time()
-paramnames <- c(
- "alpha.1","alpha.2","alpha.3","alpha.4",
- "sigma.1","sigma.2","sigma.3",
- "tau",
- "x1.0","x2.0"
- )
-simdat.Cvect <- simulate(ou2.Cvect,params=theta[paramnames],nsim=100000,states=T)
-toc <- Sys.time()
-etime.Cvect <- toc-tic
-@
-Note that we've had to rearrange the order of parameters here to ensure that they arrive at the native codes in the right order.
-<<echo=F>>=
-units(etime.Cvect) <- units(etime.Rplug)
-speedup <- as.numeric(etime.Rplug/dim(simdat.Rplug)[2])/as.numeric(etime.Cvect/dim(simdat.Cvect)[2])
-@
-Doing \Sexpr{dim(simdat.Cvect)[2]} simulations of \code{ou2.Cvect} took \Sexpr{round(etime.Cvect,2)}~\Sexpr{units(etime.Cvect)}, a \Sexpr{round(speedup)}-fold speed-up relative to \code{ou2.Rplug}.
-
-
-\subsection{More on native codes and plug-ins}
-
-It's possible to use native codes for \code{dprocess} and for the measurement model portions of the \code{pomp} as well.
-In the ``Introduction to pomp'' vignette, we looked at the SIR model, which we implemented using an Euler-multinomial approximation to the continuous-time Markov process.
-Here is the same model implemented using native C codes:
-<<sir-def,keep.source=T>>=
- pomp(
- data=data.frame(
- time=seq(from=1/52,to=4,by=1/52),
- reports=NA
- ),
- times="time",
- t0=0,
- ## native routine for the process simulator:
- rprocess=euler.sim(
- step.fun="_sir_euler_simulator",
- delta.t=1/52/20,
- PACKAGE="pomp"
- ),
- ## native routine for the skeleton:
- skeleton.type="vectorfield",
- skeleton="_sir_ODE",
- ## native measurement-model routines:
- rmeasure="_sir_binom_rmeasure",
- dmeasure="_sir_binom_dmeasure",
- ## name of the shared-object library containing the
- PACKAGE="pomp",
- ## the order of the observable assumed in the native routines:
- obsnames = c("reports"),
- ## the order of the state variables assumed in the native routines:
- statenames=c("S","I","R","cases","W"),
- ## the order of the parameters assumed in the native routines:
- paramnames=c(
- "gamma","mu","iota",
- "beta1","beta.sd","pop","rho",
- "S.0","I.0","R.0"
- ),
- nbasis=3L, # three seasonal basis functions
- degree=3L, # use cubic B-splines
- period=1.0, # seasonality has period 1yr
- ## designate 'cases' as an accumulator variable
- ## i.e., set it to zero after each observation
- zeronames=c("cases"),
- ## parameter transformations in native routines:
- parameter.transform="_sir_par_trans",
- parameter.inv.transform="_sir_par_untrans",
- ## some variables to be used in the initializer
- comp.names=c("S","I","R"),
- ic.names=c("S.0","I.0","R.0"),
- ## parameterization of the initial conditions:
- initializer=function(params, t0, comp.names, ic.names, ...) {
- snames <- c("S","I","R","cases","W")
- fracs <- params[ic.names]
- x0 <- numeric(length(snames))
- names(x0) <- snames
- x0[comp.names] <- round(params['pop']*fracs/sum(fracs))
- x0
- }
- ) -> sir
-@
-The source code for the native routines \verb+_sir_euler_simulator+, \verb+_sir_ODE+, \verb+_sir_binom_rmeasure+, and \verb+_sir_binom_dmeasure+ is provided with the package (in the \code{examples} directory).
-To see the source code, do
-<<view-sir-source,eval=F>>=
-file.show(file=system.file("examples/sir.c",package="pomp"))
-@
-In the \code{demo} directory is an \R\ script that shows how to compile \verb+sir.c+ into a shared-object library and link it with \R.
-Do \code{demo(sir)} to run and view this script.
-Note that the native routines for this model are included in the package, which is why we give the \verb+PACKAGE="pomp"+ argument to \code{pomp}.
-When you write your own model using native routines, you'll compile them into a dynamically-loadable library.
-In this case, you'll want to specify the name of that library using the \code{PACKAGE} argument.
-Again, refer to the SIR example included in the \code{examples} directory to see how this is done.
-
-You can also use the R package \pkg{inline} to put C or FORTRAN codes directly into your R functions.
-
-There is an important issue that arises when using native codes.
-This has to do with the order in which parameters, states, and observables are passed to these codes.
-\pkg{pomp} relies on the names (also row-names and column-names) attributes to identify variables in vectors and arrays.
-When you write a C or FORTRAN version of \code{rprocess} or \code{dmeasure} for example, you write a routine that takes parameters, state variables, and/or observables in the form of a vector.
-However, you have no control over the order in which these are given to you.
-Without some means of knowing which element of each vector corresponds to which variable, you cannot write the codes correctly.
-This is where the \code{paramnames}, \code{statenames}, \code{covarnames}, and \code{obsnames} arguments to \code{pomp} come in: use these arguments to specify the order in which your C code expects to see the parameters, state variables, covariates, and observables (data variables).
-\code{pomp} will match these names against the corresponding names attributes of vectors.
-It will then pass to your native routines index vectors you can use to locate the correct variables.
-See the source code to see how this is done.
-
-Let's specify some parameters, simulate, and compute a deterministic trajectory:
-<<sir-sim>>=
-params <- c(
- gamma=26,mu=0.02,iota=0.01,
- beta1=400,beta2=480,beta3=320,
- beta.sd=1e-3,
- pop=2.1e6,
- rho=0.6,
- S.0=26/400,I.0=0.001,R.0=1-0.001-26/400
- )
-
-sir <- simulate(sir,params=c(params,nbasis=3,degree=3,period=1),seed=3493885L)
-sims <- simulate(sir,nsim=10,obs=T)
-traj <- trajectory(sir,hmax=1/52)
-@
-
-\section{Accumulator variables}
-
-Recall the SIR example discussed in the ``Introduction to \pkg{pomp}'' vignette.
-In this example, the data consist of reported cases, which are modeled as binomial draws from the true number of recoveries having occurred since the last observation.
-In particular, suppose the zero time for the process is $t_0$ and let $t_1, t_2, \dots, t_n$ be the times at which the data $y_1, y_2, \dots, y_n$ are recorded.
-Then the $k$-th observation $y_k = C(t_{k-1},t_{k})$ is the observed number of cases in time interval $[t_{k-1},t_{k})$.
-If $\Delta_{I{\to}R}(t_{k-1},t_{k})$ is the accumulated number of recoveries (I to R transitions) in the same interval, then the model assumes
-\begin{equation*}
- y_{k}=C(t_{k-1},t_{k})\;\sim\;\mathrm{binomial}({\Delta}_{I{\to}R}(t_{k-1},t_{k}),\rho)
-\end{equation*}
-where $\rho$ is the probability a given case is actually recorded.
-
-Now, it is easy to keep track of the cumulative number of recoveries when simulating the continuous-time SIR state process;
-one simply has to add each recovery to an accumulator variable when it occurs.
-The SIR simulator codes in the ``Introduction to \pkg{pomp}'' vignette do this, storing the cumulative number of recoveries in a state variable \code{cases}, so that at any time $t$,
-\begin{equation*}
- \mathtt{cases}(t) = \text{cumulative number of recoveries having occurred in the interval $[t_0,t)$}.
-\end{equation*}
-It follows that $\Delta_{I{\to}R}(t_{k-1},t_{k})=\mathtt{cases}(t_{k})-\mathtt{cases}(t_{k-1})$.
-Does this not violate the Markov assumption upon which all the algorithms in \pkg{pomp} are based?
-Not really.
-Straightforwardly, one could augment the state process, adding $\mathtt{cases}(t_{k-1})$ to the state vector at time $t_k$.
-The state process would then become a \emph{hybrid} process, with one component (the $S$, $I$, $R$, and $\mathtt{cases}$ variables) evolving in continuous time, while the retarded $\mathtt{cases}$ variable would update discretely.
-
-It would, of course, be relatively easy to code up the model in this way, but because the need for accumulator variables is so common, \pkg{pomp} provides an easier work-around.
-Specifically, in the \code{pomp}-object constructing call to \code{pomp}, any variables named in the \code{zeronames} argument are assumed to be accumulator variables.
-At present, however, only the \code{rprocess} plug-ins and the deterministic-skeleton trajectory codes take this into account;
-setting \code{zeronames} will have no effect on custom \code{rprocess} codes.
-
-\clearpage
-\section{The low-level interface}
-
-There is a low-level interface to \code{pomp} objects, primarily designed for package developers.
-Ordinary users should have little reason to use this interface.
-In this section, each of the methods that make up this interface will be introduced.
-
-\subsection{Getting initial states}
-
-The \code{init.state} method is called to initialize the state (unobserved) process.
-It takes a vector or matrix of parameters and returns a matrix of initial states.
-<<>>=
-data(ou2)
-true.p <- coef(ou2)
-x0 <- init.state(ou2)
-x0
-new.p <- cbind(true.p,true.p,true.p)
-new.p["x1.0",] <- 1:3
-init.state(ou2,params=new.p)
-@
-
-\subsection{Simulating the process model}
-
-The \code{rprocess} method gives access to the process model simulator.
-It takes initial conditions (which need not correspond to the zero-time \code{t0} specified when the \code{pomp} object was constructed), a set of times, and a set of parameters.
-The initial states and parameters must be matrices, and they are checked for commensurability.
-The method returns a rank-3 array containing simulated state trajectories, sampled at the times specified.
-<<>>=
-x <- rprocess(ou2,xstart=x0,times=time(ou2,t0=T),params=true.p)
-dim(x)
-x[,,1:5]
-@
-Note that the dimensions of \code{x} are \verb+nvars x nreps x ntimes+, where \code{nvars} is the number of state variables, \code{nreps} is the number of simulated trajectories (which is the number of columns in the \code{params} and \code{xstart} matrices), and \code{ntimes} is the length of the \code{times} argument.
-Note also that \verb+x[,,1]+ is identical to \verb+xstart+.
-
-\subsection{Simulating the measurement model}
-
-The \code{rmeasure} method gives access to the measurement model simulator:
-<<>>=
-x <- x[,,-1,drop=F]
-y <- rmeasure(ou2,x=x,times=time(ou2),params=true.p)
-dim(y)
-y[,,1:5]
-@
-
-\subsection{Process and measurement model densities}
-
-The \code{dmeasure} and \code{dprocess} methods give access to the measurement and process model densities, respectively.
-<<>>=
-fp <- dprocess(ou2,x=x,times=time(ou2),params=true.p)
-dim(fp)
-fp[,36:40]
-@
-<<>>=
-fm <- dmeasure(ou2,y=y[,1,],x=x,times=time(ou2),params=true.p)
-dim(fm)
-fm[,36:40]
-@
-All of these are to be preferred to direct access to the slots of the \code{pomp} object, because they do error checking on the inputs and outputs.
-
-\section{Other examples}
-
-There are a number of example \code{pomp} objects included with the package.
-These can be found by running
-<<all-data-loadable,eval=F>>=
-data(package="pomp")
-@
-The \R\ scripts that generated these are included in the \code{data-R} directory of the installed package.
-The majority of these use compiled code, which can be found in the package source.
-
-\section{Pomp Builder}
-
-<<pomp-builder,eval=F>>=
-rmeas <- "
- double size = 1.0/sigma/sigma;
- double prob = 1.0/(1.0+rho*cases/size);
- reports = rnbinom(size,prob);
-"
-
-dmeas <- "
- double size = 1.0/sigma/sigma;
- double prob = 1.0/(1.0+rho*cases/size);
- lik = dnbinom(reports,size,prob,give_log);
-"
-
-stepfn <- "
- int nrate = 6;
- int nbasis = 3;
- int degree = 3;
- double period = 1.0;
- double rate[nrate]; // transition rates
- double trans[nrate]; // transition numbers
- double dW;
- double seasonality[nbasis];
- double beta;
- int k;
-
- dW = rgammawn(beta_sd,dt); // gamma noise, mean=dt, variance=(beta_sd^2 dt)
- periodic_bspline_basis_eval(t,period,degree,nbasis,seasonality);
- beta = beta1*seasonality[0]+beta2*seasonality[1]+beta3*seasonality[2];
-
- // compute the transition rates
- rate[0] = mu*popsize; // birth into susceptible class
- rate[1] = (iota+beta*I*dW/dt)/popsize; // force of infection
- rate[2] = mu; // death from susceptible class
- rate[3] = gamma; // recovery
- rate[4] = mu; // death from infectious class
- rate[5] = mu; // death from recovered class
-
- // compute the transition numbers
- trans[0] = rpois(rate[0]*dt); // births are Poisson
- reulermultinom(2,S,&rate[1],dt,&trans[1]);
- reulermultinom(2,I,&rate[3],dt,&trans[3]);
- reulermultinom(1,R,&rate[5],dt,&trans[5]);
-
- // balance the equations
- S += trans[0]-trans[1]-trans[2];
- I += trans[1]-trans[3]-trans[4];
- R += trans[3]-trans[5];
- cases += trans[3]; // cases are cumulative recoveries
- if (beta_sd > 0.0) W += (dW-dt)/beta_sd; // mean = 0, variance = dt
-"
-
-skel <- "
- int nrate = 6;
- int nbasis = 3;
- int degree = 3; // degree of seasonal basis functions
- double period = 1.0;
- double rate[nrate]; // transition rates
- double term[nrate]; // terms in the equations
- double beta;
- double seasonality[nbasis];
-
- // compute transmission rate from seasonality
- periodic_bspline_basis_eval(t,period,degree,nbasis,seasonality);
- beta = exp(log(beta1)*seasonality[0]+log(beta2)*seasonality[1]+log(beta3)*seasonality[2]);
-
- // compute the transition rates
- rate[0] = mu*popsize; // birth into susceptible class
- rate[1] = (iota+beta*I)/popsize; // force of infection
- rate[2] = mu; // death from susceptible class
- rate[3] = gamma; // recovery
- rate[4] = mu; // death from infectious class
- rate[5] = mu; // death from recovered class
-
- // compute the several terms
- term[0] = rate[0];
- term[1] = rate[1]*S;
- term[2] = rate[2]*S;
- term[3] = rate[3]*I;
- term[4] = rate[4]*I;
- term[5] = rate[5]*R;
-
- // assemble the differential equations
- DS = term[0]-term[1]-term[2];
- DI = term[1]-term[3]-term[4];
- DR = term[3]-term[5];
- Dcases = term[3]; // accumulate the new I->R transitions
- DW = 0;
-"
-
-pompBuilder(
- data=data.frame(
- reports=NA,
- time=seq(0,10,by=1/52)
- ),
- times="time",
- t0=-1/52,
- name="SIR",
- step.fn.delta.t=1/52/20,
- paramnames=c(
- "beta1","beta2","beta3","gamma","mu",
- "beta.sd","rho","popsize","iota","sigma"
- ),
- statenames=c("S","I","R","W","cases"),
- zeronames="cases",
- rmeasure=rmeas,
- dmeasure=dmeas,
- step.fn=stepfn,
- skeleton=skel,
- skeleton.type="vectorfield"
- ) -> sir
-
-simulate(
- sir,
- params=c(gamma=26,mu=0.05,beta.sd=0.1,
- rho=0.6,sigma=0.1,popsize=1e5,iota=10,
- beta1=100,beta2=120,beta3=80,
- S.0=26000,I.0=0,R.0=74000,W.0=0,cases.0=0
- )
- ) -> sir
-@
-
-<<pomp-builder-eval,echo=F,eval=T,results=hide>>=
-if (Sys.getenv("POMP_BUILD_VIGNETTES")=="yes") {
- require(pomp)
-<<pomp-builder>>
- simulate(sir) -> x
- pfilter(sir,Np=500) -> pf
- dyn.unload("SIR.so")
- print(logLik(pf))
-}
-@
-
-
-
-<<restore-opts,echo=F,results=hide>>=
-options(glop)
-@
-
\end{document}
Modified: branches/mif2/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Deleted: branches/mif2/inst/doc/bsmc-ricker-flat-prior.rda
===================================================================
(Binary files differ)
Deleted: branches/mif2/inst/doc/fullnat.bst
===================================================================
--- branches/mif2/inst/doc/fullnat.bst 2013-01-07 16:07:24 UTC (rev 789)
+++ branches/mif2/inst/doc/fullnat.bst 2013-01-07 16:18:08 UTC (rev 790)
@@ -1,1429 +0,0 @@
-%%
-%% This is file `fullnat.bst',
-%% generated with the docstrip utility.
-%%
-%% The original source files were:
-%%
-%% merlin.mbs (with options: `ay,nat,vonx,nm-init,ed-au,keyxyr,yr-par,note-yr,vol-bf,vnum-x,num-xser,jnm-x,pub-date,pp,ed,abr,xedn,etal-it,nfss,')
-%% ----------------------------------------
-%% *** ***
-%%
-%% Copyright 1994-2004 Patrick W Daly
- % ===============================================================
- % IMPORTANT NOTICE:
- % This bibliographic style (bst) file has been generated from one or
- % more master bibliographic style (mbs) files, listed above.
- %
- % This generated file can be redistributed and/or modified under the terms
- % of the LaTeX Project Public License Distributed from CTAN
- % archives in directory macros/latex/base/lppl.txt; either
- % version 1 of the License, or any later version.
- % ===============================================================
- % Name and version information of the main mbs file:
- % \ProvidesFile{merlin.mbs}[2004/02/09 4.13 (PWD, AO, DPC)]
- % For use with BibTeX version 0.99a or later
- %-------------------------------------------------------------------
- % This bibliography style file is intended for texts in ENGLISH
- % This is an author-year citation style bibliography. As such, it is
- % non-standard LaTeX, and requires a special package file to function properly.
- % Such a package is natbib.sty by Patrick W. Daly
- % The form of the \bibitem entries is
- % \bibitem[Jones et al.(1990)]{key}...
- % \bibitem[Jones et al.(1990)Jones, Baker, and Smith]{key}...
- % The essential feature is that the label (the part in brackets) consists
- % of the author names, as they should appear in the citation, with the year
- % in parentheses following. There must be no space before the opening
- % parenthesis!
- % With natbib v5.3, a full list of authors may also follow the year.
- % In natbib.sty, it is possible to define the type of enclosures that is
- % really wanted (brackets or parentheses), but in either case, there must
- % be parentheses in the label.
- % The \cite command functions as follows:
- % \citet{key} ==>> Jones et al. (1990)
- % \citet*{key} ==>> Jones, Baker, and Smith (1990)
- % \citep{key} ==>> (Jones et al., 1990)
- % \citep*{key} ==>> (Jones, Baker, and Smith, 1990)
- % \citep[chap. 2]{key} ==>> (Jones et al., 1990, chap. 2)
- % \citep[e.g.][]{key} ==>> (e.g. Jones et al., 1990)
- % \citep[e.g.][p. 32]{key} ==>> (e.g. Jones et al., p. 32)
- % \citeauthor{key} ==>> Jones et al.
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 790
More information about the pomp-commits
mailing list