[Pomp-commits] r142 - in pkg: data inst/doc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 2 23:50:08 CEST 2009
Author: kingaa
Date: 2009-06-02 23:50:05 +0200 (Tue, 02 Jun 2009)
New Revision: 142
Added:
pkg/inst/doc/advanced_topics_in_pomp.Rnw
pkg/inst/doc/advanced_topics_in_pomp.pdf
Removed:
pkg/inst/doc/compiled_code_in_pomp.Rnw
pkg/inst/doc/compiled_code_in_pomp.pdf
Modified:
pkg/data/euler.sir.rda
pkg/data/ou2.rda
pkg/inst/doc/Makefile
pkg/inst/doc/intro_to_pomp.Rnw
pkg/inst/doc/intro_to_pomp.pdf
Log:
more changes to documentation
reconfigure vignettes somewhat
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/Makefile
===================================================================
--- pkg/inst/doc/Makefile 2009-06-02 20:17:37 UTC (rev 141)
+++ pkg/inst/doc/Makefile 2009-06-02 21:50:05 UTC (rev 142)
@@ -6,7 +6,7 @@
default: vignettes clean
-vignettes: compiled_code_in_pomp.pdf intro_to_pomp.pdf
+vignettes: advanced_topics_in_pomp.pdf intro_to_pomp.pdf
%.pdf: %.tex
$(PDFLATEX) $*
@@ -16,4 +16,4 @@
$(EPSTOPDF) $*.eps --outfile=$*.pdf
clean:
- $(RM) *.tex *.log *.aux *.out *.Rout *.toc *.lof *.lot *-???.pdf Rplots.ps
+ $(RM) *.tex *.log *.aux *.out *.Rout *.toc *.lof *.lot *-???.pdf Rplots.ps Rplots.pdf
Copied: pkg/inst/doc/advanced_topics_in_pomp.Rnw (from rev 137, pkg/inst/doc/compiled_code_in_pomp.Rnw)
===================================================================
--- pkg/inst/doc/advanced_topics_in_pomp.Rnw (rev 0)
+++ pkg/inst/doc/advanced_topics_in_pomp.Rnw 2009-06-02 21:50:05 UTC (rev 142)
@@ -0,0 +1,276 @@
+\documentclass[10pt,reqno,final]{amsart}
+%\VignetteIndexEntry{Advanced topics in pomp}
+\usepackage{times}
+\usepackage[utf8]{inputenc}
+\usepackage{graphicx}
+\usepackage{natbib}
+
+\setlength{\textwidth}{6.25in}
+\setlength{\textheight}{8.75in}
+\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}}
+
+\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://www.umich.edu/\~{}kingaa}
+
+%% \date{\today}
+
+\newcommand\code[1]{\texttt{#1}}
+
+\begin{document}
+
+\SweaveOpts{echo=T,results=verbatim,print=F,eps=F,pdf=T}
+
+\maketitle
+
+\tableofcontents
+
+<<echo=F,results=hide>>=
+ options(keep.source=TRUE,width=60,continue=" ",prompt=" ")
+ library(pomp)
+ set.seed(5384959)
+@
+
+\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.
+
+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)
+@
+
+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=as.matrix(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+.
+
+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=as.matrix(true.p))
+dim(y)
+y[,,1:5]
+@
+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))
+@
+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.}
+
+Since many of the methods we will use require us to simulate the process and/or measurement models many times, it is a good idea to use native (compiled) codes for the computational heavy lifting.
+This can result in many-fold speedup.
+The \code{pomp} package includes some examples that use C codes.
+Here, we'll have a look at how the discrete-time 2-D Ornstein-Uhlenbeck process with normal measurement error is implemented.
+
+Recall that the unobserved Ornstein-Uhlenbeck (OU) process $X_{t}\in\mathbb{R}^2$ satisfies
+\begin{equation*}
+ X_{t} = A\,X_{t-1}+\xi_{t}.
+\end{equation*}
+The observation process is
+\begin{equation*}
+ Y_{t} = B\,X_{t}+\varepsilon_{t}.
+\end{equation*}
+In these equations, $A$ and and $B$ are 2$\times$2 constant matrices; $\xi_{t}$ and $\varepsilon_{t}$ are mutually-independent families of i.i.d. bivariate normal random variables.
+We let $\sigma\sigma^T$ be the variance-covariance matrix of $\xi_{t}$, where $\sigma$ is lower-triangular;
+likewise, we let $\tau\tau^T$ be that of $\varepsilon_{t}$.
+
+You can load a \code{pomp} object for this model with the command
+<<eval=F>>=
+data(ou2)
+@
+Here we'll examine how this object is put together.
+
+The process model simulator and density functions are as follows:
+<<>>=
+ ou2.rprocess <- function (xstart, times, params, paramnames, ...) {
+ nvar <- nrow(xstart)
+ npar <- nrow(params)
+ nrep <- ncol(xstart)
+ ntimes <- length(times)
+ ## get indices of the various parameters in the 'params' matrix
+ ## C uses zero-based indexing!
+ parindex <- match(paramnames,rownames(params))-1
+ 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)),
+ parindex = as.integer(parindex),
+ DUP = FALSE,
+ NAOK = TRUE,
+ PACKAGE = "pomp"
+ )$X,
+ dim=c(nvar,nrep,ntimes),
+ dimnames=list(rownames(xstart),NULL,NULL)
+ )
+ }
+
+ ou2.dprocess <- function (x, times, params, log, paramnames, ...) {
+ nvar <- nrow(x)
+ npar <- nrow(params)
+ nrep <- ncol(x)
+ ntimes <- length(times)
+ parindex <- match(paramnames,rownames(params))-1
+ array(
+ .C("ou2_pdf",
+ d = double(nrep*(ntimes-1)),
+ X = as.double(x),
+ par = as.double(params),
+ times = as.double(times),
+ n = as.integer(c(nvar,npar,nrep,ntimes)),
+ parindex = as.integer(parindex),
+ give_log=as.integer(log),
+ DUP = FALSE,
+ NAOK = TRUE,
+ PACKAGE = "pomp"
+ )$d,
+ dim=c(nrep,ntimes-1)
+ )
+ }
+@
+
+The call that constructs the \code{pomp} object is:
+<<>>=
+ou2 <- pomp(
+ times=seq(1,100),
+ data=rbind(
+ y1=rep(0,100),
+ y2=rep(0,100)
+ ),
+ t0=0,
+ rprocess = ou2.rprocess,
+ dprocess = ou2.dprocess,
+ dmeasure = "normal_dmeasure",
+ rmeasure = "normal_rmeasure",
+ paramnames=c(
+ "alpha.1","alpha.2","alpha.3","alpha.4",
+ "sigma.1","sigma.2","sigma.3",
+ "tau"
+ ),
+ statenames = c("x1","x2")
+ )
+@
+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.
+Read the source (file `ou2.c') to see the definitions of these functions.
+
+We'll specify some parameters:
+<<>>=
+p <- c(
+ alpha.1=0.9,alpha.2=0,alpha.3=0,alpha.4=0.99,
+ sigma.1=1,sigma.2=0,sigma.3=2,
+ tau=1,x1.0=50,x2.0=-50
+ )
+@
+
+<<>>=
+tic <- Sys.time()
+x <- simulate(ou2,params=p,nsim=500,seed=800733088)
+toc <- Sys.time()
+print(toc-tic)
+@
+
+In this example, we've written our simulators and density functions ``from scratch''.
+\code{pomp} provides ``plug-in'' facilities to make it easier to define certain kinds of models.
+These plug-ins can be used with native codes as well, as we'll see in the next example.
+
+\section{A more complex example: a seasonal SIR model}
+
+<<>>=
+euler.sir <- pomp(
+ times=seq(1/52,4,by=1/52),
+ data=rbind(measles=numeric(52*4)),
+ t0=0,
+ tcovar=seq(0,25,by=1/52),
+ covar=matrix(
+ periodic.bspline.basis(seq(0,25,by=1/52),nbasis=3,period=1,degree=3),
+ ncol=3,
+ dimnames=list(NULL,paste("seas",1:3,sep=''))
+ ),
+ delta.t=1/52/20,
+ statenames=c("S","I","R","cases","W","B","dW"),
+ paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho"),
+ covarnames=c("seas1"),
+ zeronames=c("cases"),
+ comp.names=c("S","I","R"),
+ step.fun="sir_euler_simulator",
+ rprocess=euler.simulate,
+ dens.fun="sir_euler_density",
+ dprocess=onestep.density,
+ skeleton.vectorfield="sir_ODE",
+ rmeasure="binom_rmeasure",
+ dmeasure="binom_dmeasure",
+ PACKAGE="pomp",
+ initializer=function(params, t0, comp.names, ...){
+ p <- exp(params)
+ snames <- c(
+ "S","I","R","cases","W","B",
+ "SI","SD","IR","ID","RD","dW"
+ )
+ fracs <- p[paste(comp.names,"0",sep=".")]
+ x0 <- numeric(length(snames))
+ names(x0) <- snames
+ x0[comp.names] <- round(p['pop']*fracs/sum(fracs))
+ x0
+ }
+ )
+@
+
+<<>>=
+coef(euler.sir) <- log(
+ c(
+ gamma=26,mu=0.02,iota=0.01,
+ beta1=1200,beta2=1800,beta3=600,
+ beta.sd=1e-3,
+ pop=2.1e6,
+ rho=0.6,
+ S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200
+ )
+ )
+
+<<>>=
+euler.sir <- simulate(euler.sir,nsim=1,seed=329348545L)
+@
+
+This example can be loaded via
+<<eval=F>>=
+data(euler.sir)
+@
+
+\end{document}
Copied: pkg/inst/doc/advanced_topics_in_pomp.pdf (from rev 137, pkg/inst/doc/compiled_code_in_pomp.pdf)
===================================================================
(Binary files differ)
Deleted: pkg/inst/doc/compiled_code_in_pomp.Rnw
===================================================================
--- pkg/inst/doc/compiled_code_in_pomp.Rnw 2009-06-02 20:17:37 UTC (rev 141)
+++ pkg/inst/doc/compiled_code_in_pomp.Rnw 2009-06-02 21:50:05 UTC (rev 142)
@@ -1,170 +0,0 @@
-\documentclass[10pt,reqno,final]{amsart}
-%\VignetteIndexEntry{Using compiled code in pomp}
-\usepackage{times}
-\usepackage[utf8]{inputenc}
-\usepackage{graphicx}
-\usepackage{natbib}
-
-\setlength{\textwidth}{6.25in}
-\setlength{\textheight}{8.75in}
-\setlength{\evensidemargin}{0in}
-\setlength{\oddsidemargin}{0in}
-\setlength{\topmargin}{-.35in}
-\setlength{\parskip}{.1in}
-\setlength{\parindent}{0.0in}
-
-%% $Revision: 1.11 $
-
-\title[Using compiled code]{Using compiled code in \texttt{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://www.umich.edu/\~{}kingaa}
-
-%% \date{\today}
-
-\newcommand\code[1]{\texttt{#1}}
-
-\begin{document}
-
-\SweaveOpts{echo=T,results=verbatim,print=F,eps=F,pdf=T}
-
-\maketitle
-
-\tableofcontents
-
-<<echo=F,results=hide>>=
- options(keep.source=TRUE,width=60,continue=" ",prompt=" ")
- library(pomp)
- set.seed(5384959)
-@
-
-\section{A two-dimensional Ornstein-Uhlenbeck process.}
-
-Let's look again at our example of the discrete-time 2-D Ornstein-Uhlenbeck process with normal measurement error.
-Recall that the unobserved Ornstein-Uhlenbeck (OU) process $X_{t}\in\mathbb{R}^2$ satisfies
-\begin{equation*}
- X_{t} = A\,X_{t-1}+\xi_{t}.
-\end{equation*}
-The observation process is
-\begin{equation*}
- Y_{t} = B\,X_{t}+\varepsilon_{t}.
-\end{equation*}
-In these equations, $A$ and and $B$ are 2$\times$2 constant matrices; $\xi_{t}$ and $\varepsilon_{t}$ are mutually-independent families of i.i.d. bivariate normal random variables.
-We let $\sigma\sigma^T$ be the variance-covariance matrix of $\xi_{t}$, where $\sigma$ is lower-triangular;
-likewise, we let $\tau\tau^T$ be that of $\varepsilon_{t}$.
-
-Since many of the methods we will use require us to simulate the process and/or measurement models many times, it is a good idea to use native (compiled) codes for the computational heavy lifting.
-This can result in many-fold speedup.
-The package includes some C codes that were written to implement the OU example.
-Read the source (file `ou2.c') for details.
-
-<<>>=
- ou2.rprocess <- function (xstart, times, params, paramnames, ...) {
- nvar <- nrow(xstart)
- npar <- nrow(params)
- nrep <- ncol(xstart)
- ntimes <- length(times)
- ## get indices of the various parameters in the 'params' matrix
- ## C uses zero-based indexing!
- parindex <- match(paramnames,rownames(params))-1
- 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)),
- parindex = as.integer(parindex),
- DUP = FALSE,
- NAOK = TRUE,
- PACKAGE = "pomp"
- )$X,
- dim=c(nvar,nrep,ntimes),
- dimnames=list(rownames(xstart),NULL,NULL)
- )
- }
-
- ou2.dprocess <- function (x, times, params, log, ...) {
- nvar <- nrow(x)
- npar <- nrow(params)
- nrep <- ncol(x)
- ntimes <- length(times)
- parindex <- match(paramnames,rownames(params))-1
- array(
- .C("ou2_pdf",
- d = double(nrep*(ntimes-1)),
- X = as.double(x),
- par = as.double(params),
- times = as.double(times),
- n = as.integer(c(nvar,npar,nrep,ntimes)),
- parindex = as.integer(parindex),
- give_log=as.integer(log),
- DUP = FALSE,
- NAOK = TRUE,
- PACKAGE = "pomp"
- )$d,
- dim=c(nrep,ntimes-1)
- )
- }
-
-ou2 <- pomp(
- times=seq(1,100),
- data=rbind(
- y1=rep(0,100),
- y2=rep(0,100)
- ),
- t0=0,
- rprocess = ou2.rprocess,
- dprocess = ou2.dprocess,
- dmeasure = "normal_dmeasure",
- rmeasure = "normal_rmeasure",
- obsnames=c("y1","y2"),
- paramnames=c(
- "alpha.1","alpha.2","alpha.3","alpha.4",
- "sigma.1","sigma.2","sigma.3",
- "tau"
- ),
- statenames=c("x1","x2")
- )
-@
-
-We'll specify some parameters:
-<<>>=
-p <- c(
- alpha.1=0.9,alpha.2=0,alpha.3=0,alpha.4=0.99,
- sigma.1=1,sigma.2=0,sigma.3=2,
- tau=1,x1.0=50,x2.0=-50
- )
-@
-
-<<>>=
- tic <- Sys.time()
- ou2 <- simulate(ou2,params=p,nsim=1000,seed=800733088)
- toc <- Sys.time()
- print(toc-tic)
-@
-
-Fig.~\ref{fig:ou_process} plots the data.
-
-\begin{figure}
-<<fig=T,echo=F>>=
- data(ou2)
- plot(ou2)
-@
- \caption{One realization of the two-dimensional OU process.}
- \label{fig:ou_process}
-\end{figure}
-
-The \code{pomp} object we just created is included in the package: use \code{data(ou2)} to retrieve it.
-
-\end{document}
Deleted: pkg/inst/doc/compiled_code_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.Rnw
===================================================================
--- pkg/inst/doc/intro_to_pomp.Rnw 2009-06-02 20:17:37 UTC (rev 141)
+++ pkg/inst/doc/intro_to_pomp.Rnw 2009-06-02 21:50:05 UTC (rev 142)
@@ -12,9 +12,9 @@
\setlength{\topmargin}{-.35in}
\setlength{\parskip}{.1in}
\setlength{\parindent}{0.0in}
+\setcounter{secnumdepth}{1}
+\setcounter{tocdepth}{1}
-%% $Revision: 1.12 $
-
\title[Introduction to \texttt{pomp}]{Introduction to \texttt{pomp} by example}
\author[A. A. King]{Aaron A. King}
@@ -39,14 +39,15 @@
\tableofcontents
<<echo=F,results=hide>>=
- options(keep.source=TRUE,width=60,continue=" ",prompt=" ")
+ options(keep.source=TRUE,continue=" ",prompt=" ")
library(pomp)
set.seed(5384959)
@
\section{A first example: the two-dimensional Ornstein-Uhlenbeck process.}
-To keep things simple, we will study a discrete-time process.
+To begin with, for simplicity, we will study a discrete-time process.
+Later we'll look at a continuous-time model.
The \code{pomp} package is designed with continuous-time processes in mind, but the associated computational effort is typically greater, and the additional complexities are best postponed until the structure and usage of the package is understood.
The unobserved Ornstein-Uhlenbeck (OU) process $X_{t}\in\mathbb{R}^2$ satisfies
\begin{equation*}
@@ -60,9 +61,25 @@
We let $\sigma\sigma^T$ be the variance-covariance matrix of $\xi_{t}$, where $\sigma$ is lower-triangular;
likewise, we let $\tau\tau^T$ be that of $\varepsilon_{t}$.
-In order to specify this partially-observed Markov process, we must implement both the process model (i.e., the unobserved process) and the measurement model (the observation process).
-In particular, we will need to be able to both simulate from and compute the p.d.f. of each of these models.
-In \code{pomp}, one constructs an object of class \code{pomp} to hold functions that do all of these things, together with data and other information.
+\subsection{Defining a partially observed Markov process.}
+
+In order to fully specify this partially-observed Markov process, we must implement both the process model (i.e., the unobserved process) and the measurement model (the observation process).
+That is, we would like to be able to:
+\begin{enumerate}
+\item \label{it:rproc} simulate from the process model, i.e., make a random draw from $X_{t+1}\,\vert\,X_{t}=x$ for arbitrary $x$ and $t$,
+\item \label{it:dproc} compute the probability density function (pdf) of state transitions, i.e., compute $f(X_{t+1}=x'\,\vert\,X_{t}=x)$ for arbitrary $x$, $x'$, and $t$,
+\item \label{it:rmeas} simulate from the measurement model, i.e., make a random draw from $Y_{t}\,\vert\,X_{t}=x$ for arbitrary $x$ and $t$, and
+\item \label{it:dmeas} compute the measurement model pdf, i.e., $f(Y_{t}=y\,\vert\,X_{t}=x)$ for arbitrary $x$, $y$, and $t$.
+\end{enumerate}
+For this simple model, all this is easy enough.
+In general, it will be difficult to do some of these things.
+Depending on what we wish to accomplish, however, we may not need all of these capabilities.
+For example, to simulate data, all we need is \ref{it:rproc} and \ref{it:rmeas}.
+To run a particle filter (and hence to use iterated filtering, \code{mif}), one needs \ref{it:rproc} and \ref{it:rmeas}.
+To do MCMC, one needs \ref{it:dproc} and \ref{it:dmeas}.
+Nonlinear forecasting (\code{nlf}) requires \ref{it:rproc} and \ref{it:rmeas}.
+In \code{pomp}, one constructs an object of class \code{pomp} by specifying functions to do some or all of \ref{it:rproc}--\ref{it:dmeas}, along with data and other information.
+The package provides algorithms for fitting the models to the data, for simulating the models, studying deterministic skeletons, and so on.
The documentation (\code{?pomp}) spells out the usage of the \code{pomp} constructor, including detailed specifications for all its arguments and a worked example.
\subsection{Building the \code{pomp} object}
@@ -186,6 +203,10 @@
This is available: one can optionally specify an alternative initializer.
See \code{?pomp} for details.
+<<echo=F>>=
+data(ou2)
+@
+
The pomp object \code{ou2} we just constructed has no data.
If we simulate the model, we'll obtain another \code{pomp} object just like \code{ou2}, but with the \code{data} slot filled with simulated data:
<<>>=
@@ -195,7 +216,7 @@
Here, we actually ran 1000 simulations: the default behavior of \code{simulate} is to return a list of \code{pomp} objects.
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+.
+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:
<<eval=F>>=
as(ou2,'data.frame')
@@ -205,7 +226,13 @@
<<eval=F>>=
data.array(ou2)
time(ou2)
+time(ou2,t0=TRUE)
@
+One can read and change parameters associated with the \code{pomp} object using
+<<eval=F>>=
+coef(ou2)
+coef(ou2,c("sigma.1","sigma.2")) <- c(1,0)
+@
One can also plot a \code{pomp} object (Fig.~\ref{fig:ou2}).
\begin{figure}
@@ -217,61 +244,15 @@
\caption{
One can plot a \code{pomp} object.
This shows the result of \code{plot(ou2)}.
- \label{fig:ou2}
}
+ \label{fig:ou2}
\end{figure}
-\subsection{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.
-Here, I'll just introduce each of the methods that make up this interface.
-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.
-<<>>=
-x0 <- init.state(ou2,params=true.p)
-x0
-@
+\section{Particle filtering.}
-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=as.matrix(x0),times=c(0,time(ou2)),params=as.matrix(true.p))
-dim(x)
-@
-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+.
-
-The \code{rmeasure} method gives access to the measurement model simulator:
-<<>>=
-y <- rmeasure(ou2,x=x[,,-1,drop=F],times=time(ou2),params=as.matrix(true.p))
-dim(y)
-@
-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[,,2:5,drop=F],times=time(ou2)[1:4],params=as.matrix(true.p))
-@
-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.
-
-Since the codes above that implement the model are written in \R, they're not particularly fast.
-To maximize computational efficiency, you'll frequently want to instead use compiled codes.
-For more information on this topic, see the vignette, \emph{Using compiled code in \code{pomp}}.
-For the remainder of this vignette, we'll use a version of \code{ou2} in which each of the four functions --- \code{dprocess}, \code{rprocess}, \code{dmeasure}, \code{rmeasure} --- is written in C.
-This \code{pomp} object is provided in the package; just do
-<<>>=
-data(ou2)
-@
-to retrieve it.
-The C codes that make it up are included in the source distribution of the package (look in the file `ou2.c').
-Another example, an SIR model of disease transmission, is included in the package `examples' directory.
-
-\section{Particle filter.}
-
<<echo=F>>=
set.seed(74094853)
@
@@ -332,9 +313,9 @@
tau <- diag(true.p['tau'],2,2)
fit2 <- kalman.filter(y,x0,a,b,sigma,tau)
@
-In this case, the Kalman filter gives us a log likelihood of \code{fit2\$loglik=\Sexpr{round(fit2$loglik,1)}}, while the particle filter gives us \code{fit1\$loglik=\Sexpr{round(fit1$loglik,1)}}.
+In this case, the Kalman filter gives us a log likelihood of \code{fit2\$loglik=\Sexpr{round(fit2$loglik,2)}}, while the particle filter gives us \code{fit1\$loglik=\Sexpr{round(fit1$loglik,2)}}.
-\section{The maximum likelihood by iterated filtering (MIF) algorithm}
+\section{Iterated filtering: the MIF algorithm}
The MIF algorithm works by modifying the model slightly.
It replaces the model we are interested in fitting --- which has time-invariant parameters --- with a model that is just the same except that its parameters take a random walk in time.
@@ -409,4 +390,6 @@
\label{fig:mifsim}
\end{figure}
+\section{Nonlinear forecasting}
+
\end{document}
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
More information about the pomp-commits
mailing list