[Pomp-commits] r219 - in pkg: . R data inst inst/doc man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 6 14:12:29 CEST 2010
Author: kingaa
Date: 2010-05-06 14:12:28 +0200 (Thu, 06 May 2010)
New Revision: 219
Added:
pkg/R/pomp-fun.R
pkg/man/pomp-fun.Rd
Modified:
pkg/DESCRIPTION
pkg/R/aaa.R
pkg/R/pfilter.R
pkg/R/pomp-methods.R
pkg/R/pomp.R
pkg/data/euler.sir.rda
pkg/data/ou2.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/ChangeLog
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/man/pfilter.Rd
pkg/man/pomp-class.Rd
pkg/tests/rw2.Rout.save
pkg/tests/sir.Rout.save
Log:
- add a new slot, 'call', to the pomp class
- improve the way that pomps are printed and 'show'n
- streamline pomp construction using pomp.fun
- save Np and tol in the list returned by 'pfilter'
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-04-29 22:49:23 UTC (rev 218)
+++ pkg/DESCRIPTION 2010-05-06 12:12:28 UTC (rev 219)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.28-4
-Date: 2010-04-29
+Version: 0.28-5
+Date: 2010-05-06
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall
Maintainer: Aaron A. King <kingaa at umich.edu>
Description: Inference methods for partially-observed Markov processes
Modified: pkg/R/aaa.R
===================================================================
--- pkg/R/aaa.R 2010-04-29 22:49:23 UTC (rev 218)
+++ pkg/R/aaa.R 2010-05-06 12:12:28 UTC (rev 219)
@@ -34,7 +34,8 @@
paramnames = 'character',
covarnames = 'character',
PACKAGE = 'character',
- userdata = 'list'
+ userdata = 'list',
+ call = "call"
)
)
Modified: pkg/R/pfilter.R
===================================================================
--- pkg/R/pfilter.R 2010-04-29 22:49:23 UTC (rev 218)
+++ pkg/R/pfilter.R 2010-05-06 12:12:28 UTC (rev 219)
@@ -261,6 +261,8 @@
cond.loglik=loglik,
states=xparticles,
seed=seed,
+ Np=Np,
+ tol=tol,
nfail=nfail,
loglik=sum(loglik)
)
Added: pkg/R/pomp-fun.R
===================================================================
--- pkg/R/pomp-fun.R (rev 0)
+++ pkg/R/pomp-fun.R 2010-05-06 12:12:28 UTC (rev 219)
@@ -0,0 +1,58 @@
+pomp.fun <- function (f = NULL, PACKAGE, proto = NULL) {
+ if (missing(PACKAGE))
+ PACKAGE <- character(0)
+ if (is.function(f)) {
+ if (!is.null(proto)) {
+ prototype <- strsplit(as.character(proto),split="\\(|\\)|\\,")
+ fname <- prototype[1]
+ args <- prototype[-1]
+ if (!all(args%in%names(formals(f))))
+ stop(sQuote(fname)," must be a function of prototype ",deparse(proto),call.=FALSE)
+ }
+ retval <- new(
+ "pomp.fun",
+ R.fun=f,
+ native.fun=character(0),
+ PACKAGE=PACKAGE,
+ use=1L
+ )
+ } else if (is.character(f)) {
+ retval <- new(
+ "pomp.fun",
+ R.fun=function(...)stop("unreachable error: please report this bug!"),
+ native.fun=f,
+ PACKAGE=PACKAGE,
+ use=2L
+ )
+ } else {
+ retval <- new(
+ "pomp.fun",
+ R.fun=function(...)stop(sQuote(fname)," not specified"),
+ native.fun=character(0),
+ PACKAGE=PACKAGE,
+ use=1L
+ )
+ }
+ retval
+}
+
+setMethod(
+ 'show',
+ 'pomp.fun',
+ function (object) {
+ if (object at use==1) {
+ show(object at R.fun)
+ } else {
+ cat("native function ",sQuote(object at native.fun),sep="")
+ if (length(object at PACKAGE)>0)
+ cat(", dynamically loaded from ",sQuote(object at PACKAGE),sep="")
+ cat ("\n")
+ }
+ }
+ )
+
+setMethod(
+ 'print',
+ 'pomp.fun',
+ function (x, ...) show(x)
+ )
Modified: pkg/R/pomp-methods.R
===================================================================
--- pkg/R/pomp-methods.R 2010-04-29 22:49:23 UTC (rev 218)
+++ pkg/R/pomp-methods.R 2010-05-06 12:12:28 UTC (rev 219)
@@ -157,7 +157,10 @@
'print',
'pomp',
function (x, ...) {
+ cat("data and states:\n")
print(as(x,'data.frame'))
+ cat("\ncall:\n")
+ print(x at call)
invisible(x)
}
)
@@ -175,31 +178,23 @@
cat ("parameter(s) unspecified\n");
}
cat("process model simulator, rprocess = \n")
- print(object at rprocess)
+ show(object at rprocess)
cat("process model density, dprocess = \n")
- print(object at dprocess)
+ show(object at dprocess)
cat("measurement model simulator, rmeasure = \n")
- if (object at rmeasure@use==2) {
- cat("native function, '",object at rmeasure@native.fun,"'",sep="")
- if (length(object at rmeasure@PACKAGE)>0)
- cat(", PACKAGE = '",object at rmeasure@PACKAGE,"'",sep="")
- cat ("\n")
- } else {
- print(object at rmeasure@R.fun)
- }
+ show(object at rmeasure)
cat("measurement model density, dmeasure = \n")
- if (object at dmeasure@use==2) {
- cat("native function, '",object at dmeasure@native.fun,"'",sep="")
- if (length(object at dmeasure@PACKAGE)>0)
- cat(", PACKAGE = '",object at dmeasure@PACKAGE,"'",sep="")
- cat ("\n")
- } else {
- print(object at dmeasure@R.fun)
+ show(object at dmeasure)
+ if (!is.na(object at skeleton.type)) {
+ cat("skeleton (",object at skeleton.type,") = \n")
+ show(object at skeleton)
}
cat("initializer = \n")
- print(object at initializer)
- cat("userdata = \n")
- show(object at userdata)
+ show(object at initializer)
+ if (length(object at userdata)>0) {
+ cat("userdata = \n")
+ show(object at userdata)
+ }
invisible(NULL)
}
)
Modified: pkg/R/pomp.R
===================================================================
--- pkg/R/pomp.R 2010-04-29 22:49:23 UTC (rev 218)
+++ pkg/R/pomp.R 2010-05-06 12:12:28 UTC (rev 219)
@@ -13,6 +13,9 @@
skeleton.map, skeleton.vectorfield, initializer, covar, tcovar,
statenames, paramnames, covarnames,
PACKAGE) {
+ ## save the call
+ this.call <- match.call()
+
## check the data
if (is.data.frame(data)) {
if (!is.character(times) || length(times)!=1 || !(times%in%names(data)))
@@ -48,7 +51,10 @@
if (!missing(measurement.model)) {
if (!(missing(dmeasure)&&missing(rmeasure))) {
- warning("specifying ",sQuote("measurement.model")," overrides specification of ",sQuote("rmeasure")," and ",sQuote("dmeasure"))
+ warning(
+ "specifying ",sQuote("measurement.model"),
+ " overrides specification of ",sQuote("rmeasure")," and ",sQuote("dmeasure")
+ )
}
mm <- measform2pomp(measurement.model)
rmeasure <- mm$rmeasure
@@ -58,65 +64,20 @@
if (missing(rmeasure))
rmeasure <- function(x,t,params,covars,...)stop(sQuote("rmeasure")," not specified")
if (missing(dmeasure))
- dmeasure <- function(y,x,t,params,log=FALSE,covars,...)stop(sQuote("dmeasure")," not specified")
+ dmeasure <- function(y,x,t,params,log,covars,...)stop(sQuote("dmeasure")," not specified")
if (missing(skeleton.map)) {
if (missing(skeleton.vectorfield)) {# skeleton is unspecified
skeleton.type <- as.character(NA)
- skeleton <- new(
- "pomp.fun",
- R.fun=function(x,t,params,covars,...)stop(sQuote("skeleton")," not specified"),
- use=as.integer(1)
- )
+ skeleton <- pomp.fun(f=function(x,t,params,covars,...)stop(sQuote("skeleton")," not specified"))
} else { # skeleton is a vectorfield (ordinary differential equation)
skeleton.type <- "vectorfield"
- if (is.function(skeleton.vectorfield)) {
- if (!all(c('x','t','params','...')%in%names(formals(skeleton.vectorfield))))
- stop(
- sQuote("skeleton.vectorfield"),
- " must be a function of prototype ",
- sQuote("skeleton.vectorfield(x,t,params,...)"),
- call.=TRUE
- )
- skeleton <- new(
- "pomp.fun",
- R.fun=skeleton.vectorfield,
- use=as.integer(1)
- )
- } else if (is.character(skeleton.vectorfield)) {
- skeleton <- new(
- "pomp.fun",
- R.fun=function(x,t,params,...)stop("very bad: skel.fun"),
- native.fun=skeleton.vectorfield,
- PACKAGE=PACKAGE,
- use=as.integer(2)
- )
- } else {
- stop("pomp error: ",sQuote("skeleton.vectorfield")," must be either a function or the name of a compiled routine")
- }
+ skeleton <- pomp.fun(f=skeleton.vectorfield,PACKAGE=PACKAGE,proto="skeleton.vectorfield(x,t,params,...)")
}
} else {
if (missing(skeleton.vectorfield)) { # skeleton is a map (discrete-time system)
skeleton.type <- "map"
- if (is.function(skeleton.map)) {
- if (!all(c('x','t','params','...')%in%names(formals(skeleton.map))))
- stop("pomp error: ",sQuote("skeleton.map")," must be a function of prototype ",sQuote("skeleton.map(x,t,params,...)"))
- skeleton <- new(
- "pomp.fun",
- R.fun=skeleton.map,
- use=as.integer(1)
- )
- } else if (is.character(skeleton.map)) {
- skeleton <- new(
- "pomp.fun",
- R.fun=function(x,t,params,...)stop("very bad: skel.fun"),
- native.fun=skeleton.map,
- PACKAGE=PACKAGE,
- use=as.integer(2)
- )
- } else {
- stop("pomp error: ",sQuote("skeleton.map")," must be either a function or the name of a compiled routine")
- }
+ skeleton <- pomp.fun(f=skeleton.map,PACKAGE=PACKAGE,proto="skeleton.map(x,t,params,...)")
} else { # a dynamical system cannot be both a map and a vectorfield
stop("pomp error: it is not permitted to specify both ",sQuote("skeleton.map")," and ",sQuote("skeleton.vectorfield"))
}
@@ -139,46 +100,9 @@
call.=TRUE
)
- if (is.function(rmeasure)) {
- if (!all(c('x','t','params','...')%in%names(formals(rmeasure))))
- stop("pomp error: ",sQuote("rmeasure")," must be a function of prototype ",sQuote("rmeasure(x,t,params,...)"))
- rmeasure <- new(
- "pomp.fun",
- R.fun=rmeasure,
- use=as.integer(1)
- )
- } else if (is.character(rmeasure)) {
- rmeasure <- new(
- "pomp.fun",
- R.fun=function(x,t,params,...)stop("very bad: rmeasure.fun"),
- native.fun=rmeasure,
- PACKAGE=PACKAGE,
- use=as.integer(2)
- )
- } else {
- stop("pomp error: ",sQuote("rmeasure")," must be either a function or the name of a compiled routine")
- }
+ rmeasure <- pomp.fun(f=rmeasure,PACKAGE=PACKAGE,proto="rmeasure(x,t,params,...)")
+ dmeasure <- pomp.fun(f=dmeasure,PACKAGE=PACKAGE,proto="dmeasure(y,x,t,params,log,...)")
- if (is.function(dmeasure)) {
- if (!all(c('y','x','t','params','log','...')%in%names(formals(dmeasure))))
- stop("pomp error: ",sQuote("dmeasure")," must be a function of prototype ",sQuote("dmeasure(y,x,t,params,log,...)"))
- dmeasure <- new(
- "pomp.fun",
- R.fun=dmeasure,
- use=as.integer(1)
- )
- } else if (is.character(dmeasure)) {
- dmeasure <- new(
- "pomp.fun",
- R.fun=function(y,x,t,params,log,...)stop("very bad: rmeasure.fun"),
- native.fun=dmeasure,
- PACKAGE=PACKAGE,
- use=as.integer(2)
- )
- } else {
- stop("pomp error: ",sQuote("dmeasure")," must be either a function or the name of a compiled routine")
- }
-
if (!is.function(initializer))
stop(
"pomp error: ",sQuote("initializer")," must be a function",
@@ -280,7 +204,8 @@
paramnames = paramnames,
covarnames = covarnames,
PACKAGE = PACKAGE,
- userdata = list(...)
+ userdata = list(...),
+ call=this.call
)
}
Modified: pkg/data/euler.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/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2010-04-29 22:49:23 UTC (rev 218)
+++ pkg/inst/ChangeLog 2010-05-06 12:12:28 UTC (rev 219)
@@ -1,5 +1,8 @@
2010-04-29 kingaa
+ * [r218] tests/pfilter.R, tests/pfilter.Rout.save: - tests for
+ 'pfilter'
+ * [r217] inst/ChangeLog: - version 0.28-4
* [r216] DESCRIPTION, R/pfilter.R, R/trajmatch.R,
data/euler.sir.rda, data/ou2.rda, data/rw2.rda,
data/verhulst.rda, inst/ChangeLog, inst/data-R/ou2.R,
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/man/pfilter.Rd
===================================================================
--- pkg/man/pfilter.Rd 2010-04-29 22:49:23 UTC (rev 218)
+++ pkg/man/pfilter.Rd 2010-05-06 12:12:28 UTC (rev 219)
@@ -94,8 +94,8 @@
If the argument \code{seed} was specified, this is a copy;
if not, this is the internal state of the random number generator at the time of call.
}
- \item{nfail}{
- The number of filtering failures encountered.
+ \item{Np, tol, nfail}{
+ The number of particles used, failure tolerance, and number of filtering failures, respectively.
}
\item{loglik}{
The estimated log-likelihood.
Modified: pkg/man/pomp-class.Rd
===================================================================
--- pkg/man/pomp-class.Rd 2010-04-29 22:49:23 UTC (rev 218)
+++ pkg/man/pomp-class.Rd 2010-05-06 12:12:28 UTC (rev 219)
@@ -63,6 +63,9 @@
A list containing any objects the user desires.
Using this mechanism, the user can store additional information necessary for the definition of the model.
}
+ \item{call}{
+ The call that created the \code{pomp} object.
+ }
}
}
\section{Methods}{
Added: pkg/man/pomp-fun.Rd
===================================================================
--- pkg/man/pomp-fun.Rd (rev 0)
+++ pkg/man/pomp-fun.Rd 2010-05-06 12:12:28 UTC (rev 219)
@@ -0,0 +1,35 @@
+\name{pomp-fun}
+\docType{methods}
+\alias{pomp.fun-methods}
+\alias{show,pomp.fun-method}
+\alias{show-pomp.fun}
+\alias{print,pomp.fun-method}
+\alias{print-pomp.fun}
+\alias{pomp.fun-class}
+\alias{pomp.fun}
+\keyword{internal}
+\title{Definition and methods of the "pomp.fun" class}
+\description{Definition and methods of the "pomp.fun" class}
+\usage{
+pomp.fun(f = NULL, PACKAGE, proto = NULL)
+\S4method{show}{pomp.fun}(object)
+\S4method{print}{pomp.fun}(x, \dots)
+}
+\arguments{
+ \item{f}{A function or the name of a native routine.}
+ \item{PACKAGE}{
+ optional; the name of the dynamically-loadable library in which the native function \code{f} can be found.
+ }
+ \item{proto}{
+ optional string; a prototype against which \code{f} will be checked.
+ }
+ \item{object, x}{The \code{pomp.fun} object.}
+}
+\details{
+ The sQuote{pomp.fun} class helps to settle common issues associated with user-defined functions which can be defined either via R code or by a native, compiled routine.
+ It is not exported to userland.
+}
+\author{Aaron A. King \email{kingaa at umich dot edu}}
+\seealso{
+ \code{\link{pomp}}
+}
Modified: pkg/tests/rw2.Rout.save
===================================================================
--- pkg/tests/rw2.Rout.save 2010-04-29 22:49:23 UTC (rev 218)
+++ pkg/tests/rw2.Rout.save 2010-05-06 12:12:28 UTC (rev 219)
@@ -105,6 +105,7 @@
+ )
>
> show(rw2)
+data and states:
time y1 y2
1 1 0 0
2 2 0 0
@@ -206,6 +207,12 @@
98 98 0 0
99 99 0 0
100 100 0 0
+
+call:
+pomp(data = rbind(y1 = rep(0, 100), y2 = rep(0, 100)), times = 1:100,
+ t0 = 0, useless = 23, rprocess = rw.rprocess, dprocess = rw.dprocess,
+ measurement.model = list(y1 ~ norm(mean = x1, sd = tau),
+ y2 ~ norm(mean = x2, sd = tau)), initializer = bad.initializer)
zero time, t0 = 0
parameter(s) unspecified
process model simulator, rprocess =
@@ -251,7 +258,7 @@
}
y
}
-<environment: 0x189dbb0>
+<environment: 0x160cbd8>
measurement model density, dmeasure =
function (y, x, t, params, log, covars, ...)
{
@@ -264,7 +271,7 @@
f
else exp(f)
}
-<environment: 0x189dbb0>
+<environment: 0x160cbd8>
initializer =
function (params, t0, ...)
{
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2010-04-29 22:49:23 UTC (rev 218)
+++ pkg/tests/sir.Rout.save 2010-05-06 12:12:28 UTC (rev 219)
@@ -144,6 +144,7 @@
+ )
>
> show(po)
+data and states:
time measles
1 0.01923077 0
2 0.03846154 0
@@ -353,6 +354,71 @@
206 3.96153846 0
207 3.98076923 0
208 4.00000000 0
+
+call:
+pomp(data = rbind(measles = numeric(52 * 4)), times = 1/52 *
+ seq.int(length = 4 * 52), t0 = 0, delta.t = 1/52/20, zeronames = c("cases"),
+ step.fun = function(t, x, params, covars, delta.t, ...) {
+ params <- exp(params)
+ with(as.list(c(x, params)), {
+ beta <- exp(sum(log(c(beta1, beta2, beta3)) * covars))
+ beta.var <- beta.sd^2
+ dW <- rgamma(n = 1, shape = delta.t/beta.var, scale = beta.var)
+ foi <- (iota + beta * I * dW/delta.t)/pop
+ trans <- c(rpois(n = 1, lambda = mu * pop * delta.t),
+ reulermultinom(n = 1, size = S, rate = c(foi,
+ mu), dt = delta.t), reulermultinom(n = 1, size = I,
+ rate = c(gamma, mu), dt = delta.t), reulermultinom(n = 1,
+ size = R, rate = c(mu), dt = delta.t))
+ c(S = S + trans[1] - trans[2] - trans[3], I = I +
+ trans[2] - trans[4] - trans[5], R = R + trans[4] -
+ trans[6], cases = cases + trans[4], W = if (beta.sd >
+ 0) W + (dW - delta.t)/beta.sd else W, B = trans[1],
+ SI = trans[2], SD = trans[3], IR = trans[4],
+ ID = trans[5], RD = trans[6], dW = dW)
+ })
+ }, dens.fun = function(t1, t2, params, x1, x2, covars, ...) {
+ params <- exp(params)
+ with(as.list(params), {
+ dt <- t2 - t1
+ beta <- exp(sum(log(c(beta1, beta2, beta3)) * covars))
+ beta.var <- beta.sd^2
+ dW <- x2["dW"]
+ foi <- (iota + beta * x1["I"] * dW/dt)/pop
+ probs <- c(dpois(x = x2["B"], lambda = mu * pop *
+ dt, log = T), deulermultinom(x = x2[c("SI", "SD")],
+ size = x1["S"], rate = c(foi, mu), dt = dt, log = T),
+ deulermultinom(x = x2[c("IR", "ID")], size = x1["I"],
+ rate = c(gamma, mu), dt = dt, log = T), deulermultinom(x = x2["RD"],
+ size = x1["R"], rate = c(mu), dt = dt, log = T),
+ dgamma(x = dW, shape = dt/beta.var, scale = beta.var,
+ log = T))
+ sum(probs)
+ })
+ }, rprocess = euler.simulate, dprocess = onestep.density,
+ measurement.model = measles ~ binom(size = cases, prob = exp(rho)),
+ skeleton.vectorfield = function(x, t, params, covars, ...) {
+ xdot <- rep(0, length(x))
+ params <- exp(params)
+ with(as.list(c(x, params)), {
+ beta <- exp(sum(log(c(beta1, beta2, beta3)) * covars))
+ foi <- (iota + beta * I)/pop
+ terms <- c(mu * pop, foi * S, mu * S, gamma * I,
+ mu * I, mu * R)
+ xdot[1:4] <- c(terms[1] - terms[2] - terms[3], terms[2] -
+ terms[4] - terms[5], terms[4] - terms[6], terms[4])
+ xdot
+ })
+ }, initializer = function(params, t0, ...) {
+ p <- exp(params)
+ with(as.list(p), {
+ fracs <- c(S.0, I.0, R.0)
+ x0 <- c(round(pop * fracs/sum(fracs)), rep(0, 9))
+ names(x0) <- c("S", "I", "R", "cases", "W", "B",
+ "SI", "SD", "IR", "ID", "RD", "dW")
+ x0
+ })
+ }, covar = basis, tcovar = tbasis)
zero time, t0 = 0
parameter(s) unspecified
process model simulator, rprocess =
@@ -423,7 +489,7 @@
}
y
}
-<environment: 0x325f890>
+<environment: 0x116d468>
measurement model density, dmeasure =
function (y, x, t, params, log, covars, ...)
{
@@ -436,7 +502,22 @@
f
else exp(f)
}
-<environment: 0x325f890>
+<environment: 0x116d468>
+skeleton ( vectorfield ) =
+function (x, t, params, covars, ...)
+{
+ xdot <- rep(0, length(x))
+ params <- exp(params)
+ with(as.list(c(x, params)), {
+ beta <- exp(sum(log(c(beta1, beta2, beta3)) * covars))
+ foi <- (iota + beta * I)/pop
+ terms <- c(mu * pop, foi * S, mu * S, gamma * I, mu *
+ I, mu * R)
+ xdot[1:4] <- c(terms[1] - terms[2] - terms[3], terms[2] -
+ terms[4] - terms[5], terms[4] - terms[6], terms[4])
+ xdot
+ })
+}
initializer =
function (params, t0, ...)
{
@@ -507,7 +588,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 2.060826 secs
+Time difference of 4.335154 secs
>
> pdf(file='sir.pdf')
>
@@ -524,7 +605,7 @@
> X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 4.699151 secs
+Time difference of 10.29815 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
@@ -584,7 +665,7 @@
> x <- simulate(po,nsim=100)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 2.14945 secs
+Time difference of 8.903283 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
@@ -592,7 +673,7 @@
> X4 <- trajectory(po,times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.734072 secs
+Time difference of 1.693917 secs
> plot(t3,X4['I',1,],type='l')
>
> g2 <- dmeasure(
More information about the pomp-commits
mailing list