[Pomp-commits] r112 - in pkg: . R data inst inst/doc inst/examples man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Apr 28 14:02:32 CEST 2009
Author: kingaa
Date: 2009-04-28 14:02:31 +0200 (Tue, 28 Apr 2009)
New Revision: 112
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/euler.R
pkg/data/euler.sir.rda
pkg/data/ou2.rda
pkg/inst/ChangeLog
pkg/inst/doc/compiled_code_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/examples/euler_sir.R
pkg/inst/examples/logistic.R
pkg/inst/examples/rw2.R
pkg/man/euler.Rd
pkg/src/euler.c
pkg/tests/logistic.R
pkg/tests/logistic.Rout.save
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
add onestep plugin to facilitate the construction of models for which explicit formulae are available for one-step transitions.
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2009-04-28 11:56:18 UTC (rev 111)
+++ pkg/DESCRIPTION 2009-04-28 12:02:31 UTC (rev 112)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.23-3
-Date: 2009-04-26
+Version: 0.23-4
+Date: 2009-04-28
Author: Aaron A. King, Edward L. Ionides, Carles Martinez Breto, Steve Ellner, Bruce Kendall
Maintainer: Aaron A. King <kingaa at umich.edu>
Description: Inference methods for partially-observed Markov processes
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2009-04-28 11:56:18 UTC (rev 111)
+++ pkg/NAMESPACE 2009-04-28 12:02:31 UTC (rev 112)
@@ -35,7 +35,8 @@
reulermultinom,
deulermultinom,
euler.simulate,
- euler.density,
+ onestep.simulate,
+ onestep.density,
sobol,
bspline.basis,
periodic.bspline.basis,
Modified: pkg/R/euler.R
===================================================================
--- pkg/R/euler.R 2009-04-28 11:56:18 UTC (rev 111)
+++ pkg/R/euler.R 2009-04-28 12:02:31 UTC (rev 112)
@@ -1,3 +1,45 @@
+onestep.simulate <- function (xstart, times, params,
+ step.fun, ...,
+ statenames = character(0),
+ paramnames = character(0),
+ covarnames = character(0),
+ zeronames = character(0),
+ tcovar, covar, PACKAGE)
+{
+ if (is.character(step.fun)) {
+ efun <- try(
+ getNativeSymbolInfo(step.fun,PACKAGE)$address,
+ silent=FALSE
+ )
+ if (inherits(efun,'try-error')) {
+ stop("no symbol named ",step.fun," in package ",PACKAGE)
+ }
+ } else if (is.function(step.fun)) {
+ if (!all(c('x','t','params','delta.t','...')%in%names(formals(step.fun))))
+ stop(sQuote("step.fun")," must be a function of prototype ",sQuote("step.fun(x,t,params,delta.t,...)"))
+ efun <- step.fun
+ } else {
+ stop(sQuote("step.fun")," must be either a function or the name of a compiled routine")
+ }
+
+ .Call(
+ euler_model_simulator,
+ func=efun,
+ xstart=xstart,
+ times=times,
+ params=params,
+ dt=0,
+ method=1L,
+ statenames=statenames,
+ paramnames=paramnames,
+ covarnames=covarnames,
+ zeronames=zeronames,
+ tcovar=tcovar,
+ covar=covar,
+ args=pairlist(...)
+ )
+}
+
euler.simulate <- function (xstart, times, params,
step.fun, delta.t, ...,
statenames = character(0),
@@ -24,28 +66,29 @@
.Call(
euler_model_simulator,
- efun,
- xstart,
- times,
- params,
- delta.t,
- statenames,
- paramnames,
- covarnames,
- zeronames,
- tcovar,
- covar,
+ 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.density <- function (x, times, params,
- dens.fun, ...,
- statenames = character(0),
- paramnames = character(0),
- covarnames = character(0),
- tcovar, covar, log = FALSE,
- PACKAGE)
+onestep.density <- function (x, times, params,
+ dens.fun, ...,
+ statenames = character(0),
+ paramnames = character(0),
+ covarnames = character(0),
+ tcovar, covar, log = FALSE,
+ PACKAGE)
{
if (is.character(dens.fun)) {
efun <- try(
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2009-04-28 11:56:18 UTC (rev 111)
+++ pkg/inst/ChangeLog 2009-04-28 12:02:31 UTC (rev 112)
@@ -1,3 +1,15 @@
+2009-04-28 kingaa
+
+ * [r110] src/dprocess.c, src/rprocess.c: explicitly cast 'times' to
+ double to avoid potential problems when R casts a sequence to
+ integer storage type
+ * [r109] man/pomp-package.Rd:
+
+2009-04-27 kingaa
+
+ * [r108] inst/ChangeLog, inst/doc/Makefile,
+ inst/doc/compiled_code_in_pomp.pdf, inst/doc/intro_to_pomp.pdf:
+
2009-04-26 kingaa
* [r107] DESCRIPTION, inst/ChangeLog: version 0.23-3
Modified: pkg/inst/doc/compiled_code_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/examples/euler_sir.R
===================================================================
--- pkg/inst/examples/euler_sir.R 2009-04-28 11:56:18 UTC (rev 111)
+++ pkg/inst/examples/euler_sir.R 2009-04-28 12:02:31 UTC (rev 112)
@@ -103,7 +103,7 @@
)
},
rprocess=euler.simulate,
- dprocess=euler.density,
+ dprocess=onestep.density,
measurement.model=measles~binom(size=cases,prob=exp(rho)),
initializer=function(params,t0,...){
p <- exp(params)
@@ -150,7 +150,7 @@
step.fun="sir_euler_simulator",
rprocess=euler.simulate,
dens.fun="sir_euler_density",
- dprocess=euler.density,
+ dprocess=onestep.density,
skeleton.vectorfield="sir_ODE",
rmeasure="binom_rmeasure",
dmeasure="binom_dmeasure",
Modified: pkg/inst/examples/logistic.R
===================================================================
--- pkg/inst/examples/logistic.R 2009-04-28 11:56:18 UTC (rev 111)
+++ pkg/inst/examples/logistic.R 2009-04-28 12:02:31 UTC (rev 112)
@@ -15,7 +15,7 @@
)
)
},
- dprocess=euler.density,
+ dprocess=onestep.density,
dens.fun=function(x1,x2,t1,t2,params,log,...){
delta.t <- t2-t1
with(
Modified: pkg/inst/examples/rw2.R
===================================================================
--- pkg/inst/examples/rw2.R 2009-04-28 11:56:18 UTC (rev 111)
+++ pkg/inst/examples/rw2.R 2009-04-28 12:02:31 UTC (rev 112)
@@ -1,5 +1,41 @@
require(pomp)
+## using the "onestep" plugins
+
+rw2 <- pomp(
+ rprocess = onestep.simulate,
+ dprocess = onestep.density,
+ step.fun = function(x, t, params, delta.t, ...) {
+ c(
+ x1=rnorm(n=1,mean=x['x1'],sd=params['s1']*delta.t),
+ x2=rnorm(n=1,mean=x['x2'],sd=params['s2']*delta.t)
+ )
+ },
+ dens.fun = function (x1, t1, x2, t2, params, ...) {
+ sum(
+ dnorm(
+ x=x2[c('x1','x2')],
+ mean=x1[c('x1','x2')],
+ sd=params[c('s1','s2')]*(t2-t1),
+ log=TRUE
+ ),
+ na.rm=TRUE
+ )
+ },
+ measurement.model=list(
+ y1 ~ norm(mean=x1,sd=tau),
+ y2 ~ norm(mean=x2,sd=tau)
+ ),
+ times=1:100,
+ data=rbind(
+ y1=rep(0,100),
+ y2=rep(0,100)
+ ),
+ t0=0
+ )
+
+## writing rprocess and dprocess from scratch
+
rw.rprocess <- function (params, xstart, times, ...) {
## this function simulates two independent random walks with intensities s1, s2
nvars <- nrow(xstart)
@@ -72,34 +108,3 @@
useless=23
)
-po <- pomp(
- rprocess = euler.simulate,
- dprocess = euler.density,
- delta.t = 1,
- step.fun = function(x, t, params, dt, ...) {
- c(
- y1=rnorm(n=1,mean=x['x1'],sd=params['s1']),
- y2=rnorm(n=1,mean=x['x2'],sd=params['s2'])
- )
- },
- dens.fun = function (x1, t1, x2, t2, params, ...) {
- sum(
- dnorm(
- x=x2[c('x1','x2')],
- mean=x1[c('x1','x2')],
- sd=params[c('s1','s2')]
- ),
- na.rm=TRUE
- )
- },
- measurement.model=list(
- y1 ~ norm(mean=x1,sd=tau),
- y2 ~ norm(mean=x2,sd=tau)
- ),
- times=1:100,
- data=rbind(
- y1=rep(0,100),
- y2=rep(0,100)
- ),
- t0=0
- )
Modified: pkg/man/euler.Rd
===================================================================
--- pkg/man/euler.Rd 2009-04-28 11:56:18 UTC (rev 111)
+++ pkg/man/euler.Rd 2009-04-28 12:02:31 UTC (rev 112)
@@ -1,7 +1,8 @@
\name{euler}
\alias{euler}
\alias{euler.simulate}
-\alias{euler.density}
+\alias{onestep.simulate}
+\alias{onestep.density}
\title{Dynamical models based on stochastic Euler algorithms}
\description{
Facilities for simulating discrete-time Markov processes and continuous-time Markov processes using the Euler algorithm.
@@ -11,10 +12,14 @@
statenames = character(0), paramnames = character(0),
covarnames = character(0), zeronames = character(0),
tcovar, covar, PACKAGE)
-euler.density(x, times, params, dens.fun, \dots,
- statenames = character(0), paramnames = character(0),
- covarnames = character(0), tcovar, covar, log = FALSE,
- PACKAGE)
+onestep.simulate(xstart, times, params, step.fun, \dots,
+ statenames = character(0), paramnames = character(0),
+ covarnames = character(0), zeronames = character(0),
+ tcovar, covar, PACKAGE)
+onestep.density(x, times, params, dens.fun, \dots,
+ statenames = character(0), paramnames = character(0),
+ covarnames = character(0), tcovar, covar, log = FALSE,
+ PACKAGE)
}
\arguments{
\item{xstart}{
@@ -37,8 +42,8 @@
For details on how to write such codes, see Details.
}
\item{dens.fun}{
- This can be either an R function or a compiled, dynamically loaded native function containing the model transition probability density function.
- This function will be called to compute the probability of the actual Euler steps.
+ This can be either an R function or a compiled, dynamically loaded native function containing the model transition log probability density function.
+ This function will be called to compute the log likelihood of the actual Euler steps.
It must be of type "euler\_step\_pdf" as defined in the header "pomp.h", which is included with the package.
For details on how to write such codes, see Details.
}
@@ -46,13 +51,13 @@
Time interval of Euler steps.
}
\item{statenames}{
- Names of state variables, in the order they will be expected by the routine named in \code{euler.step.fun} and \code{euler.dens.fun}.
+ Names of state variables, in the order they will be expected by the routine named in \code{step.fun} and \code{dens.fun}.
}
\item{paramnames}{
- Names of parameters, in the order they will be expected by the routine named in \code{euler.step.fun} and \code{euler.dens.fun}.
+ Names of parameters, in the order they will be expected by the routine named in \code{step.fun} and \code{dens.fun}.
}
\item{covarnames}{
- Names of columns of the matrix of covariates \code{covar}, in the order they will be expected by the routine named in \code{euler.step.fun} and \code{euler.dens.fun}.
+ Names of columns of the matrix of covariates \code{covar}, in the order they will be expected by the routine named in \code{step.fun} and \code{dens.fun}.
}
\item{zeronames}{
Names of additional variables which will be zeroed before each time in \code{times}.
@@ -78,6 +83,11 @@
}
}
\details{
+ \code{onestep.simulate} assumes that a single call to \code{step.fun} will advance the state process from one time to the next.
+ \code{euler.simulate} will take multiple Euler steps, each of size at most \code{delta.t} (see below for information on how the actual Euler step size is chosen) to get from one time to the next.
+
+ \code{onestep.density} assumes that no state transitions occure between consecutive times.
+
If \code{step.fun} is written as an R function, it must have at least the arguments \code{x}, \code{t}, \code{params}, \code{delta.t}, and \code{\dots}.
On a call to this function, \code{x} will be a named vector of state variables, \code{t} a scalar time, and \code{params} a named vector of parameters.
The length of the Euler step will be \code{delta.t}.
@@ -97,10 +107,10 @@
If \code{dens.fun} is written in a native language, it must be a function of type "euler\_step\_pdf" as defined in the header "pomp.h" included with the package (see the directory "include" in the installed package directory).
}
\value{
- \code{euler.simulate} returns a \code{nvar} x \code{nrep} x \code{ntimes} array, where \code{nvar} is the number of state variables, \code{nrep} is the number of replicate simulations (= number of columns of \code{xstart} and \code{params}), and \code{ntimes} is the length of \code{times}.
+ \code{euler.simulate} and \code{onestep.simulate} each return a \code{nvar} x \code{nrep} x \code{ntimes} array, where \code{nvar} is the number of state variables, \code{nrep} is the number of replicate simulations (= number of columns of \code{xstart} and \code{params}), and \code{ntimes} is the length of \code{times}.
If \code{x} is this array, \code{x[,,1]} will be identical to \code{xstart}; the rownames of \code{x} and \code{xstart} will also coincide.
- \code{euler.density} returns a \code{nrep} x \code{ntimes-1} array.
+ \code{onestep.density} returns a \code{nrep} x \code{ntimes-1} array.
If \code{f} is this array, \code{f[i,j]} is the likelihood of a transition from \code{x[,i,j]} to \code{x[,i,j+1]} in exactly one Euler step of duration \code{times[j+1]-times[j]}.
}
\author{Aaron A. King (kingaa at umich dot edu)}
Modified: pkg/src/euler.c
===================================================================
--- pkg/src/euler.c 2009-04-28 11:56:18 UTC (rev 111)
+++ pkg/src/euler.c 2009-04-28 12:02:31 UTC (rev 112)
@@ -92,6 +92,62 @@
}
}
+// take one step from t1 to t2
+static void onestep_simulator (euler_step_sim *estep,
+ double *x, double *xstart, double *times, double *params,
+ int *ndim,
+ int *stateindex, int *parindex, int *covindex, int *zeroindex,
+ double *time_table, double *covar_table)
+{
+ double t, *xp, *pp;
+ int nvar = ndim[0];
+ int npar = ndim[1];
+ int nrep = ndim[2];
+ int ntimes = ndim[3];
+ int covlen = ndim[4];
+ int covdim = ndim[5];
+ int nzero = ndim[6];
+ double covar_fn[covdim];
+ int j, k, p, step, neuler;
+ double dt, tol;
+
+ struct lookup_table covariate_table = {covlen, covdim, 0, time_table, covar_table};
+
+ // copy the start values into the result array
+ for (p = 0; p < nrep; p++)
+ for (k = 0; k < nvar; k++)
+ x[k+nvar*p] = xstart[k+nvar*p];
+
+ // loop over times
+ for (step = 1; step < ntimes; step++) {
+
+ R_CheckUserInterrupt();
+
+ t = times[step-1];
+ dt = times[step]-t;
+
+ // interpolate the covar functions for the covariates
+ if (covdim > 0)
+ table_lookup(&covariate_table,t,covar_fn,0);
+
+ for (p = 0; p < nrep; p++) {
+ xp = &x[nvar*(p+nrep*step)];
+ // copy in the previous values of the state variables
+ for (k = 0; k < nvar; k++)
+ xp[k] = x[k+nvar*(p+nrep*(step-1))];
+ // set some variables to zero
+ for (k = 0; k < nzero; k++)
+ xp[zeroindex[k]] = 0.0;
+
+ pp = ¶ms[npar*p];
+ xp = &x[nvar*(p+nrep*step)];
+
+ (*estep)(xp,pp,stateindex,parindex,covindex,covdim,covar_fn,t,dt);
+
+ }
+ }
+}
+
// these global objects will pass the needed information to the user-defined function (see 'default_euler_step_fn')
// each of these is allocated once, globally, and refilled many times
static SEXP euler_step_Xvec; // state variable vector
@@ -174,7 +230,7 @@
SEXP euler_model_simulator (SEXP func,
SEXP xstart, SEXP times, SEXP params,
- SEXP dt,
+ SEXP dt, SEXP method,
SEXP statenames, SEXP paramnames, SEXP covarnames, SEXP zeronames,
SEXP tcovar, SEXP covar, SEXP args)
{
@@ -191,12 +247,15 @@
SEXP X, pindex, sindex, cindex, zindex;
int *sidx, *pidx, *cidx, *zidx;
SEXP fn, Pnames, Cnames;
+ int do_euler = 1, *meth;
dim = INTEGER(GET_DIM(xstart)); nvar = dim[0]; nrep = dim[1];
dim = INTEGER(GET_DIM(params)); npar = dim[0];
dim = INTEGER(GET_DIM(covar)); covlen = dim[0]; covdim = dim[1];
ntimes = LENGTH(times);
+ if (*(INTEGER(AS_INTEGER(method)))) do_euler = 0;
+
PROTECT(VNAMES = GET_ROWNAMES(GET_DIMNAMES(xstart))); nprotect++;
PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(covar))); nprotect++;
@@ -277,10 +336,16 @@
if (use_native) GetRNGstate();
- euler_simulator(ff,REAL(X),REAL(xstart),REAL(times),REAL(params),
- ndim,REAL(dt),sidx,pidx,cidx,zidx,
- REAL(tcovar),REAL(covar));
-
+ if (do_euler) {
+ euler_simulator(ff,REAL(X),REAL(xstart),REAL(times),REAL(params),
+ ndim,REAL(dt),sidx,pidx,cidx,zidx,
+ REAL(tcovar),REAL(covar));
+ } else {
+ onestep_simulator(ff,REAL(X),REAL(xstart),REAL(times),REAL(params),
+ ndim,sidx,pidx,cidx,zidx,
+ REAL(tcovar),REAL(covar));
+ }
+
if (use_native) PutRNGstate();
if (VINDEX != 0) Free(VINDEX);
Modified: pkg/tests/logistic.R
===================================================================
--- pkg/tests/logistic.R 2009-04-28 11:56:18 UTC (rev 111)
+++ pkg/tests/logistic.R 2009-04-28 12:02:31 UTC (rev 112)
@@ -15,7 +15,7 @@
)
)
},
- dprocess=euler.density,
+ dprocess=onestep.density,
dens.fun=function(x1,x2,t1,t2,params,log,...){
delta.t <- t2-t1
with(
Modified: pkg/tests/logistic.Rout.save
===================================================================
--- pkg/tests/logistic.Rout.save 2009-04-28 11:56:18 UTC (rev 111)
+++ pkg/tests/logistic.Rout.save 2009-04-28 12:02:31 UTC (rev 112)
@@ -1,5 +1,5 @@
-R version 2.7.1 (2008-06-23)
+R version 2.8.1 (2008-12-22)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
@@ -16,7 +16,8 @@
Type 'q()' to quit R.
> library(pomp)
-Loading required package: odesolve
+Loading required package: deSolve
+Loading required package: subplex
>
> po <- pomp(
+ data=rbind(obs=rep(0,1000)),
@@ -33,7 +34,7 @@
+ )
+ )
+ },
-+ dprocess=euler.density,
++ dprocess=onestep.density,
+ dens.fun=function(x1,x2,t1,t2,params,log,...){
+ delta.t <- t2-t1
+ with(
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2009-04-28 11:56:18 UTC (rev 111)
+++ pkg/tests/sir.R 2009-04-28 12:02:31 UTC (rev 112)
@@ -104,7 +104,7 @@
)
},
rprocess=euler.simulate,
- dprocess=euler.density,
+ dprocess=onestep.density,
measurement.model=measles~binom(size=cases,prob=exp(rho)),
initializer=function(params,t0,...){
p <- exp(params)
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2009-04-28 11:56:18 UTC (rev 111)
+++ pkg/tests/sir.Rout.save 2009-04-28 12:02:31 UTC (rev 112)
@@ -123,7 +123,7 @@
+ )
+ },
+ rprocess=euler.simulate,
-+ dprocess=euler.density,
++ dprocess=onestep.density,
+ measurement.model=measles~binom(size=cases,prob=exp(rho)),
+ initializer=function(params,t0,...){
+ p <- exp(params)
@@ -148,7 +148,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 4.005935 secs
+Time difference of 2.599395 secs
>
> pdf(file='sir.pdf')
>
@@ -165,7 +165,7 @@
> X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 9.746132 secs
+Time difference of 5.949128 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
@@ -223,7 +223,7 @@
> x <- simulate(euler.sir,nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.2860382 secs
+Time difference of 0.1701131 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
@@ -231,7 +231,7 @@
> X4 <- trajectory(euler.sir,times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 7.802866 secs
+Time difference of 4.644418 secs
> plot(t3,X4['I',1,],type='l')
>
> f2 <- dprocess(
More information about the pomp-commits
mailing list