[Pomp-commits] r246 - in pkg: . R data inst/doc man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue May 18 20:36:34 CEST 2010
Author: kingaa
Date: 2010-05-18 20:36:33 +0200 (Tue, 18 May 2010)
New Revision: 246
Modified:
pkg/NAMESPACE
pkg/R/plugins.R
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/ou2.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.Rnw
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/doc/ou2-first-mif.rda
pkg/inst/doc/ou2-trajmatch.rda
pkg/man/plugins.Rd
Log:
- add the rprocess plugin 'discrete.time.sim', which is essentially identical to 'euler.sim'
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-05-17 21:48:41 UTC (rev 245)
+++ pkg/NAMESPACE 2010-05-18 18:36:33 UTC (rev 246)
@@ -41,6 +41,7 @@
onestep.simulate,
onestep.density,
euler.sim,
+ discrete.time.sim,
onestep.sim,
onestep.dens,
gillespie.sim,
Modified: pkg/R/plugins.R
===================================================================
--- pkg/R/plugins.R 2010-05-17 21:48:41 UTC (rev 245)
+++ pkg/R/plugins.R 2010-05-18 18:36:33 UTC (rev 246)
@@ -29,6 +29,37 @@
}
}
+discrete.time.sim <- function (step.fun, delta.t = 1, PACKAGE) {
+ efun <- pomp.fun(
+ f=step.fun,
+ PACKAGE=PACKAGE,
+ proto="step.fun(x,t,params,delta.t,...)"
+ )
+ function (xstart, times, params, ...,
+ statenames = character(0),
+ paramnames = character(0),
+ covarnames = character(0),
+ zeronames = character(0),
+ tcovar, covar) {
+ .Call(
+ euler_model_simulator,
+ func=efun,
+ xstart=xstart,
+ times=times,
+ params=params,
+ dt=delta.t,
+ method=0L,
+ statenames=statenames,
+ paramnames=paramnames,
+ covarnames=covarnames,
+ zeronames=zeronames,
+ tcovar=tcovar,
+ covar=covar,
+ args=pairlist(...)
+ )
+ }
+}
+
euler.sim <- function (step.fun, delta.t, PACKAGE) {
efun <- pomp.fun(
f=step.fun,
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/verhulst.rda
===================================================================
(Binary files differ)
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 2010-05-17 21:48:41 UTC (rev 245)
+++ pkg/inst/doc/intro_to_pomp.Rnw 2010-05-18 18:36:33 UTC (rev 246)
@@ -220,6 +220,7 @@
When this function is called, the unobserved states at time \code{t} will be in the named numeric vector \code{x} and the parameters in \code{params} as before.
The function returns a named numeric vector that represents a single draw from the observation process \eqref{eq:ou-obs}.
+\clearpage
\section{Simulating the model}
With the two functions above, we already have all we need to simulate the entire model.
@@ -233,7 +234,7 @@
y2=NA
),
times="time",
- rprocess=euler.sim(
+ rprocess=discrete.time.sim(
step.fun=ou2.proc.sim,
delta.t=1
),
@@ -246,7 +247,7 @@
In a moment, however, we'll simulate some data so we can explore \pomp's various methods for fitting models to data.
The second argument (\code{times}) specifies which of the columns of \code{data} is the time variable.
The third argument (\code{rprocess}) specifies that the process model simulator will be in discrete-time, one step at a time.
-The function \code{euler.sim} belongs to the package \pomp.
+The function \code{discrete.time.sim} belongs to the package \pomp.
It takes several arguments, one of which, \code{step.fun}, is the particular function that actually takes the step.
Its second argument, \code{delta.t}, sets the size of the time-step.
The argument \code{rmeasure} specifies the measurement model simulator function.
@@ -316,6 +317,8 @@
\label{fig:ou2-first-simulation-plot}
\end{figure}
+\section{Computing likelihood using particle filtering}
+
Since some parameter estimation algorithms in the \pomp\ package only require simulations of the full process, we are already in a position to use them.
At present, the only such algorithm is nonlinear forecasting (\code{nlf}).
In the near future, probe-matching algorithms that require only \code{rprocess} and \code{rmeasure} will be included in the package.
@@ -346,7 +349,7 @@
ou2 <- pomp(
data=dat[c("time","y1","y2")],
times="time",
- rprocess=euler.sim(ou2.proc.sim,delta.t=1),
+ rprocess=discrete.time.sim(ou2.proc.sim,delta.t=1),
rmeasure=ou2.meas.sim,
dmeasure=ou2.meas.dens,
t0=0
@@ -360,8 +363,6 @@
Note that we've first extracted the data from our old \code{ou2} and set up the new one with the same data and parameters.
The calls to \code{coef} and \code{coef<-} in the lines above make sure the parameters have the same values they had before.
-\section{Computing likelihood using particle filtering}
-
To compute the likelihood of the data, we can use the function \code{pfilter}.
This runs a plain vanilla particle filter (AKA sequential Monte Carlo) algorithm and results in an unbiased estimate of the likelihood.
See \citet{Arulampalam2002} for an excellent tutorial on particle filtering.
@@ -690,7 +691,7 @@
new.ou2 <- pomp(
data=dat[c("time","y1","y2")],
times="time",
- rprocess=euler.sim(ou2.proc.sim,delta.t=1),
+ rprocess=discrete.time.sim(ou2.proc.sim,delta.t=1),
rmeasure=ou2.meas.sim,
dmeasure=ou2.meas.dens,
skeleton.map=ou2.skel,
@@ -849,7 +850,8 @@
@
The specification of \code{data}, \code{times}, and \code{t0} should be familiar.
The covariates are specified using the arguments \code{tcovar} and \code{covar}.
-Again, we use \code{euler.sim} to specify the process simulator (\code{rprocess}): we use a time-step of 1/20 of a week.
+We use \code{euler.sim} to specify the process simulator (\code{rprocess}).
+We are approximating the continuous-time process using an Euler simulator with a time-step of 1/20 of a week.
Both \code{rmeasure} and \code{dmeasure} can be specified at once using the \code{measurement.model} argument.
Here, we model the observation process using a binomial process, where the \emph{reporting rate}, \code{rho}, is the probability that a case is reported.
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ou2-first-mif.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ou2-trajmatch.rda
===================================================================
(Binary files differ)
Modified: pkg/man/plugins.Rd
===================================================================
--- pkg/man/plugins.Rd 2010-05-17 21:48:41 UTC (rev 245)
+++ pkg/man/plugins.Rd 2010-05-18 18:36:33 UTC (rev 246)
@@ -2,6 +2,7 @@
\alias{plugins}
\alias{onestep.sim}
\alias{euler.sim}
+\alias{discrete.time.sim}
\alias{onestep.dens}
\alias{gillespie.sim}
\title{Plug-ins for dynamical models based on stochastic Euler algorithms}
@@ -12,6 +13,7 @@
\usage{
onestep.sim(step.fun, PACKAGE)
euler.sim(step.fun, delta.t, PACKAGE)
+discrete.time.sim(step.fun, delta.t = 1, PACKAGE)
gillespie.sim(rate.fun, v, d, PACKAGE)
onestep.dens(dens.fun, PACKAGE)
}
@@ -61,21 +63,23 @@
}
}
\details{
- \code{onestep.sim} is the appropriate choice when it is possible to simulate the change in state from one time to another, regardless of how large the difference in times is.
+ \code{onestep.sim} is the appropriate choice when it is possible to simulate the change in state from one time to another, regardless of how large the interval between them is.
To use \code{onestep.sim}, you must write a function \code{step.fun} that will advance the state process from one arbitrary time to another.
- \code{euler.sim} is appropriate when one cannot do this but can compute the change in state via repeated simulation of smaller steps.
- There are two principal circumstances under which the latter approach is desirable.
- First, if one is simulating a continuous time process but is willing to approximate it using an Euler approach.
- Second, if one is simulating a discrete-time process.
- To use \code{euler.sim}, you must write a function \code{step.fun} that will take a single Euler step, of size at most \code{delta.t}.
- \code{euler.sim} will take as many Euler steps as it needs to in order to get from one time to another.
+ \code{euler.sim} is appropriate when one cannot do this but can compute the change in state via a sequence of smaller steps.
+ This is desirable, for example, if one is simulating a continuous time process but is willing to approximate it using an Euler approach.
+ \code{discrete.time.sim} is appropriate when the process evolves in discrete time.
+ In this case, it is assumed that the intervals between observations are integer multiples of \code{delta.t}.
+ Beware, however: this assumption is not checked.
+
+ To use \code{euler.sim} or \code{discrete.time.sim}, you must write a function \code{step.fun} that will take a single Euler step, of size at most \code{delta.t}.
+ \code{euler.sim} and \code{discrete.time.sim} will create simulators that take as many steps as needed to get from one time to another.
See below for information on how \code{euler.sim} chooses the actual step size it uses.
- If one uses \code{euler.sim} to simulate a discrete-time process, \code{delta.t} should be set to the discrete time-step.
+
\code{gillespie.sim} allows exact simulation of a continuous-time, discrete-state Markov process using Gillespie's algorithm.
- This is an \dQuote{event-driven} approach: correspondingly, to use \code{gillespie.sim}, you must write a function \code{rate.fun} that computes the rates of each elementary event and specify two matrices that describe the dependencies of each rate and the consequences of each event.
+ This is an \dQuote{event-driven} approach: correspondingly, to use \code{gillespie.sim}, you must write a function \code{rate.fun} that computes the rates of each elementary event and specify two matrices (\code{d,v}) that describe, respectively, the dependencies of each rate and the consequences of each event.
\code{onestep.dens} will generate a suitable \code{dprocess} function when one can compute the likelihood of a given state transition simply by knowing the states at two times under the assumption that the state has not changed between the times.
- This is typically possible, for instance, when the \code{rprocess} function is implemented using \code{euler.sim}.
+ This is typically possible, for instance, when the \code{rprocess} function is implemented using \code{onestep.sim}, \code{euler.sim}, or \code{discrete.time.sim}.
[NB: currently, there are no high-level algorithms in \pkg{pomp} that use \code{dprocess}.
This function is provided for completeness only, and with an eye toward future development.]
@@ -85,7 +89,7 @@
If the argument \code{covars} is included and a covariate table has been included in the \code{pomp} object, then on a call to this function, \code{covars} will be filled with the values, at time \code{t}, of the covariates.
This is accomplished via interpolation of the covariate table.
Additional arguments may be given: these will be filled by the correspondingly-named elements in the \code{userdata} slot of the \code{pomp} object (see \code{\link{pomp}}).
- If \code{step.fun} is written in a native language, it must be a function of type "pomp_onestep_sim" as specified in the header "pomp.h" included with the package (see the directory "include" in the installed package directory).
+ If \code{step.fun} is written in a native language, it must be a function of type \dQuote{pomp_onestep_sim} as specified in the header \dQuote{pomp.h} included with the package (see the directory \dQuote{include} in the installed package directory).
If \code{rate.fun} is written as an \R function, it must have at least the arguments \code{j}, \code{x}, \code{t}, \code{params}, and \code{\dots}.
Here, \code{j} is the an integer that indicates which specific elementary event we desire the rate of.
@@ -95,8 +99,7 @@
This is accomplished via interpolation of the covariate table.
If \code{rate.fun} is a native function, it must be of type \dQuote{pomp_ssa_rate_fn} as defined in the header \dQuote{pomp.h}, which is included with the package.
- In writing \code{onestep.dens}, you must assume that no state transitions have occurred between \code{t1} and \code{t2}.
-
+ In writing \code{dens.fun}, you must assume that no state transitions have occurred between \code{t1} and \code{t2}.
If \code{dens.fun} is written as an \R function, it must have at least the arguments \code{x1}, \code{x2}, \code{t1}, \code{t2}, \code{params}, and \code{\dots}.
On a call to this function, \code{x1} and \code{x2} will be named vectors of state variables at times \code{t1} and \code{t2}, respectively.
The named vector \code{params} contains the parameters.
@@ -104,10 +107,10 @@
If the argument \code{covars} is included and a covariate table has been included in the \code{pomp} object, then on a call to this function, \code{covars} will be filled with the values, at time \code{t1}, of the covariates.
This is accomplished via interpolation of the covariate table.
As above, any additional arguments will be filled by the correspondingly-named elements in the \code{userdata} slot of the \code{pomp} object (see \code{\link{pomp}}).
- If \code{dens.fun} is written in a native language, it must be a function of type "pomp_onestep_pdf" as defined in the header "pomp.h" included with the package (see the directory "include" in the installed package directory).
+ If \code{dens.fun} is written in a native language, it must be a function of type \dQuote{pomp_onestep_pdf} as defined in the header \dQuote{pomp.h} included with the package (see the directory \dQuote{include} in the installed package directory).
}
\value{
- \code{onestep.sim}, \code{euler.sim}, and \code{gillespie.sim} each return functions suitable for use as the argument \code{rprocess} argument in \code{\link{pomp}}.
+ \code{onestep.sim}, \code{euler.sim}, \code{discrete.time.sim}, and \code{gillespie.sim} each return functions suitable for use as the argument \code{rprocess} argument in \code{\link{pomp}}.
\code{onestep.dens} returns a function suitable for use as the argument \code{dprocess} in \code{\link{pomp}}.
}
More information about the pomp-commits
mailing list