[Pomp-commits] r627 - in pkg: . R data inst inst/data-R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 15 18:17:35 CET 2012
Author: kingaa
Date: 2012-03-15 18:17:34 +0100 (Thu, 15 Mar 2012)
New Revision: 627
Added:
pkg/inst/data-R/sir.R
pkg/tests/dacca.R
pkg/tests/dacca.Rout.save
Removed:
pkg/inst/data-R/euler.sir.R
pkg/inst/data-R/gillespie.sir.R
Modified:
pkg/DESCRIPTION
pkg/R/simulate-pomp.R
pkg/R/trajectory-pomp.R
pkg/data/blowflies.rda
pkg/data/dacca.rda
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/gompertz.rda
pkg/data/ou2.rda
pkg/data/ricker.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/NEWS
pkg/inst/TODO
pkg/man/simulate-pomp.Rd
pkg/man/sir.Rd
pkg/man/trajectory-pomp.Rd
pkg/tests/ricker.R
pkg/tests/ricker.Rout.save
Log:
- add new 'as.data.frame' argument to 'simulate' and 'trajectory'
- consolidate both SIR models in one inst/data-R file
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-03-09 16:00:59 UTC (rev 626)
+++ pkg/DESCRIPTION 2012-03-15 17:17:34 UTC (rev 627)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.40-8
-Date: 2012-03-09
+Version: 0.40-9
+Date: 2012-03-15
Revision: $Rev$
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing, Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
Modified: pkg/R/simulate-pomp.R
===================================================================
--- pkg/R/simulate-pomp.R 2012-03-09 16:00:59 UTC (rev 626)
+++ pkg/R/simulate-pomp.R 2012-03-15 17:17:34 UTC (rev 627)
@@ -2,7 +2,7 @@
simulate.internal <- function (object, nsim = 1, seed = NULL, params,
states = FALSE, obs = FALSE,
- times, t0, ...) {
+ times, t0, as.data.frame = FALSE, ...) {
if (missing(times))
times <- time(object,t0=FALSE)
@@ -14,6 +14,10 @@
else
t0 <- as.numeric(t0)
+ obs <- as.logical(obs)
+ states <- as.logical(states)
+ as.data.frame <- as.logical(as.data.frame)
+
if (missing(params))
params <- coef(object)
@@ -52,6 +56,66 @@
assign('.Random.seed',save.seed,envir=.GlobalEnv)
}
+ if (as.data.frame) {
+ if (obs && states) {
+ nsim <- ncol(retval$obs)
+ retval <- lapply(
+ seq_len(nsim),
+ function (k) {
+ nm <- rownames(retval$obs)
+ dm <- dim(retval$obs)
+ dim(retval$obs) <- c(dm[1],prod(dm[-1]))
+ rownames(retval$obs) <- nm
+ nm <- rownames(retval$states)
+ dm <- dim(retval$states)
+ dim(retval$states) <- c(dm[1],prod(dm[-1]))
+ rownames(retval$states) <- nm
+ cbind(
+ time=times,
+ as.data.frame(t(retval$obs)),
+ as.data.frame(t(retval$states)),
+ sim=as.integer(k)
+ )
+ }
+ )
+ retval <- do.call(rbind,retval)
+ } else if (obs || states) {
+ nsim <- ncol(retval)
+ retval <- lapply(
+ seq_len(nsim),
+ function (k) {
+ nm <- rownames(retval)
+ dm <- dim(retval)
+ dim(retval) <- c(dm[1],prod(dm[-1]))
+ rownames(retval) <- nm
+ cbind(
+ time=times,
+ as.data.frame(t(retval)),
+ sim=as.integer(k)
+ )
+ }
+ )
+ retval <- do.call(rbind,retval)
+ } else {
+ nsim <- length(retval)
+ if (nsim > 1) {
+ retval <- lapply(
+ seq_len(nsim),
+ function (k) {
+ x <- as.data.frame(retval[[k]])
+ x$sim <- as.integer(k)
+ x
+ }
+ )
+ retval <- do.call(rbind,retval)
+ } else {
+ retval <- as.data.frame(retval)
+ retval$sim <- 1
+ }
+ }
+
+ }
+
retval
}
Modified: pkg/R/trajectory-pomp.R
===================================================================
--- pkg/R/trajectory-pomp.R 2012-03-09 16:00:59 UTC (rev 626)
+++ pkg/R/trajectory-pomp.R 2012-03-15 17:17:34 UTC (rev 627)
@@ -3,14 +3,16 @@
setGeneric('trajectory')
-trajectory.internal <- function (object, params, times, t0, ...) {
+trajectory.internal <- function (object, params, times, t0, as.data.frame = FALSE, ...) {
if (missing(times))
times <- time(object,t0=FALSE)
else
times <- as.numeric(times)
- if (length(times)==0)
+ as.data.frame <- as.logical(as.data.frame)
+
+if (length(times)==0)
stop("if ",sQuote("times")," is empty, there is no work to do",call.=FALSE)
if (any(diff(times)<0))
@@ -21,7 +23,7 @@
else
t0 <- as.numeric(t0)
- if (t0>times[1])
+ if (t0>times[1])
stop("the zero-time ",sQuote("t0")," must occur no later than the first observation",call.=FALSE)
ntimes <- length(times)
@@ -62,7 +64,7 @@
if (type=="map") {
x <- .Call(iterate_map,object,times,t0,x0,params)
-
+
} else if (type=="vectorfield") {
skel <- get.pomp.fun(object at skeleton)
@@ -107,7 +109,7 @@
warning("abnormal exit from ODE integrator, istate = ",attr(X,'istate'),call.=FALSE)
x <- array(data=t(X[-1,-1]),dim=c(nvar,nrep,ntimes),dimnames=list(statenames,NULL,NULL))
-
+
if (length(znames)>0)
x[znames,,-1] <- apply(x[znames,,,drop=FALSE],c(1,2),diff)
@@ -117,6 +119,23 @@
}
+ if (as.data.frame) {
+ x <- lapply(
+ seq_len(ncol(x)),
+ function (k) {
+ nm <- rownames(x)
+ y <- x[,k,,drop=FALSE]
+ dim(y) <- dim(y)[c(1,3)]
+ y <- as.data.frame(t(y))
+ names(y) <- nm
+ y$time <- times
+ y$traj <- as.integer(k)
+ y
+ }
+ )
+ x <- do.call(rbind,x)
+ }
+
x
}
Modified: pkg/data/blowflies.rda
===================================================================
(Binary files differ)
Modified: pkg/data/dacca.rda
===================================================================
(Binary files differ)
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gompertz.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ricker.rda
===================================================================
(Binary files differ)
Modified: pkg/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/verhulst.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2012-03-09 16:00:59 UTC (rev 626)
+++ pkg/inst/NEWS 2012-03-15 17:17:34 UTC (rev 627)
@@ -1,4 +1,7 @@
NEWS
+0.40-9
+ o Setting the new argument 'as.data.frame' to 'TRUE' in 'simulate' and 'trajectory' causes the results to be returned as a data-frame.
+
0.40-8
o A new method 'partrans' allows transformation of vectors or matrices of parameters.
The parameter transformations have been pushed into C for speed.
Modified: pkg/inst/TODO
===================================================================
--- pkg/inst/TODO 2012-03-09 16:00:59 UTC (rev 626)
+++ pkg/inst/TODO 2012-03-15 17:17:34 UTC (rev 627)
@@ -1,11 +1,10 @@
-* as.data.frame option for 'simulate'
+* native C routine for vectorfield evaluation inside lsoda
+* parameter transformations: put 'transform' option into each estimation routine
+ (bsmc, mif, pmcmc, nlf, *.match)
+
* unit tests for 'sannbox'
-* Create objective functions for 'probe.match', 'traj.match', and
- 'spect.match' so that users can roll their own optimization
- routines. Use 'nloptr' as the target optimizer package?
-
* Improve 'pmcmc' to take advantage of all particles.
* Partial rejection control for 'pfilter'?
@@ -13,11 +12,10 @@
* Adaptive particle numbers in pfilter.
* It might be possible to make the writing of basic model functions
- using R expressions rather than functions.
+ using R expressions rather than functions....
-* Write 'plot' method for 'pfilterd.pomp' objects. This would display
- effective sample size, conditional log likelihood, filter means,
- prediction variances, prediction means?
+* Write 'plot' and 'print' methods for 'pfilterd.pomp' objects.
+ This would display effective sample size, conditional log likelihood, filter means, prediction variances, prediction means?
* Add plugin for adaptive tau leaping algorithm.
@@ -25,9 +23,10 @@
* Add LPA model examples.
+* Add BBS example.
+
* SDE examples.
* Extended Kalman filter.
* Plugin for compartmental models.
-
Deleted: pkg/inst/data-R/euler.sir.R
===================================================================
--- pkg/inst/data-R/euler.sir.R 2012-03-09 16:00:59 UTC (rev 626)
+++ pkg/inst/data-R/euler.sir.R 2012-03-15 17:17:34 UTC (rev 627)
@@ -1,65 +0,0 @@
-require(pomp)
-
-po <- pomp(
- data=data.frame(
- time=seq(from=1/52,to=4,by=1/52),
- reports=NA
- ),
- times="time",
- t0=0,
- rprocess=euler.sim(
- step.fun="_sir_euler_simulator",
- delta.t=1/52/20,
- PACKAGE="pomp"
- ),
- skeleton.type="vectorfield",
- skeleton="_sir_ODE",
- rmeasure="_sir_binom_rmeasure",
- dmeasure="_sir_binom_dmeasure",
- PACKAGE="pomp",
- obsnames = c("reports"),
- statenames=c("S","I","R","cases","W"),
- paramnames=c(
- "gamma","mu","iota",
- "beta1","beta.sd","pop","rho",
- "nbasis","degree","period"
- ),
- zeronames=c("cases"),
- comp.names=c("S","I","R"),
- to.log.transform=c(
- "gamma","mu","iota",
- "beta1","beta2","beta3","beta.sd",
- "rho",
- "S.0","I.0","R.0"
- ),
- parameter.transform=function (params, to.log.transform, ...) {
- params[to.log.transform] <- log(params[to.log.transform])
- params
- },
- parameter.inv.transform=function (params, to.log.transform, ...) {
- params[to.log.transform] <- exp(params[to.log.transform])
- params
- },
- initializer=function(params, t0, comp.names, ...) {
- ic.names <- paste(comp.names,".0",sep="")
- snames <- c("S","I","R","cases","W")
- fracs <- exp(params[ic.names])
- x0 <- numeric(length(snames))
- names(x0) <- snames
- x0[comp.names] <- round(params['pop']*fracs/sum(fracs))
- x0["cases"] <- 0
- x0
- }
- )
-coef(po,transform=TRUE) <- c(gamma=26,mu=0.02,iota=0.01,
- nbasis=3,degree=3,period=1,
- beta1=1200,beta2=1800,beta3=600,
- beta.sd=1e-3,
- pop=2.1e6,
- rho=0.6,
- S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200
- )
-
-simulate(po,nsim=1,seed=329348545L) -> euler.sir
-
-save(euler.sir,file="euler.sir.rda",compress="xz")
Deleted: pkg/inst/data-R/gillespie.sir.R
===================================================================
--- pkg/inst/data-R/gillespie.sir.R 2012-03-09 16:00:59 UTC (rev 626)
+++ pkg/inst/data-R/gillespie.sir.R 2012-03-15 17:17:34 UTC (rev 627)
@@ -1,91 +0,0 @@
-library(pomp)
-
-po <- pomp(
- data=data.frame(
- time=seq(from=0,to=10,by=1/52),
- reports=NA
- ),
- times="time",
- t0=0,
- rprocess=gillespie.sim(
- rate.fun="_sir_rates",
- PACKAGE="pomp",
- v=cbind(
- birth=c(1,0,0,1,0),
- sdeath=c(-1,0,0,-1,0),
- infection=c(-1,1,0,0,0),
- ideath=c(0,-1,0,-1,0),
- recovery=c(0,-1,1,0,1),
- rdeath=c(0,0,-1,-1,0)
- ),
- d=cbind(
- birth=c(0,0,0,1,0),
- sdeath=c(1,0,0,0,0),
- infection=c(1,1,0,1,0),
- ideath=c(0,1,0,0,0),
- recovery=c(0,1,0,0,0),
- rdeath=c(0,0,1,0,0)
- )
- ),
- obsnames="reports",
- statenames=c("S","I","R","N","cases"),
- paramnames=c(
- "gamma","mu","iota",
- "beta1","nu",
- "nbasis","degree","period"
- ),
- zeronames=c("cases"),
- measurement.model=reports~binom(size=cases,prob=exp(rho)),
- to.log.transform=c(
- "gamma","nu","mu","iota",
- "beta1","beta2","beta3",
- "rho",
- "S.0","I.0","R.0"
- ),
- parameter.transform=function (params, to.log.transform, ...) {
- params[to.log.transform] <- log(params[to.log.transform])
- params
- },
- parameter.inv.transform=function (params, to.log.transform, ...) {
- params[to.log.transform] <- exp(params[to.log.transform])
- params
- },
- initializer=function(params, t0, ...){
- comp.names <- c("S","I","R")
- icnames <- paste(comp.names,"0",sep=".")
- snames <- c("S","I","R","N","cases")
- fracs <- exp(params[icnames])
- x0 <- numeric(length(snames))
- names(x0) <- snames
- x0["N"] <- params["pop"]
- x0[comp.names] <- round(params["pop"]*fracs/sum(fracs))
- x0
- }
- )
-
-coef(po,transform=TRUE) <- c(
- gamma=24,
- nu=1/70,
- mu=1/70,
- iota=0.1,
- beta1=330,
- beta2=410,
- beta3=490,
- rho=0.1,
- S.0=0.05,
- I.0=1e-4,
- R.0=0.95,
- pop=1000000,
- nbasis=3,
- degree=3,
- period=1
- )
-
-
-simulate(
- po,
- nsim=1,
- seed=1165270654L
- ) -> gillespie.sir
-
-save(gillespie.sir,file="gillespie.sir.rda",compress="xz")
Copied: pkg/inst/data-R/sir.R (from rev 626, pkg/inst/data-R/euler.sir.R)
===================================================================
--- pkg/inst/data-R/sir.R (rev 0)
+++ pkg/inst/data-R/sir.R 2012-03-15 17:17:34 UTC (rev 627)
@@ -0,0 +1,154 @@
+require(pomp)
+
+po <- pomp(
+ data=data.frame(
+ time=seq(from=1/52,to=4,by=1/52),
+ reports=NA
+ ),
+ times="time",
+ t0=0,
+ rprocess=euler.sim(
+ step.fun="_sir_euler_simulator",
+ delta.t=1/52/20,
+ PACKAGE="pomp"
+ ),
+ skeleton.type="vectorfield",
+ skeleton="_sir_ODE",
+ rmeasure="_sir_binom_rmeasure",
+ dmeasure="_sir_binom_dmeasure",
+ PACKAGE="pomp",
+ obsnames = c("reports"),
+ statenames=c("S","I","R","cases","W"),
+ paramnames=c(
+ "gamma","mu","iota",
+ "beta1","beta.sd","pop","rho",
+ "nbasis","degree","period"
+ ),
+ zeronames=c("cases"),
+ comp.names=c("S","I","R"),
+ to.log.transform=c(
+ "gamma","mu","iota",
+ "beta1","beta2","beta3","beta.sd",
+ "rho",
+ "S.0","I.0","R.0"
+ ),
+ parameter.transform=function (params, to.log.transform, ...) {
+ params[to.log.transform] <- log(params[to.log.transform])
+ params
+ },
+ parameter.inv.transform=function (params, to.log.transform, ...) {
+ params[to.log.transform] <- exp(params[to.log.transform])
+ params
+ },
+ initializer=function(params, t0, comp.names, ...) {
+ ic.names <- paste(comp.names,".0",sep="")
+ snames <- c("S","I","R","cases","W")
+ fracs <- exp(params[ic.names])
+ x0 <- numeric(length(snames))
+ names(x0) <- snames
+ x0[comp.names] <- round(params['pop']*fracs/sum(fracs))
+ x0["cases"] <- 0
+ x0
+ }
+ )
+coef(po,transform=TRUE) <- c(gamma=26,mu=0.02,iota=0.01,
+ nbasis=3,degree=3,period=1,
+ beta1=1200,beta2=1800,beta3=600,
+ beta.sd=1e-3,
+ pop=2.1e6,
+ rho=0.6,
+ S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200
+ )
+
+simulate(po,nsim=1,seed=329348545L) -> euler.sir
+
+save(euler.sir,file="euler.sir.rda",compress="xz")
+
+po <- pomp(
+ data=data.frame(
+ time=seq(from=0,to=10,by=1/52),
+ reports=NA
+ ),
+ times="time",
+ t0=0,
+ rprocess=gillespie.sim(
+ rate.fun="_sir_rates",
+ PACKAGE="pomp",
+ v=cbind(
+ birth=c(1,0,0,1,0),
+ sdeath=c(-1,0,0,-1,0),
+ infection=c(-1,1,0,0,0),
+ ideath=c(0,-1,0,-1,0),
+ recovery=c(0,-1,1,0,1),
+ rdeath=c(0,0,-1,-1,0)
+ ),
+ d=cbind(
+ birth=c(0,0,0,1,0),
+ sdeath=c(1,0,0,0,0),
+ infection=c(1,1,0,1,0),
+ ideath=c(0,1,0,0,0),
+ recovery=c(0,1,0,0,0),
+ rdeath=c(0,0,1,0,0)
+ )
+ ),
+ obsnames="reports",
+ statenames=c("S","I","R","N","cases"),
+ paramnames=c(
+ "gamma","mu","iota",
+ "beta1","nu",
+ "nbasis","degree","period"
+ ),
+ zeronames=c("cases"),
+ measurement.model=reports~binom(size=cases,prob=exp(rho)),
+ to.log.transform=c(
+ "gamma","nu","mu","iota",
+ "beta1","beta2","beta3",
+ "rho",
+ "S.0","I.0","R.0"
+ ),
+ parameter.transform=function (params, to.log.transform, ...) {
+ params[to.log.transform] <- log(params[to.log.transform])
+ params
+ },
+ parameter.inv.transform=function (params, to.log.transform, ...) {
+ params[to.log.transform] <- exp(params[to.log.transform])
+ params
+ },
+ initializer=function(params, t0, ...){
+ comp.names <- c("S","I","R")
+ icnames <- paste(comp.names,"0",sep=".")
+ snames <- c("S","I","R","N","cases")
+ fracs <- exp(params[icnames])
+ x0 <- numeric(length(snames))
+ names(x0) <- snames
+ x0["N"] <- params["pop"]
+ x0[comp.names] <- round(params["pop"]*fracs/sum(fracs))
+ x0
+ }
+ )
+
+coef(po,transform=TRUE) <- c(
+ gamma=24,
+ nu=1/70,
+ mu=1/70,
+ iota=0.1,
+ beta1=330,
+ beta2=410,
+ beta3=490,
+ rho=0.1,
+ S.0=0.05,
+ I.0=1e-4,
+ R.0=0.95,
+ pop=1000000,
+ nbasis=3,
+ degree=3,
+ period=1
+ )
+
+simulate(
+ po,
+ nsim=1,
+ seed=1165270654L
+ ) -> gillespie.sir
+
+save(gillespie.sir,file="gillespie.sir.rda",compress="xz")
Modified: pkg/man/simulate-pomp.Rd
===================================================================
--- pkg/man/simulate-pomp.Rd 2012-03-09 16:00:59 UTC (rev 626)
+++ pkg/man/simulate-pomp.Rd 2012-03-15 17:17:34 UTC (rev 627)
@@ -8,7 +8,8 @@
}
\usage{
\S4method{simulate}{pomp}(object, nsim = 1, seed = NULL, params,
- states = FALSE, obs = FALSE, times, t0, \dots)
+ states = FALSE, obs = FALSE, times, t0,
+ as.data.frame = FALSE, \dots)
}
\arguments{
\item{object}{An object of class \code{pomp}.}
@@ -33,6 +34,9 @@
\code{t0} specifies the start time (the time at which the initial conditions hold).
The default for \code{times} is is \code{times=time(object,t0=FALSE)} and \code{t0=timezero(object)}, respectively.
}
+ \item{as.data.frame}{
+ logical; if \code{TRUE}, return the result as a data-frame.
+ }
\item{\dots}{further arguments that are currently ignored.}
}
\value{
Modified: pkg/man/sir.Rd
===================================================================
--- pkg/man/sir.Rd 2012-03-09 16:00:59 UTC (rev 626)
+++ pkg/man/sir.Rd 2012-03-15 17:17:34 UTC (rev 627)
@@ -1,4 +1,5 @@
\name{sir}
+\alias{sir}
\alias{euler.sir}
\alias{gillespie.sir}
\docType{data}
Modified: pkg/man/trajectory-pomp.Rd
===================================================================
--- pkg/man/trajectory-pomp.Rd 2012-03-09 16:00:59 UTC (rev 626)
+++ pkg/man/trajectory-pomp.Rd 2012-03-15 17:17:34 UTC (rev 627)
@@ -11,7 +11,7 @@
In the case of a continuous-time system, the deterministic skeleton is a vector-field; \code{trajectory} integrates the vectorfield to obtain a trajectory.
}
\usage{
-\S4method{trajectory}{pomp}(object, params, times, t0, \dots)
+\S4method{trajectory}{pomp}(object, params, times, t0, as.data.frame = FALSE, \dots)
}
\arguments{
\item{object}{an object of class \code{pomp}.}
@@ -24,6 +24,9 @@
\code{t0} specifies the start time (the time at which the initial conditions hold).
The default for \code{times} is \code{times=time(object,t0=FALSE)} and \code{t0=timezero(object)}, respectively.
}
+ \item{as.data.frame}{
+ logical; if \code{TRUE}, return the result as a data-frame.
+ }
\item{\dots}{
additional arguments are passed to the ODE integrator if the skeleton is a vectorfield and ignored if it is a map.
See \code{\link[deSolve]{ode}} for a description of the additional arguments accepted.
@@ -54,6 +57,8 @@
x <- trajectory(euler.sir)
plot(time(euler.sir),x["I",1,],type='l',xlab='time',ylab='I')
lines(time(euler.sir),x["cases",1,],col='red')
+
+x <- trajectory(euler.sir,as.data.frame=TRUE)
}
\seealso{\code{\link{pomp}}, \code{\link{traj.match}}, \code{\link[deSolve]{ode}}}
\keyword{models}
Added: pkg/tests/dacca.R
===================================================================
--- pkg/tests/dacca.R (rev 0)
+++ pkg/tests/dacca.R 2012-03-15 17:17:34 UTC (rev 627)
@@ -0,0 +1,13 @@
+library(pomp)
+
+set.seed(1420306530L)
+
+data(dacca)
+
+x <- as.data.frame(dacca)
+print(names(x))
+print(dim(x))
+
+x <- simulate(dacca,nsim=3,as.data.frame=TRUE)
+print(names(x))
+print(dim(x))
Added: pkg/tests/dacca.Rout.save
===================================================================
--- pkg/tests/dacca.Rout.save (rev 0)
+++ pkg/tests/dacca.Rout.save 2012-03-15 17:17:34 UTC (rev 627)
@@ -0,0 +1,46 @@
+
+R version 2.14.2 (2012-02-29)
+Copyright (C) 2012 The R Foundation for Statistical Computing
+ISBN 3-900051-07-0
+Platform: x86_64-unknown-linux-gnu (64-bit)
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+> library(pomp)
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: deSolve
+>
+> set.seed(1420306530L)
+>
+> data(dacca)
+>
+> x <- as.data.frame(dacca)
+> print(names(x))
+ [1] "time" "cholera.deaths" "seas.1" "seas.2"
+ [5] "seas.3" "seas.4" "seas.5" "seas.6"
+ [9] "pop" "dpopdt" "trend"
+> print(dim(x))
+[1] 600 11
+>
+> x <- simulate(dacca,nsim=3,as.data.frame=TRUE)
+> print(names(x))
+ [1] "time" "cholera.deaths" "S" "I"
+ [5] "Rs" "R1" "R2" "R3"
+ [9] "M" "W" "count" "seas.1"
+[13] "seas.2" "seas.3" "seas.4" "seas.5"
+[17] "seas.6" "pop" "dpopdt" "trend"
+[21] "sim"
+> print(dim(x))
+[1] 1800 21
+>
Modified: pkg/tests/ricker.R
===================================================================
--- pkg/tests/ricker.R 2012-03-09 16:00:59 UTC (rev 626)
+++ pkg/tests/ricker.R 2012-03-15 17:17:34 UTC (rev 627)
@@ -10,6 +10,35 @@
lines(30:50,tj.2[1,,],col='red',lwd=2)
max(abs(tj.1[,,time(ricker)>=30]-tj.2[,,]))
+tj.3 <- trajectory(ricker,as.data.frame=TRUE)
+plot(tj.3)
+tj.3 <- trajectory(ricker,as.data.frame=TRUE,params=parmat(coef(ricker),3),times=1:100)
+plot(N~time,data=tj.3,subset=traj==3,type='l')
+
+sm <- simulate(ricker,seed=343995,as.data.frame=TRUE)
+sm1 <- as.data.frame(simulate(ricker,seed=343995))
+stopifnot(identical(sm[names(sm1)],sm1))
+
+sm <- simulate(ricker,nsim=3,as.data.frame=TRUE)
+print(names(sm))
+print(dim(sm))
+
+sm <- simulate(ricker,nsim=3,obs=T,as.data.frame=TRUE)
+print(names(sm))
+print(dim(sm))
+
+sm <- simulate(ricker,nsim=3,states=T,as.data.frame=TRUE)
+print(names(sm))
+print(dim(sm))
+
+sm <- simulate(ricker,nsim=3,states=T,obs=T,as.data.frame=TRUE)
+print(names(sm))
+print(dim(sm))
+
+sm <- simulate(ricker,nsim=1,states=T,obs=T,as.data.frame=TRUE)
+print(names(sm))
+print(dim(sm))
+
po <- ricker
try(
coef(po,"r")
Modified: pkg/tests/ricker.Rout.save
===================================================================
--- pkg/tests/ricker.Rout.save 2012-03-09 16:00:59 UTC (rev 626)
+++ pkg/tests/ricker.Rout.save 2012-03-15 17:17:34 UTC (rev 627)
@@ -1,6 +1,6 @@
-R version 2.13.1 (2011-07-08)
-Copyright (C) 2011 The R Foundation for Statistical Computing
+R version 2.14.2 (2012-02-29)
+Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-linux-gnu (64-bit)
@@ -32,6 +32,45 @@
> max(abs(tj.1[,,time(ricker)>=30]-tj.2[,,]))
[1] 0
>
+> tj.3 <- trajectory(ricker,as.data.frame=TRUE)
+> plot(tj.3)
+> tj.3 <- trajectory(ricker,as.data.frame=TRUE,params=parmat(coef(ricker),3),times=1:100)
+> plot(N~time,data=tj.3,subset=traj==3,type='l')
+>
+> sm <- simulate(ricker,seed=343995,as.data.frame=TRUE)
+> sm1 <- as.data.frame(simulate(ricker,seed=343995))
+> stopifnot(identical(sm[names(sm1)],sm1))
+>
+> sm <- simulate(ricker,nsim=3,as.data.frame=TRUE)
+> print(names(sm))
+[1] "time" "y" "N" "e" "sim"
+> print(dim(sm))
+[1] 153 5
+>
+> sm <- simulate(ricker,nsim=3,obs=T,as.data.frame=TRUE)
+> print(names(sm))
+[1] "time" "y" "sim"
+> print(dim(sm))
+[1] 459 3
+>
+> sm <- simulate(ricker,nsim=3,states=T,as.data.frame=TRUE)
+> print(names(sm))
+[1] "time" "N" "e" "sim"
+> print(dim(sm))
+[1] 459 4
+>
+> sm <- simulate(ricker,nsim=3,states=T,obs=T,as.data.frame=TRUE)
+> print(names(sm))
+[1] "time" "y" "N" "e" "sim"
+> print(dim(sm))
+[1] 459 5
+>
+> sm <- simulate(ricker,nsim=1,states=T,obs=T,as.data.frame=TRUE)
+> print(names(sm))
+[1] "time" "y" "N" "e" "sim"
+> print(dim(sm))
+[1] 51 5
+>
> po <- ricker
> try(
+ coef(po,"r")
More information about the pomp-commits
mailing list