[Pomp-commits] r776 - in www: . content vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Aug 22 10:26:17 CEST 2012
Author: kingaa
Date: 2012-08-22 10:26:17 +0200 (Wed, 22 Aug 2012)
New Revision: 776
Added:
www/content/vignettes.htm
www/vignettes/
www/vignettes/Makefile
www/vignettes/advanced_topics_in_pomp.Rnw
www/vignettes/advanced_topics_in_pomp.pdf
www/vignettes/bsmc-ricker-flat-prior.rda
www/vignettes/bsmc-ricker-normal-prior.rda
www/vignettes/complex-sir-def.rda
www/vignettes/fullnat.bst
www/vignettes/gompertz-multi-mif.rda
www/vignettes/gompertz-pfilter-guess.rda
www/vignettes/gompertz-trajmatch.rda
www/vignettes/index.html
www/vignettes/intro_to_pomp.Rnw
www/vignettes/intro_to_pomp.pdf
www/vignettes/nlf-block-boot.rda
www/vignettes/nlf-boot.rda
www/vignettes/nlf-fit-from-truth.rda
www/vignettes/nlf-fits.rda
www/vignettes/nlf-lag-tests.rda
www/vignettes/nlf-multi-short.rda
www/vignettes/plugin-C-code.rda
www/vignettes/plugin-R-code.rda
www/vignettes/pomp.bib
www/vignettes/ricker-comparison.rda
www/vignettes/ricker-first-probe.rda
www/vignettes/ricker-mif.rda
www/vignettes/ricker-probe-match.rda
www/vignettes/ricker-probe.rda
www/vignettes/sim-sim.rda
www/vignettes/sir-pomp-def.rda
www/vignettes/vectorized-C-code.rda
www/vignettes/vectorized-R-code.rda
Modified:
www/index.php
Log:
- move vignettes to R-Forge
Added: www/content/vignettes.htm
===================================================================
--- www/content/vignettes.htm (rev 0)
+++ www/content/vignettes.htm 2012-08-22 08:26:17 UTC (rev 776)
@@ -0,0 +1,8 @@
+<h3>Vignettes of package pomp</h3>
+<dl>
+<dt><a href="vignettes/intro_to_pomp.pdf">intro_to_pomp.pdf</a>:</dt>
+<dd>Introduction to <b>pomp</b></dd>
+<dt><a href=
+"vignettes/advanced_topics_in_pomp.pdf">advanced_topics_in_pomp.pdf</a>:</dt>
+<dd>Advanced topics in <b>pomp</b></dd>
+</dl>
Modified: www/index.php
===================================================================
--- www/index.php 2012-08-21 16:25:29 UTC (rev 775)
+++ www/index.php 2012-08-22 08:26:17 UTC (rev 776)
@@ -52,6 +52,7 @@
<li><a target="_blank" href="http://cran.at.r-project.org/web/packages/pomp/pomp.pdf"><i>pomp</i> manual (PDF)</a></li>
<li><a href="http://lists.r-forge.r-project.org/pipermail/pomp-announce/"><i>pomp-announce</i> mailing list archives</a></li>
<li><a href="http://cran.at.r-project.org/web/packages/pomp/NEWS">Package NEWS file</a></li>
+<li><a href="./index.php?nav=vignettes">Tutorial vignettes</a></li>
<li><a href="./index.php?nav=bib">References to the literature</a></li>
<li><a href="./index.php?nav=authors">Authors' homepages</a></li>
</ul>
@@ -68,6 +69,9 @@
<?php
$nav = $_REQUEST["nav"];
switch ($nav) {
+ case "vignettes":
+ $dfile = "content/vignettes.htm";
+ break;
case "bib":
$dfile = "content/refs.htm";
break;
Copied: www/vignettes/Makefile (from rev 774, pkg/pomp/inst/doc/Makefile)
===================================================================
--- www/vignettes/Makefile (rev 0)
+++ www/vignettes/Makefile 2012-08-22 08:26:17 UTC (rev 776)
@@ -0,0 +1,26 @@
+RSCRIPT = $(R_HOME)/bin/Rscript --vanilla
+REXE = $(R_HOME)/bin/R --vanilla
+LATEX = latex
+BIBTEX = bibtex
+PDFLATEX = pdflatex
+
+default: vignettes clean
+
+vignettes: advanced_topics_in_pomp.pdf intro_to_pomp.pdf
+
+%.tex: %.Rnw
+ $(REXE) CMD Sweave $*
+
+%.pdf: %.tex
+ $(PDFLATEX) $*
+ -$(BIBTEX) $*
+ $(PDFLATEX) $*
+ $(PDFLATEX) $*
+ $(RSCRIPT) -e "tools::compactPDF(\"$*.pdf\")";
+
+clean:
+ $(RM) *.tex *.log *.aux *.blg *.bbl *.out *.Rout *.toc *.lof *.lot
+ $(RM) *.c *.o *.so
+ $(RM) Rplots.pdf
+ $(RM) intro_to_pomp-*.pdf advanced_topics_in_pomp-*.pdf
+ $(RM) intro_to_pomp-*.png advanced_topics_in_pomp-*.png
Copied: www/vignettes/advanced_topics_in_pomp.Rnw (from rev 774, pkg/pomp/inst/doc/advanced_topics_in_pomp.Rnw)
===================================================================
--- www/vignettes/advanced_topics_in_pomp.Rnw (rev 0)
+++ www/vignettes/advanced_topics_in_pomp.Rnw 2012-08-22 08:26:17 UTC (rev 776)
@@ -0,0 +1,687 @@
+\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}
+
+\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,eval=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>>=
+binary.file <- "plugin-R-code.rda"
+if (file.exists(binary.file)) {
+ load(binary.file)
+} else {
+ tic <- Sys.time()
+simdat.Rplug <- simulate(ou2.Rplug,params=coef(ou2),nsim=1000,states=T)
+toc <- Sys.time()
+etime.Rplug <- toc-tic
+n.Rplug <- dim(simdat.Rplug)[2]
+save(etime.Rplug,n.Rplug,file=binary.file,compress='xz')
+}
+@
+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,eval=F>>=
+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,eval=F>>=
+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,eval=F>>=
+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,eval=F>>=
+tic <- Sys.time()
+simdat.Rvect <- simulate(ou2.Rvect,params=theta,states=T,nsim=1000)
+toc <- Sys.time()
+etime.Rvect <- toc-tic
+@
+<<vectorized-R-code-eval,eval=T>>=
+binary.file <- "vectorized-R-code.rda"
+if (file.exists(binary.file)) {
+ load(binary.file)
+} else {
+<<vectorized-R-code>>
+<<vectorized-R-pomp>>
+<<theta>>
+<<vectorized-R-code-sim>>
+units(etime.Rvect) <- units(etime.Rplug)
+n.Rvect <- dim(simdat.Rvect)[2]
+save(etime.Rvect,n.Rvect,file=binary.file,compress='xz')
+}
+@
+Doing \Sexpr{n.Rvect} 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{n.Rplug} simulations of \code{ou2.Rplug}, this is a \Sexpr{round(as.numeric(etime.Rplug)*n.Rvect/as.numeric(etime.Rvect)/n.Rplug)}-fold speed-up.
+
+\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,eval=F>>=
+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,eval=F>>=
+tic <- Sys.time()
+simdat.Cplug <- simulate(ou2.Cplug,params=theta,states=T,nsim=100000)
+toc <- Sys.time()
+etime.Cplug <- toc-tic
+<<plugin-C-sim-eval,eval=T,echo=F>>=
+binary.file <- "plugin-C-code.rda"
+if (file.exists(binary.file)) {
+ load(binary.file)
+} else {
+<<plugin-C-code>>
+<<plugin-C-sim>>
+ n.Cplug <- dim(simdat.Cplug)[2]
+units(etime.Cplug) <- units(etime.Rplug)
+speedup <- as.numeric(etime.Rplug/n.Rplug)/as.numeric(etime.Cplug/n.Cplug)
+save(n.Cplug,etime.Cplug,speedup,file=binary.file,compress='xz')
+}
+@
+Note that \code{ou2\_step} is written in such a way that we must specify \code{paramnames}, \code{statenames}, and \code{obsnames}.
+These \Sexpr{n.Cplug} simulations of \code{ou2.Cplug} took \Sexpr{round(etime.Cplug,2)}~\Sexpr{units(etime.Cplug)}.
+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+ interface.
+<<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,eval=F>>=
+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
+@
+<<vectorized-C-code-eval,echo=F,eval=T>>=
+binary.file <- "vectorized-C-code.rda"
+if (file.exists(binary.file)) {
+ load(binary.file)
+} else {
+<<vectorized-C-code-sim>>
+ n.Cvect <- dim(simdat.Cvect)[2]
+units(etime.Cvect) <- units(etime.Rplug)
+speedup <- as.numeric(etime.Rplug/n.Rplug)/as.numeric(etime.Cvect/n.Cvect)
+save(n.Cvect,etime.Cvect,speedup,file=binary.file,compress='xz')
+}
+@
+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.
+Doing \Sexpr{n.Cvect} 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,eval=F>>=
+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)
+<<sir-sim-eval,echo=F,eval=T>>=
+binary.file <- "sim-sim.rda"
+if (file.exists(binary.file)) {
+ load(binary.file)
+} else {
+<<sir-sim>>
+ save(sir,sims,traj,file=binary.file,compress='xz')
+}
+@
+
+\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}
Copied: www/vignettes/advanced_topics_in_pomp.pdf (from rev 774, pkg/pomp/inst/doc/advanced_topics_in_pomp.pdf)
===================================================================
(Binary files differ)
Copied: www/vignettes/bsmc-ricker-flat-prior.rda (from rev 774, pkg/pomp/inst/doc/bsmc-ricker-flat-prior.rda)
===================================================================
(Binary files differ)
Copied: www/vignettes/bsmc-ricker-normal-prior.rda (from rev 774, pkg/pomp/inst/doc/bsmc-ricker-normal-prior.rda)
===================================================================
(Binary files differ)
Copied: www/vignettes/complex-sir-def.rda (from rev 774, pkg/pomp/inst/doc/complex-sir-def.rda)
===================================================================
(Binary files differ)
Copied: www/vignettes/fullnat.bst (from rev 774, pkg/pomp/inst/doc/fullnat.bst)
===================================================================
--- www/vignettes/fullnat.bst (rev 0)
+++ www/vignettes/fullnat.bst 2012-08-22 08:26:17 UTC (rev 776)
@@ -0,0 +1,1429 @@
+%%
+%% 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:
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 776
More information about the pomp-commits
mailing list