[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