[Pomp-commits] r144 - in pkg: data inst/doc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jun 3 18:24:00 CEST 2009
Author: kingaa
Date: 2009-06-03 18:24:00 +0200 (Wed, 03 Jun 2009)
New Revision: 144
Modified:
pkg/data/euler.sir.rda
pkg/inst/doc/advanced_topics_in_pomp.Rnw
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.Rnw
pkg/inst/doc/intro_to_pomp.pdf
Log:
improvements to vignettes
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/advanced_topics_in_pomp.Rnw
===================================================================
--- pkg/inst/doc/advanced_topics_in_pomp.Rnw 2009-06-02 21:50:50 UTC (rev 143)
+++ pkg/inst/doc/advanced_topics_in_pomp.Rnw 2009-06-03 16:24:00 UTC (rev 144)
@@ -6,15 +6,18 @@
\usepackage{natbib}
\setlength{\textwidth}{6.25in}
-\setlength{\textheight}{8.75in}
+\setlength{\textheight}{8.7in}
\setlength{\evensidemargin}{0in}
\setlength{\oddsidemargin}{0in}
\setlength{\topmargin}{-.35in}
\setlength{\parskip}{.1in}
\setlength{\parindent}{0.0in}
-\title[Using compiled code]{Using compiled code in \texttt{pomp}}
+\newcommand\code[1]{\texttt{#1}}
+\newcommand{\R}{\textsf{R}}
+\title[Advanced topics]{Advanced topics in \texttt{pomp}}
+
\author[A. A. King]{Aaron A. King}
\address{
@@ -30,22 +33,22 @@
%% \date{\today}
-\newcommand\code[1]{\texttt{#1}}
-
\begin{document}
-\SweaveOpts{echo=T,results=verbatim,print=F,eps=F,pdf=T}
+\SweaveOpts{echo=T,results=verbatim,print=F,eps=F,pdf=T,keep.source=T}
\maketitle
\tableofcontents
<<echo=F,results=hide>>=
- options(keep.source=TRUE,width=60,continue=" ",prompt=" ")
+ glop <- options(keep.source=TRUE,width=60,continue=" ",prompt=" ")
library(pomp)
set.seed(5384959)
@
+This document serves to introduce the low-level interface to \code{pomp} objects and to give some examples of the use of native (C or FORTRAN) codes in \code{pomp}.
+
\section{The low-level interface}
There is a low-level interface to \code{pomp} objects, primarily designed for package developers.
@@ -85,9 +88,15 @@
@
The \code{dmeasure} and \code{dprocess} methods give access to the measurement and process model densities, respectively.
<<>>=
-dprocess(ou2,x[,,36:41,drop=F],times=time(ou2)[35:40],params=as.matrix(true.p))
-dmeasure(ou2,y=y[,1,1:4],x=x[,,1:4,drop=F],times=time(ou2)[1:4],params=as.matrix(true.p))
+fp <- dprocess(ou2,x=x,times=time(ou2),params=as.matrix(true.p))
+dim(fp)
+fp[,36:40]
@
+<<>>=
+fm <- dmeasure(ou2,y=y[,1,],x=x,times=time(ou2),params=as.matrix(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 sanity checks on the inputs and outputs.
\section{Acceleration using native codes.}
@@ -184,7 +193,8 @@
"sigma.1","sigma.2","sigma.3",
"tau"
),
- statenames = c("x1","x2")
+ statenames = c("x1","x2"),
+ PACKAGE="pomp"
)
@
Notice that the process model is implemented using using \verb+.C+, while the measurement model is specified by giving the names of native C routines.
@@ -212,6 +222,20 @@
\section{A more complex example: a seasonal SIR model}
+The SIR model is a mainstay of theoretical epidemiology.
+It has the deterministic skeleton
+\begin{equation*}
+ \begin{aligned}
+ &\frac{dS}{dt}=\mu\,(N-S)+\beta(t)\,\frac{I}{N}\,S\\
+ &\frac{dI}{dt}=\beta(t)\,\frac{I}{N}\,S-\gamma\,I-\mu\,I\\
+ &\frac{dR}{dt}=\gamma\,I-\mu\,R\\
+ \end{aligned}
+\end{equation*}
+Here $N=S+I+R$ is the (constant) population size and $\beta$ is a time-dependent contact rate.
+We'll assume that the contact rate is periodic and implement it as a covariate.
+We'll implement a stochastic version of this model using an Euler-multinomial approximation to the continuous-time Markov process.
+As an additonal wrinkle, we'll assume that the rate of the infection process $\beta\,I/N$ is perturbed by white noise.
+
<<>>=
euler.sir <- pomp(
times=seq(1/52,4,by=1/52),
@@ -229,10 +253,10 @@
covarnames=c("seas1"),
zeronames=c("cases"),
comp.names=c("S","I","R"),
+ rprocess=euler.simulate,
step.fun="sir_euler_simulator",
- rprocess=euler.simulate,
+ dprocess=onestep.density,
dens.fun="sir_euler_density",
- dprocess=onestep.density,
skeleton.vectorfield="sir_ODE",
rmeasure="binom_rmeasure",
dmeasure="binom_dmeasure",
@@ -273,4 +297,8 @@
data(euler.sir)
@
+<<echo=F,results=hide>>=
+ options(glop)
+@
+
\end{document}
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.Rnw
===================================================================
--- pkg/inst/doc/intro_to_pomp.Rnw 2009-06-02 21:50:50 UTC (rev 143)
+++ pkg/inst/doc/intro_to_pomp.Rnw 2009-06-03 16:24:00 UTC (rev 144)
@@ -15,31 +15,31 @@
\setcounter{secnumdepth}{1}
\setcounter{tocdepth}{1}
+\newcommand\code[1]{\texttt{#1}}
+\newcommand{\R}{\textsf{R}}
+
\title[Introduction to \texttt{pomp}]{Introduction to \texttt{pomp} by example}
\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}
+\email{kingaa at umich.edu}
\urladdr{http://www.umich.edu/\~{}kingaa}
%% \date{\today}
-\newcommand\code[1]{\texttt{#1}}
-\newcommand{\R}{\textsf{R}}
-
\begin{document}
-\SweaveOpts{echo=T,results=verbatim,print=F,eps=F,pdf=T}
+\SweaveOpts{echo=T,results=verbatim,print=F,eps=F,pdf=T,keep.source=T}
\maketitle
\tableofcontents
<<echo=F,results=hide>>=
- options(keep.source=TRUE,continue=" ",prompt=" ")
+ glop <- options(keep.source=TRUE,width=60,continue=" ",prompt=" ")
library(pomp)
set.seed(5384959)
@
@@ -84,7 +84,10 @@
\subsection{Building the \code{pomp} object}
-We build a \code{pomp} object by specifying the four basic elements mentioned above.
+We build a \code{pomp} object by specifying some or all of the four basic elements mentioned above.
+We'll go through this in some detail here, writing each of these functions from scratch.
+Later, we'll look at some shortcuts that the package provides to streamline the process.
+
First, we write a function that implements the process model simulator.
In this function, we assume that:
\begin{enumerate}
@@ -103,8 +106,13 @@
for (k in 2:ntimes) {
for (j in 1:nreps) {
eps <- rnorm(2,mean=0,sd=1)
- x['x1',j,k] <- params['alpha.1',j]*x['x1',j,k-1]+params['alpha.3',j]*x['x2',j,k-1]+params['sigma.1',j]*eps[1]
- x['x2',j,k] <- params['alpha.2',j]*x['x1',j,k-1]+params['alpha.4',j]*x['x2',j,k-1]+params['sigma.2',j]*eps[1]+params['sigma.3',j]*eps[2]
+ x['x1',j,k] <- params['alpha.1',j]*x['x1',j,k-1]+
+ params['alpha.3',j]*x['x2',j,k-1]+
+ params['sigma.1',j]*eps[1]
+ x['x2',j,k] <- params['alpha.2',j]*x['x1',j,k-1]+
+ params['alpha.4',j]*x['x2',j,k-1]+
+ params['sigma.2',j]*eps[1]+
+ params['sigma.3',j]*eps[2]
}
}
x
@@ -115,26 +123,41 @@
When this function is called, in the course of any algorithm that uses it, some basic error checks will be performed.
Next, we write a function that computes the likelihoods of a set of process model state transitions.
-Again, pay special attention to the structure of the input arguments and return value.
+Critically, in writing this function, we are allowed to assume that the transition from \verb+x[,j,k]+ to \verb!x[,j,k+1]! is elementary, i.e., that no transition has occurred between times \verb+times[k]+ and \verb!times[k+1]!.
+If we weren't allowed to make this assumption, it would be very difficult indeed to write the transition probabilities.
+Indeed, were we able to do so, we'd have no need for \code{pomp}!
<<>>=
ou2.dprocess <- function (x, times, params, log, ...) {
## this function simulates two discrete-time OU processes
nreps <- ncol(x)
ntimes <- length(times)
eps <- numeric(2)
- d <- array(0,dim=c(nreps,ntimes-1))
+ f <- array(0,dim=c(nreps,ntimes-1))
for (k in 2:ntimes) {
for (j in 1:nreps) {
- eps[1] <- (x['x1',j,k]-params['alpha.1',j]*x['x1',j,k-1]-params['alpha.3',j]*x['x2',j,k-1])/params['sigma.1',j]
- eps[2] <- (x['x2',j,k]-params['alpha.2',j]*x['x1',j,k-1]-params['alpha.4',j]*x['x2',j,k-1]-params['sigma.2',j]*eps[1])/params['sigma.3',j]
- d[j,k-1] <- sum(dnorm(eps,mean=0,sd=1,log=TRUE),na.rm=T)-sum(log(params[c('sigma.1','sigma.3'),j]))
+ eps[1] <- x['x1',j,k]-params['alpha.1',j]*x['x1',j,k-1]-
+ params['alpha.3',j]*x['x2',j,k-1]
+ eps[2] <- x['x2',j,k]-params['alpha.2',j]*x['x1',j,k-1]-
+ params['alpha.4',j]*x['x2',j,k-1]-
+ params['sigma.2',j]/params['sigma.1',j]*eps[1]
+ f[j,k-1] <- sum(
+ dnorm(
+ x=eps,
+ mean=0,
+ sd=params[c("sigma.1","sigma.3"),j],
+ log=TRUE
+ ),
+ na.rm=T
+ )
}
}
- if (log) d else exp(d)
+ if (log) f else exp(f)
}
@
In this function, \code{times} and \code{params} are as before, and \code{x} is a rank-3 array.
-In practice, you can think of \code{x} as an array that might have been generated by a call to the \code{rprocess} function above.
+Notice that \code{dprocess} returns a rank-2 array (\code{f}) with dimensions \code{ncol(params)}$times$\code{(length(times)-1)}.
+Each row of \code{f} corresponds to a different parameter point; each column corresponds to a different transition.
+You can think of \code{x} as an array that might have been generated by a call to the \code{rprocess} function above.
Third, we write a measurement model simulator.
In this function, \code{x}, \code{t}, and \code{params} are states, time, and parameters, but they have a different form from those above.
@@ -203,8 +226,17 @@
This is available: one can optionally specify an alternative initializer.
See \code{?pomp} for details.
-<<echo=F>>=
+<<echo=F,results=hide>>=
+## this is a check to make sure the different implementations of 'ou2' are equivalent
+x <- simulate(ou2,params=true.p,nsim=2,states=T,obs=T)
+new.fp <- dprocess(ou2,x=x$states,params=cbind(true.p,true.p),times=time(ou2,t0=T),log=T)
+new.fm <- dmeasure(ou2,x=x$states,y=x$obs[,1,],params=cbind(true.p,true.p),times=time(ou2,t0=T),log=T)
+
data(ou2)
+old.fp <- dprocess(ou2,x=x$states,params=cbind(true.p,true.p),times=time(ou2,t0=T),log=T)
+old.fm <- dmeasure(ou2,x=x$states,y=x$obs[,1,],params=cbind(true.p,true.p),times=time(ou2,t0=T),log=T)
+
+if ((max(abs(old.fp-new.fp))>1e-12)||(max(abs(old.fm-new.fm))>1e-12)) stop("error in scratchy code")
@
The pomp object \code{ou2} we just constructed has no data.
@@ -215,6 +247,8 @@
@
Here, we actually ran 1000 simulations: the default behavior of \code{simulate} is to return a list of \code{pomp} objects.
+\subsection{Methods of the \code{pomp} class.}
+
There are a number of \emph{methods} that perform operations on \code{pomp} objects.
One can read the documentation on all of these by doing \verb+class?pomp+ and \verb+methods?pomp+.
For example, one can coerce a \code{pomp} object to a data frame:
@@ -248,12 +282,78 @@
\label{fig:ou2}
\end{figure}
+\subsection{Building the \code{pomp} object using plugins.}
+A fair amount of the difficulty in putting together the process model components above had to do with mundane bookkeeping:
+constructing arrays of the appropriate dimensions, filling them in the right order, making sure that the names lined up.
+The package provides some \emph{plug-in} facilities that allow the user to bypass these tedious and error-prone aspects.
+At the moment, there are plug-ins to facilitate implementation of discrete-time dynamical systems and continuous-time systems via an Euler-type algorithm.
+Let's see how we can implement the Ornstein-Uhlenbeck example using the discrete-time plugin.
+<<>>=
+ou2 <- pomp(
+ times=seq(1,100),
+ data=rbind(
+ y1=rep(0,100),
+ y2=rep(0,100)
+ ),
+ t0=0,
+ rprocess = onestep.simulate,
+ dprocess = onestep.density,
+ step.fun = function(x, t, params, delta.t, ...) {
+ eps <- rnorm(n=2)
+ with(
+ as.list(c(x,params)),
+ c(
+ x1=alpha.1*x1+alpha.3*x2+sigma.1*eps[1],
+ x2=alpha.2*x1+alpha.4*x2+sigma.2*eps[1]+sigma.3*eps[2]
+ )
+ )
+ },
+ dens.fun = function (x1, t1, x2, t2, params, ...) {
+ eps.1 <- x2['x1']-params['alpha.1']*x1['x1']-
+ params['alpha.3']*x1['x2']
+ eps.2 <- x2['x2']-params['alpha.2']*x1['x1']-
+ params['alpha.4']*x1['x2']-
+ params['sigma.2']/params['sigma.1']*eps.1
+ sum(
+ dnorm(
+ c(eps.1,eps.2),
+ mean=0,
+ sd=params[c('sigma.1','sigma.3')],
+ log=T
+ ),
+ na.rm=T
+ )
+ },
+ rmeasure = bvnorm.rmeasure,
+ dmeasure = bvnorm.dmeasure
+ )
+@
+The \code{rprocess} portion of the model is provided by the \code{onestep.simulate} plug-in.
+One specifies \verb+rprocess=onestep.simulate+ and provides an additional argument \code{step.fun} that simulates one step of the process model given one state and one set of parameters.
+The \code{dprocess} portion is specified using the \code{onestep.density} plug-in.
+In this bit, one specifies \verb+dprocess=onestep.density+ and provides an additional argument \code{dens.fun} that computes the log pdf of a transition from \verb+(x1,t1)+ to \verb+(x2,t2)+ given the parameter vector \verb+params+.
+Critically, in writing this function, one is allowed to assume that the transition from \verb+x1+ to \verb+x2+ is elementary, i.e., that no transition has occurred between times \verb+t1+ and \verb+t2+.
+See the help page (\verb+?onestep.simulate+) for complete documentation.
+<<echo=F,results=hide>>=
+## this is a check to make sure the different implementations of 'ou2' are equivalent
+x <- simulate(ou2,params=true.p,nsim=2,states=T,obs=T)
+new.fp <- dprocess(ou2,x=x$states,params=cbind(true.p,true.p),times=time(ou2,t0=T),log=T)
+new.fm <- dmeasure(ou2,x=x$states,y=x$obs[,1,],params=cbind(true.p,true.p),times=time(ou2,t0=T),log=T)
+
+data(ou2)
+old.fp <- dprocess(ou2,x=x$states,params=cbind(true.p,true.p),times=time(ou2,t0=T),log=T)
+old.fm <- dmeasure(ou2,x=x$states,y=x$obs[,1,],params=cbind(true.p,true.p),times=time(ou2,t0=T),log=T)
+
+if ((max(abs(old.fp-new.fp))>1e-12)||(max(abs(old.fm-new.fm))>1e-12)) stop("error in plugin code")
+@
+
\section{Particle filtering.}
<<echo=F>>=
+data(ou2)
set.seed(74094853)
@
@@ -350,7 +450,11 @@
)
fit <- continue(fit,Nmif=79,max.fail=100)
fitted.pars <- c("alpha.1","alpha.4","x1.0","x2.0")
-cbind(start=start.p[fitted.pars],mle=signif(coef(fit,fitted.pars),3),truth=true.p[fitted.pars])
+cbind(
+ start=start.p[fitted.pars],
+ mle=signif(coef(fit,fitted.pars),3),
+ truth=true.p[fitted.pars]
+ )
@
One can plot various diagnostics for the fitted \code{mif} object using
@@ -364,10 +468,13 @@
\begin{figure}
<<fig=T,echo=F>>=
-op <- par(mfrow=c(3,1))
-plot(conv.rec(fit,'loglik'),type='l')
-plot(conv.rec(fit,'alpha.1'),type='l')
-plot(conv.rec(fit,'alpha.4'),type='l')
+x <- conv.rec(fit,c("loglik","alpha.1","alpha.4"))
+op <- par(fig=c(0,1,0.66,0.99),mar=c(0,4,4,0))
+plot(x[,"loglik"],type='l',bty='l',xlab='',ylab=expression(log(L)),xaxt='n')
+par(fig=c(0,1,0.33,0.66),mar=c(2,4,2,0),new=T)
+plot(x[,"alpha.1"],type='l',bty='l',xlab='',ylab=expression(alpha[1]),xaxt='n')
+par(fig=c(0,1,0.0,0.33),mar=c(4,4,0,0),new=T)
+plot(x[,"alpha.4"],type='l',bty='l',xlab="MIF iteration",ylab=expression(alpha[4]))
par(op)
@
\caption{Convergence plots can be used to help diagnose convergence of the MIF algorithm.}
@@ -392,4 +499,8 @@
\section{Nonlinear forecasting}
+<<echo=F,results=hide>>=
+ options(glop)
+@
+
\end{document}
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
More information about the pomp-commits
mailing list