[Pomp-commits] r224 - in pkg: inst/include man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 10 15:51:49 CEST 2010
Author: kingaa
Date: 2010-05-10 15:51:49 +0200 (Mon, 10 May 2010)
New Revision: 224
Modified:
pkg/inst/include/pomp.h
pkg/man/euler.sir.Rd
pkg/man/ou2.Rd
pkg/man/pomp-package.Rd
pkg/man/pomp.Rd
pkg/man/verhulst.Rd
pkg/tests/logistic.R
pkg/tests/logistic.Rout.save
pkg/tests/ou2-kalman.R
pkg/tests/ou2-kalman.Rout.save
pkg/tests/ou2-mif.R
pkg/tests/ou2-mif.Rout.save
pkg/tests/ou2-nlf.R
pkg/tests/ou2-nlf.Rout.save
pkg/tests/ou2-trajmatch.R
pkg/tests/ou2-trajmatch.Rout.save
pkg/tests/pfilter.Rout.save
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
- tests and man pages corresponding to last changes
Modified: pkg/inst/include/pomp.h
===================================================================
--- pkg/inst/include/pomp.h 2010-05-10 13:50:37 UTC (rev 223)
+++ pkg/inst/include/pomp.h 2010-05-10 13:51:49 UTC (rev 224)
@@ -18,7 +18,28 @@
// facility for computing evaluating a basis of periodic bsplines
void periodic_bspline_basis_eval (double x, double period, int degree, int nbasis, double *y);
-// Prototype for one-step simulator, as used by "euler.simulate" and "onestep.simulate":
+// Prototype for stochastic simulation algorithm reaction-rate function, as used by "gillespie.sim":
+typedef double pomp_ssa_rate_fn(int j, double t, const double *x, const double *p,
+ int *stateindex, int *parindex, int *covindex,
+ int ncovar, double *covars);
+// Description:
+// on input:
+// j = integer specifying the number of the reaction whose rate is desired
+// t = time at which the rates are to be evaluated
+// x = vector of state variables
+// p = vector of parameters
+// stateindex = pointer to vector of integers pointing to the states in 'x' in the order specified by
+// the 'statenames' argument of 'SSA.simulator'
+// parindex = pointer to vector of integers pointing to the parameters in 'p' in the order specified by
+// the 'paramnames' argument of 'SSA.simulator'
+// covindex = pointer to vector of integers pointing to the covariates in 'covars' in the order
+// specified by the 'covarnames' argument of 'SSA.simulator'
+// ncovars = number of covariates
+// covars = pointer to a vector containing the values of the covariates at time t, as interpolated
+// from the covariate table supplied to 'SSA.simulator'
+// returns the rate of the j-th reaction
+
+// Prototype for one-step simulator, as used by "euler.sim" and "onestep.sim":
typedef void pomp_onestep_sim(double *x, const double *p,
const int *stateindex, const int *parindex, const int *covindex,
int ncovars, const double *covars,
@@ -46,7 +67,7 @@
// Inclusion of these calls in the user-defined function may result in significant slowdown.
-// Prototype for one-step log probability density function, as used by "onestep.density":
+// Prototype for one-step log probability density function, as used by "onestep.dens":
typedef void pomp_onestep_pdf(double *f,
double *x1, double *x2, double t1, double t2, const double *p,
const int *stateindex, const int *parindex, const int *covindex,
Modified: pkg/man/euler.sir.Rd
===================================================================
--- pkg/man/euler.sir.Rd 2010-05-10 13:50:37 UTC (rev 223)
+++ pkg/man/euler.sir.Rd 2010-05-10 13:51:49 UTC (rev 224)
@@ -1,16 +1,27 @@
-\name{euler.sir}
+\name{sir}
\alias{euler.sir}
+\alias{gillespie.sir}
\docType{data}
-\title{Seasonal SIR model implemented as an Euler-multinomial model}
+\title{Seasonal SIR model implemented using two stochastic simulation algorithms.}
\description{
\code{euler.sir} is a \code{pomp} object encoding a simple seasonal SIR model.
+ Simulation is performed using an Euler multinomial approximation (AKA tau leap method).
+ \code{gillespie.sir} has the same model implemented using Gillespie's algorithm.
}
-\usage{data(euler.sir)}
+\usage{
+data(euler.sir)
+data(gillespie.sir)
+}
\examples{
data(euler.sir)
plot(euler.sir)
x <- simulate(euler.sir,nsim=10,seed=20348585)
plot(x[[1]])
+
+data(gillespie.sir)
+plot(gillespie.sir)
+x <- simulate(gillespie.sir,nsim=1,seed=20348585)
+plot(x)
}
\seealso{\code{\link{pomp-class}} and the vignettes}
\keyword{datasets}
Modified: pkg/man/ou2.Rd
===================================================================
--- pkg/man/ou2.Rd 2010-05-10 13:50:37 UTC (rev 223)
+++ pkg/man/ou2.Rd 2010-05-10 13:51:49 UTC (rev 224)
@@ -3,28 +3,28 @@
\docType{data}
\title{Two-dimensional Ornstein-Uhlenbeck process}
\description{
- \code{ou2} is a \code{pomp} object encoding a bivariate Ornstein-Uhlenbeck process.
+ \code{ou2} is a \code{pomp} object encoding a bivariate discrete-time Ornstein-Uhlenbeck process.
}
\usage{data(ou2)}
\details{
- The OU process is fully but noisily observed.
- The functions \code{rprocess}, \code{dprocess}, \code{rmeasure}, and \code{dmeasure} are implemented using compiled C code for computational speed:
+ The process is fully but noisily observed.
+ The functions \code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, and \code{skeleton} are implemented using compiled C code for computational speed:
see the source code for details.
- The use of this object is demonstrated in the vignette.
+ The use of this object is demonstrated in the vignette "intro_to_pomp".
}
\examples{
data(ou2)
plot(ou2)
coef(ou2)
-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
- )
-x <- simulate(ou2,nsim=10,seed=20348585,params=p)
+theta <- c(
+ 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,
+ tau=1,
+ x1.0=-3, x2.0=4
+ )
+x <- simulate(ou2,nsim=10,seed=20348585,params=theta)
plot(x[[1]])
-pfilter(ou2,params=p,Np=1000)$loglik
+pfilter(ou2,params=theta,Np=1000)$loglik
}
-\seealso{\code{\link{pomp-class}} and the vignettes}
+\seealso{\code{\link{pomp}} and the vignettes}
\keyword{datasets}
Modified: pkg/man/pomp-package.Rd
===================================================================
--- pkg/man/pomp-package.Rd 2010-05-10 13:50:37 UTC (rev 223)
+++ pkg/man/pomp-package.Rd 2010-05-10 13:51:49 UTC (rev 224)
@@ -3,23 +3,29 @@
\alias{pomp-package}
\title{Partially-observed Markov processes}
\description{
- This is documentation for version \Sexpr[echo=F,stage=install,results=text]{packageDescription("pomp",fields="Version")} of \pkg{pomp}.
+%% This is documentation for version \Sexpr[echo=F,stage=install,results=text]{packageDescription("pomp",fields="Version")} of \pkg{pomp}.
- The \code{pomp} package provides facilities for inference on time series data using partially-observed Markov processes (AKA state-space models or nonlinear stochastic dynamical systems).
- The user encodes a model as a \code{\link{pomp}} object by providing functions specifying some or all of the model's process and measurement components.
- The package's algorithms for fitting models to data, simulating, etc. then call these functions.
- At the moment, algorithms are provided for
- particle filtering (AKA sequential Monte Carlo or sequential importance sampling, see \code{\link{pfilter}}),
- the likelihood maximization by iterated filtering (\acronym{MIF}) method of Ionides, Breto, and King (PNAS, 103:18438-18443, 2006, see \code{\link{mif}}),
- and the nonlinear forecasting algorithm of Kendall, Ellner, et al. (Ecol. Monog. 75:259-276, 2005, see \code{\link{nlf}}).
- Future support for a variety of other algorithms is envisioned.
- A working group of the National Center for Ecological Analysis and Synthesis (\acronym{NCEAS}), "Inference for Mechanistic Models", is currently implementing and testing additional methods for this package.
+ The \pkg{pomp} package provides facilities for inference on time series data using partially-observed Markov processes (AKA state-space models or nonlinear stochastic dynamical systems).
+ One might use \pkg{pomp} to fit a model to time-series data or simply to simulate a stochastic dynamical model.
+ The first step in using \pkg{pomp} is to encode one's model and data in an object of class \code{pomp}.
+ One does this via a call to \code{\link{pomp}}, which involves specifying the process and measurement components of the model in one or more of a variety of ways.
+ Details on this are given in the documentation for the \code{\link{pomp}} function and examples are given in the \sQuote{intro_to_pomp} vignette.
+ Currently, \pkg{pomp} provides algorithms for
+ (i) simulation of stochastic dynamical systems (see \code{\link[=simulate-pomp]{simulate}}),
+ (ii) particle filtering (AKA sequential Monte Carlo or sequential importance sampling, see \code{\link{pfilter}}),
+ (iii) the likelihood maximization by iterated filtering (\acronym{MIF}) method of Ionides, Breto, and King (PNAS, 103:18438-18443, 2006, see \code{\link{mif}}),
+ and (iv) the nonlinear forecasting algorithm of Kendall, Ellner, et al. (Ecol. Monog. 75:259-276, 2005, see \code{\link{nlf}}).
+ Future support for additional algorithms in envisioned, and implementations of
+ the Bayesian sequential Monte Carlo approach of Liu \& West and generalized probe-matching methods are currently under development.
+ Much of the work in \pkg{pomp} has been done under the auspices of a working group of the National Center for Ecological Analysis and Synthesis (\acronym{NCEAS}), "Inference for Mechanistic Models".
+
The package is provided under the \acronym{GNU} Public License (\acronym{GPL}).
Contributions are welcome, as are comments, suggestions for improvements, and bug reports.
}
\section{Classes}{
- The basic class, \code{\link{pomp}}, is provided to encode a partially-observed Markov process together with a uni- or multi-variate data set.
+ \pkg{pomp} makes extensive use of S4 classes.
+ The basic class, \code{\link{pomp}}, encodes a partially-observed Markov process together with a uni- or multi-variate data set and (possibly) parameters.
}
\section{Vignettes}{
The vignette \sQuote{intro_to_pomp} illustrates the facilities of the package using familiar stochastic processes.
@@ -27,7 +33,11 @@
}
\author{Aaron A. King \email{kingaa at umich dot edu}}
\seealso{
- \code{\link{pomp}}, \code{\link{simulate-pomp}}, \code{\link{mif}}, \code{\link{nlf}}
+ \code{\link{pomp}},
+ \code{\link{pfilter}},
+ \code{\link[=simulate-pomp]{simulate}},
+ \code{\link{mif}},
+ \code{\link{nlf}}
}
\keyword{models}
\keyword{ts}
Modified: pkg/man/pomp.Rd
===================================================================
--- pkg/man/pomp.Rd 2010-05-10 13:50:37 UTC (rev 223)
+++ pkg/man/pomp.Rd 2010-05-10 13:51:49 UTC (rev 224)
@@ -14,7 +14,7 @@
\item{data}{
An array holding the data.
This array is of dimensions \code{nobs} x \code{ntimes}, where \code{nobs} is the number of observed variables and \code{ntimes} is the number of times at which observations were made.
- It is also possible to specify \code{data} as a data-frame, in which case the argument \code{times} must be the name of the column of observation times.
+ One can also specify \code{data} as a data-frame, in which case the argument \code{times} must be the name of the column of observation times.
Note that if \code{data} is a data-frame, it will be internally coerced to an array with storage-mode \sQuote{double}.
}
\item{times}{
@@ -25,60 +25,62 @@
\item{t0}{
The zero-time.
This must be no later than the time of the first observation, \code{times[1]}.
+ The stochastic dynamical system is initialized at time \code{t0}.
}
\item{rprocess}{
optional; a function of prototype \code{rprocess(xstart,times,params,\dots)} that simulates from the unobserved process.
+ The easiest way to specify \code{rprocess} is to use one of the \code{\link{plugins}} provided as part of the \pkg{pomp} package.
See below for details.
}
\item{dprocess}{
optional; a function of prototype \code{dprocess(x,times,params,log,\dots)} that evaluates the likelihood of a sequence of consecutive state transitions.
+ The easiest way to specify \code{dprocess} is to use one of the \code{\link{plugins}} provided as part of the \pkg{pomp} package.
+ It is not typically necessary (or even feasible) to define \code{dprocess}.
See below for details.
}
\item{rmeasure}{
- optional; the measurement model simulator.
+ optional function; the measurement model simulator.
This can be specified in one of three ways:
(1) as a function of prototype \code{rmeasure(x,t,params,\dots)} that makes a draw from the observation process given states \code{x}, time \code{t}, and parameters \code{params}.
- (2) as the name of a native (compiled) routine with prototype \dQuote{pomp_measure_model_simulator} as defined in the header file \dQuote{pomp.h}.
+ (2) as the name of a native (compiled) routine with prototype \dQuote{pomp_measure_model_simulator} as defined in the header file \dQuote{examples/pomp.h}.
In the above cases, if the measurement model depends on covariates, the optional argument \code{covars} will be filled with interpolated values at each call.
(3) using the formula-based \code{measurement.model} facility (see below).
}
\item{dmeasure}{
- optional; the measurement model probability density function.
+ optional function; the measurement model probability density function.
This can be specified in one of three ways:
(1) as a function of prototype \code{dmeasure(y,x,t,params,log,\dots)} that computes the p.d.f. of \code{y} given \code{x}, \code{t}, and \code{params}.
- (2) as the name of a native (compiled) routine with prototype \dQuote{pomp_measure_model_density} as defined in the header file \dQuote{pomp.h}.
+ (2) as the name of a native (compiled) routine with prototype \dQuote{pomp_measure_model_density} as defined in the header file \dQuote{examples/pomp.h}.
In the above cases, if the measurement model depends on covariates, the optional argument \code{covars} will be filled with interpolated values at each call.
(3) using the formula-based \code{measurement.model} facility (see below).
As might be expected, if \code{log=TRUE}, this function should return the log likelihood.
}
\item{measurement.model}{
- optional; a formula or list of formulae, specifying the measurement model.
- These formulae are parsed and used to generate \code{rmeasure} and \code{dmeasure} functions.
+ optional formula; a formula or list of formulae, specifying the measurement model.
+ These formulae are parsed internally and used to generate \code{rmeasure} and \code{dmeasure} functions.
If \code{measurement.model} is given it overrides any specification of \code{rmeasure} or \code{dmeasure}.
See below for an example.
\strong{NB:} it will typically be possible to acclerate measurement model computations by writing \code{dmeasure} and/or \code{rmeasure} functions directly.
}
- \item{skeleton.map}{
- optional.
- If we are dealing with a discrete-time Markov process, its deterministic skeleton is a map and can be specified using \code{skeleton.map} in one of two ways:
+ \item{skeleton.map, skeleton.vectorfield}{
+ optional functions.
+ If we are dealing with a discrete-time Markov process, its deterministic skeleton is a map and can be specified using \code{skeleton.map}.
+ If we are dealing with a continuous-time Markov process, its deterministic skeleton is a vectorfield and can be specified using \code{skeleton.vectorfield}.
+ In either event, the skeleton function can be specified in one of two ways:
(1) as a function of prototype \code{skeleton(x,t,params,\dots)} that evaluates the deterministic skeleton at state \code{x} and time \code{t} given the parameters \code{params}, or
(2) as the name of a native (compiled) routine with prototype \dQuote{pomp_vectorfield_map} as defined in the header file \dQuote{pomp.h}.
If the deterministic skeleton depends on covariates, the optional argument \code{covars} will be filled with interpolated values of the covariates at the time \code{t}.
+ A dynamical system cannot be both discrete and continuous, so specifying both \code{skeleton.map} and \code{skeleton.vectorfield} generates an error.
}
- \item{skeleton.vectorfield}{
- optional.
- If we are dealing with a continuous-time Markov process, its deterministic skeleton is a vectorfield and can be specified using \code{skeleton.vectorfield} in either of the two ways described above, under \code{skeleton.map}.
- Note that it makes no sense to specify the deterministic skeleton both as a map and as a vectorfield: an attempt to do so will generate an error.
- }
\item{initializer}{
- optional; a function of prototype \code{initializer(params,t0,\dots)} that yields initial conditions for the state process when given a vector, \code{params}, of parameters.
- By default, i.e., if it is unspecified when \code{pomp} is called, the initializer assumes any names in \code{params} that end in \dQuote{.0} correspond to initial value parameters.
+ optional function of prototype \code{initializer(params,t0,\dots)} that yields initial conditions for the state process when given a vector, \code{params}, of parameters.
+ By default (i.e., if it is unspecified when \code{pomp} is called), the initializer assumes any parameters in \code{params} the names of which end in \dQuote{\code{.0}} are initial values.
These are simply copied over as initial conditions when \code{init.state} is called (see \code{\link{init.state-pomp}}).
- The names of the state variables are the same as the corresponding initial value parameters, but with the \dQuote{.0} dropped.
+ The names of the state variables are the same as the corresponding initial value parameters, but with the \dQuote{\code{.0}} dropped.
}
\item{covar, tcovar}{
An optional table of covariates: \code{covar} is the table (with one column per variable) and \code{tcovar} the corresponding times (one entry per row of \code{covar}).
- This can be in the form of a matrix or a data frame.
+ \code{covar} can be specified as either a matrix or a data frame.
In either case the columns are taken to be distinct covariates.
If \code{covar} is a data frame, \code{tcovar} can be either the name or the index of the time variable.
If a covariate table is supplied, then the value of each of the covariates is interpolated as needed, i.e., whenever \code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, or \code{init.state} is evaluated.
@@ -88,9 +90,10 @@
Optional character vectors specifying the names of state variables, parameters, or covariates, respectively.
These are only used in the event that one or more of the basic functions (\code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, \code{skeleton}) are defined using native routines.
In that case, these name vectors are matched against the corresponding names and the indices of the names are passed to the native routines.
+ Using this facility allows one to write one or more of \code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, \code{skeleton} in native code in a way that does not depend on the order of states, parameters, and covariates at run time.
}
\item{PACKAGE}{
- An optional string giving the name of the dynamically loaded library in which the native routines are to be found.
+ An optional string giving the name of the dynamically loaded library in which any native routines are to be found.
}
\item{\dots}{
Any additional arguments are passed as arguments to each of the functions \code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, and \code{initializer} whenever they are evaluated.
@@ -98,14 +101,17 @@
}
\details{
\strong{
- It is not typically necessary (or easy) to define all of the functions \code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, and \code{initializer} in any given problem.
+ It is not typically necessary (or desirable, or even feasible) to define all of the functions \code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, and \code{skeleton} in any given problem.
Each algorithm makes use of a different subset of these functions.
}
- The specification of process-model codes \code{rprocess} and/or \code{dprocess} can be somewhat nontrivial:
- for this reason, \emph{plug-ins} have been developed to streamline this process for the user.
- At present, if the process model evolves in discrete-time or one is willing to make such an approximation (e.g., via an Euler approximation), then the \code{\link{onestep.simulate}} and \code{\link{euler.simulate}} plug-ins for \code{rprocess} and \code{\link{onestep.density}} plug-in for \code{dprocess} are available.
- The following describes how each of these functions should be written.
- Examples are provided in the help documentation and the vignettes.
+ In general, the specification of process-model codes \code{rprocess} and/or \code{dprocess} can be somewhat nontrivial:
+ for this reason, \code{\link{plugins}} have been developed to streamline this process for the user.
+ Currently, if one's process model evolves in discrete time or one is willing to make such an approximation (e.g., via an Euler approximation), then the \code{\link{euler.sim}} or \code{\link{onestep.sim}} plugin for \code{rprocess} and \code{\link{onestep.dens}} plugin for \code{dprocess} are available.
+ For exact simulation of certain continuous-time Markov chains, an implementation of Gillespie's algorithm is available (see \code{\link{gillespie.sim}}).
+ To use the plugins, consult the help documentation (\code{?\link{plugins}}) and the vignettes.
+
+ It is anticipated that in specific cases, it may be possible to obtain increased computational efficiency by writing custom versions of \code{rprocess} and/or \code{dprocess}.
+ The following describes how each of these functions should be written in this case.
\describe{
\item{\code{rprocess}}{
if provided, must have at least the following arguments:
@@ -126,8 +132,6 @@
In particular, \code{x[,,1]} must be identical to \code{xstart}.
The rownames of \code{x} must correspond to those of \code{xstart}.
- As mentioned above, some plug-ins are available for \code{rprocess}.
-
At present, the following methods make use of \code{rprocess}:
\itemize{
\item \code{\link[=simulate-pomp]{simulate}}
@@ -153,8 +157,6 @@
\code{d[j,k]} is the probability density of the transition from state \code{x[,j,k-1]} at time \code{times[k-1]} to state \code{x[,j,k]} at time \code{times[k]}.
If \code{log=TRUE}, then the log of the pdf is returned.
- As mentioned above, some plug-ins are available for \code{dprocess}.
-
\strong{In writing this function, you may assume that the transitions are consecutive.}
It should be quite clear that, but for this assumption, it would be quite difficult in general to write the transition probabilities.
In fact, from one perspective, the algorithms in \pkg{pomp} are designed to overcome just this difficulty.
@@ -223,113 +225,17 @@
}
\value{An object of class \code{pomp}.}
\examples{
-## Simple example: a 2-D Brownian motion, observed with normal error
-## first, set up the pomp object:
-
-rw2 <- pomp(
- rprocess = function (xstart, times, params, ...) {
- ## this function simulates two independent random walks with intensities s1, s2
- nreps <- ncol(params)
- ntimes <- length(times)
- dt <- diff(times)
- x <- array(0,dim=c(2,nreps,ntimes))
- rownames(x) <- rownames(xstart)
- noise.sds <- params[c('s1','s2'),]
- x[,,1] <- xstart
- for (j in 2:ntimes) {
- ## we are mimicking a continuous-time process, so the increments have SD ~ sqrt(dt)
- ## note that we do not have to assume that 'times' are equally spaced
- x[,,j] <- rnorm(n=2*nreps,mean=x[,,j-1],sd=noise.sds*dt[j-1])
- }
- x
- },
- dprocess = function (x, times, params, log, ...) {
- ## given a sequence of consecutive states in 'x', this function computes the p.d.f.
- nreps <- ncol(params)
- ntimes <- length(times)
- dt <- diff(times)
- d <- array(0,dim=c(2,nreps,ntimes-1))
- noise.sds <- params[c('s1','s2'),]
- for (j in 2:ntimes)
- d[,,j-1] <- dnorm(x[,,j]-x[,,j-1],mean=0,sd=noise.sds*dt[j-1],log=TRUE)
- if (log) {
- apply(d,c(2,3),sum)
- } else {
- exp(apply(d,c(2,3),sum))
- }
- },
- measurement.model=list(
- y1 ~ norm(mean=x1,sd=tau),
- y2 ~ norm(mean=x2,sd=tau)
- ),
- skeleton = function (x, t, params, covars, ...) {
- ## since this is a continuous-time process, the skeleton is the vectorfield
- ## since the random walk is unbiased, the vectorfield is identically zero
- rep(0,length=length(x))
- },
- times=1:100,
- data=rbind(
- y1=rep(0,100),
- y2=rep(0,100)
- ),
- t0=0
- )
-
-## specify some parameters
-p <- rbind(s1=c(2,2,3),s2=c(0.1,1,2),tau=c(1,5,0),x1.0=c(0,0,5),x2.0=c(0,0,0))
-
-## simulate
-examples <- simulate(rw2,params=p)
-plot(examples[[1]])
-
-## construct an (almost) equivalent pomp object using a plugin:
-rw2 <- pomp(
- rprocess = euler.simulate,
- step.fun = function (x, t, params, delta.t, ...) {
- y <- rnorm(n=2,mean=x[c("x1","x2")],sd=params[c("s1","s2")]*delta.t)
- names(y) <- c("x1","x2")
- y
- },
- rmeasure = function (x, t, params, ...) {
- y <- rnorm(n=2,mean=x[c("x1","x2")],sd=params["tau"])
- names(y) <- c("y1","y2")
- y
- },
- dmeasure = function (y, x, t, params, log, ...) {
- d <- dnorm(x=y[c("y1","y2")],mean=x[c("x1","x2")],sd=params["tau"],log=TRUE)
- print(d)
- if (log) sum(d,na.rm=T) else exp(sum(d,na.rm=T))
- },
- delta.t=1,
- data=data.frame(
- time=1:100,
- y1=rep(c(1,2,3,4,NA),20),
- y2=0
- ),
- times="time",
- t0=0
- )
-rw2 <- simulate(rw2,params=c(s1=2,s2=0.1,tau=1,x1.0=0,x2.0=0))
-
-## For more examples, see the vignettes.
+## For examples, see the vignettes.
}
\section{Warning}{
Some error checking is done, but complete error checking is impossible.
If the user-specified functions do not conform to the above specifications (see Details), then the results may be invalid.
- Each algorithm that uses a \code{pomp} object uses some subset of the four basic components (\code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}).
- It is unlikely that any algorithm will use all of these.
+ Each algorithm that uses a \code{pomp} object uses some subset of the five basic components (\code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, \code{skeleton}).
}
\author{Aaron A. King \email{kingaa at umich dot edu}}
\seealso{
- \link{pomp-class},
\link{pomp-methods},
- \code{\link{rprocess}},
- \code{\link{dprocess}},
- \code{\link{rmeasure}},
- \code{\link{dmeasure}},
- \code{\link{skeleton}},
- \code{\link{init.state}},
- \link[pomp]{euler}
+ \link[pomp]{plugins}
}
\keyword{models}
\keyword{ts}
Modified: pkg/man/verhulst.Rd
===================================================================
--- pkg/man/verhulst.Rd 2010-05-10 13:50:37 UTC (rev 223)
+++ pkg/man/verhulst.Rd 2010-05-10 13:51:49 UTC (rev 224)
@@ -8,7 +8,7 @@
\usage{data(verhulst)}
\details{
The model is written as an Ito diffusion, \eqn{dn = r n (1-n/K) dt + \sigma n dW}, where \eqn{W} is a Wiener process.
- It is implemented using the \code{\link{euler.simulate}} plug-in.
+ It is implemented using the \code{\link{euler.sim}} plug-in.
}
\examples{
data(verhulst)
Modified: pkg/tests/logistic.R
===================================================================
--- pkg/tests/logistic.R 2010-05-10 13:50:37 UTC (rev 223)
+++ pkg/tests/logistic.R 2010-05-10 13:51:49 UTC (rev 224)
@@ -4,37 +4,39 @@
data=rbind(obs=rep(0,1000)),
times=0.1*seq.int(length=1000),
t0=0,
- rprocess=euler.simulate,
- step.fun=function(x,t,params,delta.t,...){
- with(
- as.list(c(x,params)),
- rnorm(
- n=1,
- mean=n+r*n*(1-n/K)*delta.t,
- sd=sigma*n*sqrt(delta.t)
- )
- )
- },
- dprocess=onestep.density,
- dens.fun=function(x1,x2,t1,t2,params,log,...){
- delta.t <- t2-t1
- with(
- as.list(c(x1,params)),
- dnorm(
- x=x2['n'],
- mean=n+r*n*(1-n/K)*delta.t,
- sd=sigma*n*sqrt(delta.t)
- )
- )
- },
+ rprocess=euler.sim(
+ step.fun=function(x,t,params,delta.t,...){
+ with(
+ as.list(c(x,params)),
+ rnorm(
+ n=1,
+ mean=n+r*n*(1-n/K)*delta.t,
+ sd=sigma*n*sqrt(delta.t)
+ )
+ )
+ },
+ delta.t=0.01
+ ),
+ dprocess=onestep.dens(
+ dens.fun=function(x1,x2,t1,t2,params,log,...){
+ delta.t <- t2-t1
+ with(
+ as.list(c(x1,params)),
+ dnorm(
+ x=x2['n'],
+ mean=n+r*n*(1-n/K)*delta.t,
+ sd=sigma*n*sqrt(delta.t)
+ )
+ )
+ }
+ ),
measurement.model=obs~lnorm(meanlog=log(n),sdlog=log(1+tau)),
skeleton.vectorfield=function(x,t,params,...){
with(
as.list(c(x,params)),
r*n*(1-n/K)
)
- },
- delta.t=0.01
+ }
)
params <- c(n.0=10000,K=10000,r=0.9,sigma=0.4,tau=0.1)
@@ -86,5 +88,6 @@
params <- cbind(c(n.0=100,K=10000,r=0.2,sigma=0.4,tau=0.1),c(n.0=1000,K=11000,r=0.1,sigma=0.4,tau=0.1))
x <- trajectory(po,params=params)
pdf(file='logistic.pdf')
+plot(po)
matplot(time(po,t0=TRUE),t(x['n',,]),type='l')
dev.off()
Modified: pkg/tests/logistic.Rout.save
===================================================================
--- pkg/tests/logistic.Rout.save 2010-05-10 13:50:37 UTC (rev 223)
+++ pkg/tests/logistic.Rout.save 2010-05-10 13:51:49 UTC (rev 224)
@@ -24,37 +24,39 @@
+ data=rbind(obs=rep(0,1000)),
+ times=0.1*seq.int(length=1000),
+ t0=0,
-+ rprocess=euler.simulate,
-+ step.fun=function(x,t,params,delta.t,...){
-+ with(
-+ as.list(c(x,params)),
-+ rnorm(
-+ n=1,
-+ mean=n+r*n*(1-n/K)*delta.t,
-+ sd=sigma*n*sqrt(delta.t)
-+ )
-+ )
-+ },
-+ dprocess=onestep.density,
-+ dens.fun=function(x1,x2,t1,t2,params,log,...){
-+ delta.t <- t2-t1
-+ with(
-+ as.list(c(x1,params)),
-+ dnorm(
-+ x=x2['n'],
-+ mean=n+r*n*(1-n/K)*delta.t,
-+ sd=sigma*n*sqrt(delta.t)
-+ )
-+ )
-+ },
++ rprocess=euler.sim(
++ step.fun=function(x,t,params,delta.t,...){
++ with(
++ as.list(c(x,params)),
++ rnorm(
++ n=1,
++ mean=n+r*n*(1-n/K)*delta.t,
++ sd=sigma*n*sqrt(delta.t)
++ )
++ )
++ },
++ delta.t=0.01
++ ),
++ dprocess=onestep.dens(
++ dens.fun=function(x1,x2,t1,t2,params,log,...){
++ delta.t <- t2-t1
++ with(
++ as.list(c(x1,params)),
++ dnorm(
++ x=x2['n'],
++ mean=n+r*n*(1-n/K)*delta.t,
++ sd=sigma*n*sqrt(delta.t)
++ )
++ )
++ }
++ ),
+ measurement.model=obs~lnorm(meanlog=log(n),sdlog=log(1+tau)),
+ skeleton.vectorfield=function(x,t,params,...){
+ with(
+ as.list(c(x,params)),
+ r*n*(1-n/K)
+ )
-+ },
-+ delta.t=0.01
++ }
+ )
>
> params <- c(n.0=10000,K=10000,r=0.9,sigma=0.4,tau=0.1)
@@ -112,6 +114,7 @@
> params <- cbind(c(n.0=100,K=10000,r=0.2,sigma=0.4,tau=0.1),c(n.0=1000,K=11000,r=0.1,sigma=0.4,tau=0.1))
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 224
More information about the pomp-commits
mailing list