[Pomp-commits] r1120 - www/vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Mar 1 04:41:10 CET 2015
Author: kingaa
Date: 2015-03-01 04:41:09 +0100 (Sun, 01 Mar 2015)
New Revision: 1120
Added:
www/vignettes/getting_started.R
www/vignettes/getting_started.Rmd
www/vignettes/getting_started.html
Modified:
www/vignettes/Makefile
www/vignettes/pomp.bib
www/vignettes/pomp.pdf
Log:
- remove advanced topics vignette
- add getting_started vignette
- update Makefile
Modified: www/vignettes/Makefile
===================================================================
--- www/vignettes/Makefile 2015-03-01 03:40:56 UTC (rev 1119)
+++ www/vignettes/Makefile 2015-03-01 03:41:09 UTC (rev 1120)
@@ -7,7 +7,7 @@
default: vignettes clean
-vignettes: advanced_topics_in_pomp.pdf advanced_topics_in_pomp.R intro_to_pomp.pdf intro_to_pomp.R
+vignettes: intro_to_pomp.pdf intro_to_pomp.R getting_started.html getting_started.R
%.R: %.Rnw
$(RSCRIPT) -e 'library(knitr); purl("$*.Rnw")'
@@ -15,6 +15,9 @@
%.tex: %.Rnw
$(RSCRIPT) -e 'library(knitr); knit("$*.Rnw")'
+%.R: %.Rmd
+ $(RSCRIPT) -e 'library(knitr); purl("$*.Rmd")'
+
%.pdf: %.tex
$(PDFLATEX) $*
-$(BIBTEX) $*
@@ -27,4 +30,3 @@
$(RM) -r cache figure
$(RM) *.c *.o *.so *.bak
$(RM) Rplots.pdf
- $(RM) -r figure
Added: www/vignettes/getting_started.R
===================================================================
--- www/vignettes/getting_started.R (rev 0)
+++ www/vignettes/getting_started.R 2015-03-01 03:41:09 UTC (rev 1120)
@@ -0,0 +1,215 @@
+## ----opts,include=FALSE--------------------------------------------------
+library(knitr)
+prefix <- "getting-started"
+opts_chunk$set(
+ progress=TRUE,
+ prompt=FALSE,tidy=FALSE,highlight=TRUE,
+ strip.white=TRUE,
+ warning=FALSE,
+ message=FALSE,
+ error=FALSE,
+ echo=TRUE,
+ cache=TRUE,
+ results='markup',
+ fig.show='asis',
+ size='small',
+ fig.lp="fig:",
+ fig.path=paste0("figure/",prefix,"-"),
+ cache.path=paste0("cache/",prefix,"-"),
+ fig.pos="h!",
+ fig.align='center',
+ fig.height=4,fig.width=6.83,
+ out.width="\\linewidth",
+ dpi=300,
+ dev='png',
+ dev.args=list(bg='transparent')
+ )
+
+## ----prelims,echo=F,cache=F----------------------------------------------
+set.seed(594709947L)
+require(ggplot2)
+require(plyr)
+require(reshape2)
+require(magrittr)
+require(pomp)
+theme_set(theme_bw())
+stopifnot(packageVersion("pomp")>="0.62-1")
+options(
+ keep.source=TRUE,
+ stringsAsFactors=FALSE,
+ encoding="UTF-8"
+ )
+
+## ----eval=FALSE----------------------------------------------------------
+## install.packages("pomp",repos="http://R-Forge.R-Project.org")
+
+## ----eval=FALSE----------------------------------------------------------
+## cat("#include <R.h>
+## void hello (void) {
+## Rprintf(\"hello world!\\n\");
+## }",file="hello.c")
+## system("R CMD SHLIB hello.c")
+## dyn.load(paste0("hello",.Platform$dynlib.ext))
+## .C("hello",PACKAGE="hello")
+
+## ----load-parus-data-----------------------------------------------------
+parus.dat <- read.csv(text="
+ year,pop
+ 1960,148
+ 1961,258
+ 1962,185
+ 1963,170
+ 1964,267
+ 1965,239
+ 1966,196
+ 1967,132
+ 1968,167
+ 1969,186
+ 1970,128
+ 1971,227
+ 1972,174
+ 1973,177
+ 1974,137
+ 1975,172
+ 1976,119
+ 1977,226
+ 1978,166
+ 1979,161
+ 1980,199
+ 1981,306
+ 1982,206
+ 1983,350
+ 1984,214
+ 1985,175
+ 1986,211"
+ )
+
+## ----parus-plot----------------------------------------------------------
+ggplot(data=parus.dat,mapping=aes(x=year,y=pop))+
+ geom_line()+geom_point()+
+ expand_limits(y=0)+
+ theme_classic()
+
+## ----logistic-step-fun---------------------------------------------------
+step.fun <- Csnippet("
+ double dW = rnorm(0,sqrt(dt));
+ N += r*N*(1-N/K)*dt+sigma*N*dW;
+ ")
+
+## ----logistic-pomp1,cache=FALSE------------------------------------------
+parus <- pomp(data=parus.dat,time="year",t0=1959,
+ rprocess=euler.sim(step.fun=step.fun,delta.t=1/365),
+ statenames="N",paramnames=c("r","K","sigma"))
+
+## ----logistic-simul1-----------------------------------------------------
+simStates <- simulate(parus,nsim=10,params=c(r=0.2,K=200,sigma=0.5,N.0=200),states=TRUE)
+
+## ----logistic-plot1,echo=FALSE-------------------------------------------
+melt(simStates) %>%
+ dcast(rep+time~variable) %>%
+ ggplot(mapping=aes(x=time,y=N,group=rep,color=factor(rep)))+
+ geom_line()+guides(color=FALSE)+
+ theme_classic()
+
+## ----logistic-rmeasure---------------------------------------------------
+rmeas <- Csnippet("
+ pop = rpois(phi*N);
+ ")
+
+## ----logistic-pomp2,cache=FALSE------------------------------------------
+parus <- pomp(parus,rmeasure=rmeas,paramnames="phi",statenames="N")
+
+## ----logistic-simul2-----------------------------------------------------
+sim <- simulate(parus,params=c(r=0.2,K=200,phi=0.5,sigma=0.5,N.0=200),
+ nsim=10,obs=TRUE,states=TRUE)
+
+## ----logistic-plot2,echo=FALSE-------------------------------------------
+sim %>% melt() %>%
+ ggplot(mapping=aes(x=time,y=value,group=rep,color=factor(rep)))+
+ geom_line()+
+ guides(color=FALSE)+scale_y_sqrt()+
+ facet_grid(variable~.,scales="free_y")
+
+sim %>% melt() %>% dcast(rep+time~variable,value.var='value') %>%
+ ggplot(mapping=aes(x=N,y=pop,color=factor(rep)))+
+ geom_point()+scale_x_sqrt()+scale_y_sqrt()+
+ coord_equal()+
+ guides(color=FALSE)
+
+## ----logistic-dmeasure---------------------------------------------------
+dmeas <- Csnippet("
+ lik = dpois(pop,phi*N,give_log);
+ ")
+
+## ----logistic-pomp3,cache=FALSE------------------------------------------
+parus <- pomp(parus,dmeasure=dmeas,paramnames="phi",statenames="N")
+
+## ----logistic-pfilter----------------------------------------------------
+pf <- pfilter(parus,Np=1000,params=c(r=0.2,K=200,phi=0.5,sigma=0.5,N.0=200))
+logLik(pf)
+
+## ----logistic-skeleton,cache=FALSE---------------------------------------
+skel <- Csnippet("
+ DN = r*N*(1-N/K);
+ ")
+
+parus <- pomp(parus,skeleton=skel,skeleton.type="vectorfield",statenames="N",paramnames=c("r","K"))
+
+## ----logistic-traj1------------------------------------------------------
+pars <- parmat(c(r=1,K=200,phi=0.5,sigma=0.5,N.0=20),5)
+pars["N.0",] <- seq(20,300,length=5)
+traj <- trajectory(parus,params=pars,times=seq(1959,1970,by=0.01))
+
+## ----logistic-plot3,echo=FALSE-------------------------------------------
+trajectory(parus,params=pars,times=seq(1959,1970,by=0.01),as.data.frame=TRUE) %>%
+ ggplot(mapping=aes(x=time,y=N,group=traj,color=traj))+
+ guides(color=FALSE)+
+ geom_line()
+
+## ----bh-stepfun----------------------------------------------------------
+bh.step <- Csnippet("
+ double eps = rlnorm(-sigma*sigma/2,sigma);
+ N = a*N/(1+b*N)*eps;
+ ")
+
+## ----bh-skeleton---------------------------------------------------------
+bh.skel <- Csnippet("
+ DN = a*N/(1+b*N);
+ ")
+
+## ----bh-pomp1,cache=FALSE------------------------------------------------
+parus.bh <- pomp(parus,rprocess=discrete.time.sim(bh.step,delta.t=1),skeleton=bh.skel,skeleton.type="map",skelmap.delta.t=1,statenames="N",paramnames=c("a","b","sigma"))
+
+## ----bh-test-------------------------------------------------------------
+coef(parus.bh) <- c(a=1.1,b=5e-4,sigma=0.5,N.0=30,phi=1)
+sim <- simulate(parus.bh)
+traj <- trajectory(parus.bh)
+pf <- pfilter(parus.bh,Np=1000)
+
+## ----logistic-partrans,cache=FALSE---------------------------------------
+partrans <- Csnippet("
+ Tr = exp(r);
+ TK = exp(K);
+ ")
+
+parinvtrans <- Csnippet("
+ Tr = log(r);
+ TK = log(K);
+ ")
+
+parus <- pomp(parus,parameter.transform=partrans,parameter.inv.transform=parinvtrans,
+ paramnames=c("r","K"))
+
+## ----logistic-partrans-test,include=FALSE,cache=FALSE--------------------
+p <- c(r=1,K=200,phi=1,N.0=200,sigma=0.5)
+coef(parus,transform=TRUE) <- partrans(parus,p,dir="inv")
+stopifnot(all.equal(p,coef(parus)))
+
+## ----tempfile-view,eval=FALSE--------------------------------------------
+## op <- options(verbose=TRUE)
+##
+## parus <- pomp(data=parus.dat,times="year",t0=1959,rprocess=euler.sim(step.fun,delta.t=1/365),
+## dmeasure=dmeas,rmeasure=rmeas,statenames="N",paramnames=c("r","phi","K","N.0","sigma"))
+##
+## options(op)
+
Added: www/vignettes/getting_started.Rmd
===================================================================
--- www/vignettes/getting_started.Rmd (rev 0)
+++ www/vignettes/getting_started.Rmd 2015-03-01 03:41:09 UTC (rev 1120)
@@ -0,0 +1,461 @@
+---
+title: "Getting started with pomp"
+author: "Aaron A. King"
+date: "02/28/2015"
+output:
+ html_document:
+ theme: flatly
+ toc: yes
+bibliography: pomp.bib
+csl: ecology.csl
+---
+
+Licensed under the [Creative Commons attribution-noncommercial license](http://creativecommons.org/licenses/by-nc/3.0).
+Please share and remix noncommercially, mentioning its origin.
+
+
+```{r opts,include=FALSE}
+library(knitr)
+prefix <- "getting-started"
+opts_chunk$set(
+ progress=TRUE,
+ prompt=FALSE,tidy=FALSE,highlight=TRUE,
+ strip.white=TRUE,
+ warning=FALSE,
+ message=FALSE,
+ error=FALSE,
+ echo=TRUE,
+ cache=TRUE,
+ results='markup',
+ fig.show='asis',
+ size='small',
+ fig.lp="fig:",
+ fig.path=paste0("figure/",prefix,"-"),
+ cache.path=paste0("cache/",prefix,"-"),
+ fig.pos="h!",
+ fig.align='center',
+ fig.height=4,fig.width=6.83,
+ out.width="\\linewidth",
+ dpi=300,
+ dev='png',
+ dev.args=list(bg='transparent')
+ )
+```
+
+
+```{r prelims,echo=F,cache=F}
+set.seed(594709947L)
+require(ggplot2)
+require(plyr)
+require(reshape2)
+require(magrittr)
+require(pomp)
+theme_set(theme_bw())
+stopifnot(packageVersion("pomp")>="0.62-1")
+options(
+ keep.source=TRUE,
+ stringsAsFactors=FALSE,
+ encoding="UTF-8"
+ )
+```
+
+## Partially observed Markov process (POMP) models
+
+As its name implies `pomp` is focused on partially observed Markov process models.
+These are also known as state-space model, hidden Markov models, and stochastic dynamical systems models.
+Such models consist of an unobserved stochastic *state process*, connected to the data via an explicit model of the observation process.
+We refer to the former as the *process model* and the latter as the *measurement model*.
+Each model is, in fact, a probability distribution.
+If we let $Y_t$ denote the measurement process at time $t$ and $X_t$ the state process, then by definition, the process model is determined by the density $f_{X_t|X_{t-1}}$ and the measurement process by the density $f_{Y_t|X_t}$.
+We think of the data, $y^*_t$, $t=1,\dots,T$ as being a single realization of the $Y$ process.
+To implement a POMP model in `pomp`, we have to specify the measurement and process distributions.
+
+Note however, that, for each of the process and the measurement model there are two distinct operations that we might wish to perform:
+
+1. we might wish to *simulate*, i.e., to draw a (pseudo)random sample from the distribution, or
+2. we might wish to *evaluate the density* itself at given values of $X_t$ and/or $Y_t$.
+
+We refer to the simulation of $f_{X_t|X_{t-1}}$ as the *rprocess* component of the POMP model, the evaluation of $f_{X_t|X_{t-1}}$ as the *dprocess* component, the simulation of $f_{Y_t|X_t}$, as the *rmeasure* component, and the evaluation of $f_{Y_t|X_t}$ as the *dmeasure* component.
+Methods that make no use of the *dprocess* component are called "plug-and-play" methods.
+At present, `pomp` is focused on such methods, so there is no reason to focus on the dprocess component any further.
+In the following, we will illustrate and explain how one specifies the rprocess, rmeasure, and dmeasure components of a model in `pomp`.
+We will illustrate this using some simple examples.
+
+-----------------------------------
+
+
+Schematic of the structure of a POMP showing causal relations and direction of information flow.
+**The key perspective to keep in mind is that the model is to be viewed as the process that generated the data.**
+
+-----------------------------------
+
+
+Another POMP model schematic, showing dependence among model variables.
+State variables, $x$, at time $t$ depend only on state variables at the previous timestep.
+Measurements, $y$, at time $t$ depend only on the state at that time.
+Formally, $\mathrm{Prob}[X_t|X_1,\dots,X_{t-1},Y_1,\dots,Y_{t-1}]=\mathrm{Prob}[X_t|X_{t-1}]$ and $\mathrm{Prob}[Y_t|X_1,\dots,X_{t},Y_1,\dots,Y_{t-1}]=\mathrm{Prob}[Y_t|X_{t}]$ for all $t=1,\dots,T$.
+
+-----------------------------------
+
+## Preliminaries
+
+### Installing the package
+
+To get started, we must install `pomp`, if it is not already installed.
+The package can be downloaded from CRAN, but the latest version is always available at the package homepage on R-Forge.
+Run the following to install this to download and install the latest version.
+```{r eval=FALSE}
+install.packages("pomp",repos="http://R-Forge.R-Project.org")
+```
+
+### Important information for Windows and Mac users.
+
+In this document, we will implement `pomp` models using the package's `Csnippet` facility.
+This allows the user to write model components using snippets of C code, which is then compiled and linked into a running `R` session.
+This typically affords a manyfold increase in computation time.
+It is possible to avoid `Csnippet`s entirely by writing model components as `R` functions, but the resulting implementations are typically too slow to be of practical use.
+
+**Note to Windows users.**
+To use `Csnippet`s, you must be able to compile C codes.
+Compilers are not typically installed on Windows systems.
+To achieve full functionality of `pomp` on a Windows computer, therefore, you must install the `Rtools` suite.
+See http://cran.r-project.org/bin/windows/Rtools for downloads and instructions.
+
+**Note to Mac users.**
+Although a variety of compilers are freely available for the Mac platform, they are not typically installed by default.
+
+The `Csnippet` functionality in essence requires that you be able to use `R CMD SHLIB`.
+To see if you can do this, execute the following in an `R` session.
+```{r eval=FALSE}
+cat("#include <R.h>
+ void hello (void) {
+ Rprintf(\"hello world!\\n\");
+ }",file="hello.c")
+system("R CMD SHLIB hello.c")
+dyn.load(paste0("hello",.Platform$dynlib.ext))
+.C("hello",PACKAGE="hello")
+```
+If you see the "hello world" message printed, your system should support `Csnippet`.
+
+
+### Loading data
+
+We will illustrate `pomp` by performing a limited data analysis on a set of bird abundance data.
+The data are from annual censuses of a population of *Parus major* (Great Tit) in Wytham Wood, Oxfordshire.
+They were retrieved as dataset #10163 from the Global Population Dynamics Database version 2 (NERC Centre for Population Biology, Imperial College, 2010).
+The original source is @McCleery1991.
+
+We load the data by doing
+```{r load-parus-data}
+parus.dat <- read.csv(text="
+ year,pop
+ 1960,148
+ 1961,258
+ 1962,185
+ 1963,170
+ 1964,267
+ 1965,239
+ 1966,196
+ 1967,132
+ 1968,167
+ 1969,186
+ 1970,128
+ 1971,227
+ 1972,174
+ 1973,177
+ 1974,137
+ 1975,172
+ 1976,119
+ 1977,226
+ 1978,166
+ 1979,161
+ 1980,199
+ 1981,306
+ 1982,206
+ 1983,350
+ 1984,214
+ 1985,175
+ 1986,211"
+ )
+```
+
+Here is a plot of these data:
+
+```{r parus-plot}
+ggplot(data=parus.dat,mapping=aes(x=year,y=pop))+
+ geom_line()+geom_point()+
+ expand_limits(y=0)+
+ theme_classic()
+```
+
+## Specifying a continuous-time process model
+
+To keep things simple, let's attempt to explain these data using a very simple model of stable but stochastic population dynamics, the logistic, or Verhulst-Pearl, equation with environmental stochasticity.
+We'll write this model as a stochastic differential equation (SDE), specifically an Ito diffusion:
+$$dN = r\,N\,\left(1-\frac{N}{K}\right)\,dt+\sigma\,N\,dW(t),$$
+where $N$ is the population size, $r$ is a fixed rate, the so-called "Malthusian parameter", $K$ is the population's "carrying capacity", $\sigma$ describes the intensity of extrinsic stochastic forcing on the system, and $dW$ is an increment of a standard Wiener process.
+[Those unfamiliar with Wiener processes and Ito diffusions will find it useful to visualize $dW(t)$, for each time $t$, as a normal random variable with mean zero and standard deviation $\sqrt{dt}$.]
+To implement this model in `pomp`, we must tell the package how to simulate this model.
+The easiest way to simulate such an SDE is via the *Euler-Maruyama* method.
+In this approximation, we take a large number of very small steps, each of duration $\Delta t$.
+At each step, we hold the right-hand side of the above equation constant, compute $\Delta N$ using that equation, and increment $N$ accordingly.
+`pomp` gives us the `euler.sim` function to assist us in implementing the Euler-Maruyama method.
+To use it, we must encode the computations that take a single step.
+The best way to do this is to write a snippet of code in the C language.
+The following is such a snippet for the stochastic logistic equation above.
+
+```{r logistic-step-fun}
+step.fun <- Csnippet("
+ double dW = rnorm(0,sqrt(dt));
+ N += r*N*(1-N/K)*dt+sigma*N*dW;
+ ")
+```
+
+This is just a snippet of C code: not all the variables are declared and the context of the snippet is not specified.
+When given this snippet, `pomp` will provide the necessary declarations and context, compile the resulting C code, dynamically link to it, and use it to simulate realizations of the process model.
+We cause all this to happen when we construct an object of class `pomp` via a call to the constructor function, which is also named `pomp`:
+
+```{r logistic-pomp1,cache=FALSE}
+parus <- pomp(data=parus.dat,time="year",t0=1959,
+ rprocess=euler.sim(step.fun=step.fun,delta.t=1/365),
+ statenames="N",paramnames=c("r","K","sigma"))
+```
+
+In the above, we've specified an Euler-Maruyama time-step of about one day.
+The `t0` argument specifies that we are going to treat the stochastic state process as being initialized in year 1959.
+We use the `statenames` and `paramnames` arguments to indicate which of the undeclared variables in the `Csnippet` `step.fun` are state variables and which are fixed parameters.
+Since `dW` is a local variable, we must provide an ordinary C declaration for it.
+Note that `dt` is a variable that is defined in the context of this snippet; it is actually provided by `pomp` and will contain the size of the Euler step.
+The `rnorm` function is part of the [`R` API](http://cran.r-project.org/doc/manuals/r-release/R-exts.html#The-R-API): see the [manual on "Writing R Extensions"](http://cran.r-project.org/doc/manuals/r-release/R-exts.html) for a description of this and the other [distribution functions provided as part of the `R` API](http://cran.r-project.org/doc/manuals/r-release/R-exts.html#Distribution-functions).
+Finally, note that the state variable `N` is over-written by this snippet: it's value when the first line is executed is overwritten by the second line.
+
+[A full set of rules for writing `pomp` C snippets is given below][rules].
+
+With the process model in place, we can simulate the process model using the `simulate` function and plot several realizations for a given set of parameters.
+
+```{r logistic-simul1}
+simStates <- simulate(parus,nsim=10,params=c(r=0.2,K=200,sigma=0.5,N.0=200),states=TRUE)
+```
+```{r logistic-plot1,echo=FALSE}
+melt(simStates) %>%
+ dcast(rep+time~variable) %>%
+ ggplot(mapping=aes(x=time,y=N,group=rep,color=factor(rep)))+
+ geom_line()+guides(color=FALSE)+
+ theme_classic()
+```
+
+
+## Specifying the measurement model
+
+Although it is the process model that is usually the focus of an investigator's interest, the measurement model is equally important because it connects the process model to the data.
+For this example, let's model the observations using a Poisson distribution:
+$$\mathrm{pop}_t \sim \mathrm{Poisson}(\phi\,N),$$
+where $\mathrm{pop}_t$ is our census population size at time $t$, $N$ is the true population size, and $\phi$ is a scaling parameter.
+
+### The rmeasure component
+
+The following snippet encodes the rmeasure component for this model.
+```{r logistic-rmeasure}
+rmeas <- Csnippet("
+ pop = rpois(phi*N);
+ ")
+```
+In this snippet of code, `N` is the state variable, `phi` is another parameter, and `pop` is the name of our observable, as defined in the dataset `parus.dat`.
+The `rpois` function is part of the [`R` API](http://cran.r-project.org/doc/manuals/r-release/R-exts.html#The-R-API): it takes a single argument---the Poisson distribution's parameter---and returns a pseudo-random draw from the Poisson distribution with that parameter.
+As with the `step.fun` snippet above, execution of this snippet causes `pop` to be overwritten.
+
+To fold this into our `pomp` object, we do
+```{r logistic-pomp2,cache=FALSE}
+parus <- pomp(parus,rmeasure=rmeas,paramnames="phi",statenames="N")
+```
+Note that, again, we must tell `pomp` which of the variables is a parameter and which is a state variable.
+We can now simulate from the full POMP model and plot the results.
+```{r logistic-simul2}
+sim <- simulate(parus,params=c(r=0.2,K=200,phi=0.5,sigma=0.5,N.0=200),
+ nsim=10,obs=TRUE,states=TRUE)
+```
+```{r logistic-plot2,echo=FALSE}
+sim %>% melt() %>%
+ ggplot(mapping=aes(x=time,y=value,group=rep,color=factor(rep)))+
+ geom_line()+
+ guides(color=FALSE)+scale_y_sqrt()+
+ facet_grid(variable~.,scales="free_y")
+
+sim %>% melt() %>% dcast(rep+time~variable,value.var='value') %>%
+ ggplot(mapping=aes(x=N,y=pop,color=factor(rep)))+
+ geom_point()+scale_x_sqrt()+scale_y_sqrt()+
+ coord_equal()+
+ guides(color=FALSE)
+```
+
+### The dmeasure component
+
+The following snippet encodes the dmeasure component.
+```{r logistic-dmeasure}
+dmeas <- Csnippet("
+ lik = dpois(pop,phi*N,give_log);
+ ")
+```
+and we fold it into the `pomp` object via
+```{r logistic-pomp3,cache=FALSE}
+parus <- pomp(parus,dmeasure=dmeas,paramnames="phi",statenames="N")
+```
+
+In the `dmeas` snippet, `dpois` again comes from the [`R` API](http://cran.r-project.org/doc/manuals/r-release/R-exts.html#The-R-API).
+It takes three arguments, the datum, the Poisson parameter, and `give_log`.
+When `give_log=0`, `dpois` returns the Poisson likelihood;
+when `give_log=1`, `dpois` returns the log of this likelihood.
+When this snippet is executed, `pomp` will provide the value of `give_log` according to its needs.
+It is the user's responsibility to make sure that the correct value is returned.
+
+Since we now have both the rprocess and the dmeasure components in place, we can perform a particle filtering operation to estimate the log likelihood at any given point in parameter space.
+For example, we can do
+```{r logistic-pfilter}
+pf <- pfilter(parus,Np=1000,params=c(r=0.2,K=200,phi=0.5,sigma=0.5,N.0=200))
+logLik(pf)
+```
+
+## A continuous-time deterministic skeleton
+
+A deterministic skeleton of the logistic model is obtained by taking $\sigma\to 0$, which leads to the deterministic differential equation
+$$\frac{dN}{dt} = r\,N\,\left(1-\frac{N}{K}\right).$$
+The following snippet encodes this deterministic skeleton and folds it into the `pomp` object.
+```{r logistic-skeleton,cache=FALSE}
+skel <- Csnippet("
+ DN = r*N*(1-N/K);
+ ")
+
+parus <- pomp(parus,skeleton=skel,skeleton.type="vectorfield",statenames="N",paramnames=c("r","K"))
+```
+Note that in this snippet, `DN` is filled with the value of the time-derivative of `N`.
+[See below for a complete set of rules for writing C snippets for `pomp`][rules].
+
+With the deterministic skeleton in place we can generate several trajectories of the skeleton using `trajectory`, as follows, and plot the result.
+```{r logistic-traj1}
+pars <- parmat(c(r=1,K=200,phi=0.5,sigma=0.5,N.0=20),5)
+pars["N.0",] <- seq(20,300,length=5)
+traj <- trajectory(parus,params=pars,times=seq(1959,1970,by=0.01))
+```
+```{r logistic-plot3,echo=FALSE}
+trajectory(parus,params=pars,times=seq(1959,1970,by=0.01),as.data.frame=TRUE) %>%
+ ggplot(mapping=aes(x=time,y=N,group=traj,color=traj))+
+ guides(color=FALSE)+
+ geom_line()
+```
+
+## Specifying a discrete-time process model and skeleton
+
+Let us demonstrate the implementation of discrete-time process models by replacing the continuous-time logistic equation used above with the stochastic Beverton-Holt model
+$$N_{t+1}=\frac{a\,N_t}{1+b\,N_t}\,\varepsilon_t,$$
+where $a$ and $b$ are parameters and $\varepsilon_t\sim\mathrm{Lognormal}(-\tfrac{1}{2}\sigma^2,\sigma)$.
+
+Since we are now dealing with a discrete-time model, `euler.sim` will no longer construct an appropriate rprocess component.
+Instead, we can use `discrete.time.sim`, which also takes a C snippet encoding the one-step simulation.
+The following snippet simulates a single iteration of the stochastic Beverton-Holt map
+```{r bh-stepfun}
+bh.step <- Csnippet("
+ double eps = rlnorm(-sigma*sigma/2,sigma);
+ N = a*N/(1+b*N)*eps;
+ ")
+```
+A corresponding skeleton is the deterministic Beverton-Holt map obtained by setting $e_t=1$ in the above equation.
+A snippet that implements this map is
+```{r bh-skeleton}
+bh.skel <- Csnippet("
+ DN = a*N/(1+b*N);
+ ")
+```
+Note that, as in the continuous case, we indicate the new value of the state variable by prepending `D` to the variable name.
+
+We create a new `pomp` object based on this process model but with the same rmeasure and dmeasure components as in `parus` by executing
+```{r bh-pomp1,cache=FALSE}
+parus.bh <- pomp(parus,rprocess=discrete.time.sim(bh.step,delta.t=1),skeleton=bh.skel,skeleton.type="map",skelmap.delta.t=1,statenames="N",paramnames=c("a","b","sigma"))
+```
+The following codes compute some simulations, a trajectory of the skeleton, and run a particle filtering operation.
+```{r bh-test}
+coef(parus.bh) <- c(a=1.1,b=5e-4,sigma=0.5,N.0=30,phi=1)
+sim <- simulate(parus.bh)
+traj <- trajectory(parus.bh)
+pf <- pfilter(parus.bh,Np=1000)
+```
+
+## Parameter transformations
+
+We can specify model-specific parameter transformations using C snippets.
+The following implements log transformation of the $r$ and $K$ parameters of the logistic model and folds them into the `pomp` object `parus`.
+```{r logistic-partrans,cache=FALSE}
+partrans <- Csnippet("
+ Tr = exp(r);
+ TK = exp(K);
+ ")
+
+parinvtrans <- Csnippet("
+ Tr = log(r);
+ TK = log(K);
+ ")
+
+parus <- pomp(parus,parameter.transform=partrans,parameter.inv.transform=parinvtrans,
+ paramnames=c("r","K"))
+```
+```{r logistic-partrans-test,include=FALSE,cache=FALSE}
+p <- c(r=1,K=200,phi=1,N.0=200,sigma=0.5)
+coef(parus,transform=TRUE) <- partrans(parus,p,dir="inv")
+stopifnot(all.equal(p,coef(parus)))
+```
+
+## Rules for writing C snippets
+
+To fully understand how `pomp` makes use of C snippets, it is informative to have a look at the C file created by any call to `pomp` involving one or more C snippets.
+Run the following in an `R` session and view the C file that is created in the session's temporary directory.
+```{r tempfile-view,eval=FALSE}
+op <- options(verbose=TRUE)
+
+parus <- pomp(data=parus.dat,times="year",t0=1959,rprocess=euler.sim(step.fun,delta.t=1/365),
+ dmeasure=dmeas,rmeasure=rmeas,statenames="N",paramnames=c("r","phi","K","N.0","sigma"))
+
+options(op)
+```
+
+### General rules
+
+1. C snippets must be valid C. They will embedded verbatim in a template file which will then be compiled by a call to `R CMD SHLIB`. If the resulting file does not compile, an error message wil be generated. No attempt is made by `pomp` to interpret this message. Typically, compilation errors are due to either the snippet's being invalid C or one or more undeclared variables.
+1. Each state variable or parameter must be left undeclared in the snippet and declared in the corresponding argument to `pomp`, i.e., in either `statenames` or `paramnames`. Compiler errors that complain about undeclared state variables or parameters are likely due to failure to include these parameters in the appropriate vector of names.
+1. Names of observables are determined by their names in the data. They must be referred to in measurement model snippets (rmeasure or dmeasure) by those names.
+1. If the `pomp` object contains a table of covariates (see `?pomp`), then the variables in the covariate table will be available by their names in the context within which the snippet is executed.
+1. `R` variables with names containing dots are replaced in the C codes by variable names in which all dots have been replaced by underscores ("_").
+1. A C snippet can declare local variables. Be careful not to use names that match those of state variables, observables, or parameters. The latter must never be declared within a snippet.
+1. The header `R.h`, provided with `R`, will be included, making all of the [`R` C API](http://cran.r-project.org/doc/manuals/r-release/R-exts.html#The-R-API) available for use in the snippet.
+1. The header `pomp.h`, provided with `pomp`, will be included, making all of the `pomp` C API available for use in every snippet. Do `file.show(system.file("include/pomp.h",package="pomp"))` to view this header file.
+1. Snippets of C code passed to the `globals` argument of `pomp` will be included at the head of the generated C file. This can be used to define global variables, define useful functions, and include header files.
+
+### C snippets implementing Euler steps for use by `euler.sim`
+
+1. The ultimate goal of such a snippet is to replace the state variables with their new random values at the end of the time interval. Accordingly, each state variable should be over-written with its new value.
+1. In addition to the states, parameters, covariates (if any), and observables, the variables `t` and `dt`, containing respectively the time at the beginning of the Euler step and the Euler step's duration, will be defined in the context in which the snippet is executed.
+
+### C snippets implementing the rmeasure component
+
+1. The ultimate goal of such a snippet is to fill the observables with random values drawn from the measurement model distribution. Accordingly, each observable should be assigned a new value.
+1. In addition to the states, parameters, covariates (if any), and observables, the variable `t`, containing the time of the observation, will be defined in the context in which the snippet is executed.
+
+### C snippets implementing the dmeasure component
+
+1. The ultimate goal of such a snippet is to fill the `lik` variable with the likelihood of the data given the state. Alternatively, if `give_log=1`, `lik` should be filled with the log likelihood.
+1. In addition to the states, parameters, covariates (if any), and observables, the variable `t`, containing the time of the observation, will be defined in the context in which the snippet is executed.
+
+### C snippets implementing the deterministic skeleton
+
+1. For each state variable, there is a corresponding component of the deterministic skeleton. The goal of such a snippet is to compute all the components.
+1. When the skeleton is a map, the component corresponding to state variable `x` is named `Dx` and is the new value of `x` after one iteration of the map.
+1. When the skeleton is a vectorfield, the component corresponding to state variable `x` is named `Dx` and is the value of $dx/dt$.
+1. As with the other C snippets, all states, parameters and covariates, as well as the current time, `t`, will be defined in the context within which the snippet is executed.
+
+### C snippets implementing parameter transformations
+
+1. The parameter transformation taking a parameter vector from the scale defined by the model codes to another scale is specified as the `parameter.inv.transform` whilst the transformation taking a parameter vector from the alternative scale to that on which the model is defined is specified as the `parameter.transform`.
+1. The goal of a either of these snippets is the computation of the values of the transformed parameters. The value of transformed parameter `x` should be assigned to variable `Tx`.
+1. Time- and covariate-dependent transformations are not allowed. Therefore, the variables `t` and any covariates will not be available in the context within which a parameter transformation snippet is executed.
+
+## References
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 1120
More information about the pomp-commits
mailing list