[Pomp-commits] r814 - in pkg/pompExamples: . R inst man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jan 13 14:30:50 CET 2013
Author: kingaa
Date: 2013-01-13 14:30:49 +0100 (Sun, 13 Jan 2013)
New Revision: 814
Added:
pkg/pompExamples/R/budmoth.R
pkg/pompExamples/R/pertussis.R
Removed:
pkg/pompExamples/inst/data-R/
Modified:
pkg/pompExamples/DESCRIPTION
pkg/pompExamples/NAMESPACE
pkg/pompExamples/inst/NEWS
pkg/pompExamples/man/budmoth.Rd
pkg/pompExamples/man/pertussis.Rd
pkg/pompExamples/tests/budmoth.R
pkg/pompExamples/tests/budmoth.Rout.save
pkg/pompExamples/tests/pertussis.R
pkg/pompExamples/tests/pertussis.Rout.save
Log:
- example pomp objects are now constructed by calls to a function rather than a call to 'data'
Modified: pkg/pompExamples/DESCRIPTION
===================================================================
--- pkg/pompExamples/DESCRIPTION 2013-01-11 19:03:31 UTC (rev 813)
+++ pkg/pompExamples/DESCRIPTION 2013-01-13 13:30:49 UTC (rev 814)
@@ -1,8 +1,8 @@
Package: pompExamples
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.20-8
-Date: 2013-01-11
+Version: 0.21-1
+Date: 2013-01-13
Author: NCEAS Working Group on Inference for Mechanistic Models: Aaron King, Steve Ellner, Bruce Kendall, Daniel C. Reuman, Matt Ferrari, Ed Ionides, Helen Wearing
Maintainer: Aaron A. King <kingaa at umich.edu>
Description: Inference methods for partially-observed Markov processes
@@ -11,4 +11,4 @@
LazyLoad: true
LazyData: false
BuildVignettes: no
-Collate: version.R
+Collate: version.R budmoth.R pertussis.R
Modified: pkg/pompExamples/NAMESPACE
===================================================================
--- pkg/pompExamples/NAMESPACE 2013-01-11 19:03:31 UTC (rev 813)
+++ pkg/pompExamples/NAMESPACE 2013-01-13 13:30:49 UTC (rev 814)
@@ -1,2 +1,3 @@
useDynLib(pompExamples)
importFrom(pomp)
+export(pertussis.sim,budmoth.sim)
Copied: pkg/pompExamples/R/budmoth.R (from rev 813, pkg/pompExamples/inst/data-R/budmoth.sim.R)
===================================================================
--- pkg/pompExamples/R/budmoth.R (rev 0)
+++ pkg/pompExamples/R/budmoth.R 2013-01-13 13:30:49 UTC (rev 814)
@@ -0,0 +1,99 @@
+budmoth.sim <- function (which) {
+ if (missing(which)) {
+ cat("available datasets:",
+ sQuote(c("food","para1","para2","tri")),"\n")
+ } else {
+ w <- as.character(substitute(which))
+ simulate(
+ pomp(
+ data=data.frame(
+ time=seq(from=0,to=60,by=1),
+ Qobs=NA,Nobs=NA,Sobs=NA
+ ),
+ time="time",
+ t0=-1,
+ params=switch(
+ w,
+ tri=c(
+ alpha=0.5, sig.alpha=0.1, gam=50, lambda=22, sig.lambda=0.25, g=0.08, delta=10,
+ a=1.7, sig.a=0.1, w=0.15, beta0=0, beta1=35, u=0.9,
+ sigQobs=0.03, sigNobs=0.5, sigSobs=0.1,
+ Q.0=0.96, N.0=0.02, S.0=0.22
+ ),
+ food=c(
+ alpha=0.5, sig.alpha=0.1, gam=20, lambda=5, sig.lambda=0.25, g=0.02, delta=10,
+ a=1, sig.a=0.1, w=0, beta0=0, beta1=35, u=0.9,
+ sigQobs=0.03, sigNobs=0.5, sigSobs=0.1,
+ Q.0=0.96, N.0=0.02, S.0=0.22
+ ),
+ para1=c(
+ alpha=0.5, sig.alpha=0.1, gam=50, lambda=22, sig.lambda=0.25, g=0.08, delta=0.5,
+ a=1.7, sig.a=0.1, w=0.15, beta0=0, beta1=35, u=0.9,
+ sigQobs=0.03, sigNobs=0.5, sigSobs=0.1,
+ Q.0=0.96, N.0=0.02, S.0=0.22
+ ),
+ para2=c(
+ alpha=0.5, sig.alpha=0.1, gam=50, lambda=10, sig.lambda=5, g=0.08, delta=0.5,
+ a=1.7, sig.a=1, w=0.15, beta0=0, beta1=35, u=0.9,
+ sigQobs=0.03, sigNobs=0.5, sigSobs=0.1,
+ Q.0=0.96, N.0=0.02, S.0=0.22
+ ),
+ stop("unrecognized dataset ",sQuote(w),call.=FALSE)
+ ),
+ rprocess=euler.sim(
+ step.fun="budmoth_map",
+ delta.t=1,
+ PACKAGE="pompExamples"
+ ),
+ dprocess=onestep.dens(
+ dens.fun="budmoth_density",
+ PACKAGE="pompExamples"
+ ),
+ rmeasure="budmoth_rmeasure",
+ dmeasure="budmoth_dmeasure",
+ skeleton.type="map",
+ skeleton="budmoth_skeleton",
+ PACKAGE="pompExamples",
+ paramnames=c(
+ "alpha","sig.alpha","gam","lambda","sig.lambda",
+ "g","delta","a","sig.a",
+ "w","beta0","beta1","u",
+ "sigQobs","sigNobs","sigSobs"
+ ),
+ statenames=c(
+ "Alpha","Lambda","A","Q","N","S"
+ ),
+ obsnames=c("Qobs","Nobs","Sobs"),
+ initializer=function (params, t0, ...) {
+ x <- c(params[c("Q.0","N.0","S.0")],c(0,0,0))
+ names(x) <- c("Q","N","S","Alpha","Lambda","A")
+ x
+ },
+ logitvar=c("alpha","Q.0","S.0","u"),
+ logvar=c(
+ "sig.alpha","gam","lambda","sig.lambda",
+ "g","delta","a","w","sig.a","beta1","sigQobs",
+ "sigNobs", "sigSobs","N.0"
+ ),
+ parameter.transform=function (params, logitvar, logvar, ...) {
+ params[logitvar] <- plogis(params[logitvar])
+ params[logvar] <- exp(params[logvar])
+ params
+ },
+ parameter.inv.transform=function (params, logitvar, logvar, ...) {
+ params[logitvar] <- qlogis(params[logitvar])
+ params[logvar] <- log(params[logvar])
+ params
+ }
+ ),
+ seed=switch(
+ w,
+ tri=1691699385L,
+ food=1054866677L,
+ para1=1116757478L,
+ para2=1361101458L,
+ stop("unrecognized dataset ",sQuote(w),call.=FALSE)
+ )
+ )
+ }
+}
Copied: pkg/pompExamples/R/pertussis.R (from rev 813, pkg/pompExamples/inst/data-R/pertussis.sim.R)
===================================================================
--- pkg/pompExamples/R/pertussis.R (rev 0)
+++ pkg/pompExamples/R/pertussis.R 2013-01-13 13:30:49 UTC (rev 814)
@@ -0,0 +1,126 @@
+## This file constructs a pomp object for one of several continuous time Markov chain models
+## for pertussis containing simulated data.
+## Simulations from the process model are approximated using an Euler approximation.
+
+pertussis.sim <- function (which) {
+ if (missing(which)) {
+ cat("available datasets:",
+ sQuote(c("SEIR.small","SEIR.big",
+ "SEIRS.small","SEIRS.big",
+ "SEIRR.small","SEIRR.big",
+ "full.small","full.big")),"\n")
+ } else {
+ which <- as.character(substitute(which))
+ simulate(
+ pomp(
+ data=data.frame(time=seq(from=0,to=20,by=1/52),reports=NA),
+ times="time",
+ t0=-1/52,
+ params=switch(
+ which,
+ SEIR.small=c(
+ birthrate=0.02, deathrate=0.02, mean.beta=450, ampl.beta=0.15,
+ imports=10, sigma=46, gamma=26, alpha=0, alpha.ratio=1,
+ report.prob=0.3, boost.prob=0, polar.prob=0, foi.mod=0,
+ popsize=5e+5, noise.sigma=0, tau=0.01,
+ S.0=0.0574148031949802, E.0=0.0004081763321755, I.0=0.00067028956509212,
+ R1.0=0.941506730907752, R2.0=0),
+ SEIR.big=c(
+ birthrate=0.02, deathrate=0.02, mean.beta=450, ampl.beta=0.15,
+ imports=10, sigma=46, gamma=26, alpha=0, alpha.ratio=1,
+ report.prob=0.3, boost.prob=0, polar.prob=0, foi.mod=0,
+ popsize=5e+6, noise.sigma=0, tau=0.01,
+ S.0=0.0515635231482973, E.0=0.000437143470487014, I.0=0.000734641109212043,
+ R1.0=0.947264692272004, R2.0=0),
+ SEIRS.small=c(
+ birthrate=0.02, deathrate=0.02, mean.beta=150, ampl.beta=0.15,
+ imports=10, sigma=46, gamma=26, alpha=0.1, alpha.ratio=1,
+ report.prob=0.1, boost.prob=0, polar.prob=0, foi.mod=0,
+ popsize=5e+5, noise.sigma=0, tau=0.01,
+ S.0=0.157360395940609, E.0=0.000837874318852172, I.0=0.00124181372794081,
+ R1.0=0.45913512973054, R2.0=0.381424786282058),
+ SEIRS.big=c(
+ birthrate=0.02, deathrate=0.02, mean.beta=150, ampl.beta=0.15,
+ imports=10, sigma=46, gamma=26, alpha=0.1, alpha.ratio=1,
+ report.prob=0.1, boost.prob=0, polar.prob=0, foi.mod=0,
+ popsize=5e+6, noise.sigma=0, tau=0.01,
+ S.0=0.157398354546347, E.0=0.00132093662562661, I.0=0.0022558671035406,
+ R1.0=0.457185201591761, R2.0=0.381839640132725),
+ SEIRR.small=c(
+ birthrate=0.02, deathrate=0.02, mean.beta=150, ampl.beta=0.15,
+ imports=10, sigma=46, gamma=26, alpha=0.1, alpha.ratio=1,
+ report.prob=0.11, boost.prob=0.75, polar.prob=0, foi.mod=0.5,
+ popsize=5e+5, noise.sigma=0, tau=0.01,
+ S.0=0.128943112158304, E.0=0.00068688724266688, I.0=0.00114414648269803,
+ R1.0=0.638074319602244, R2.0=0.231151534514087),
+ SEIRR.big=c(
+ birthrate=0.02, deathrate=0.02, mean.beta=150, ampl.beta=0.15,
+ imports=10, sigma=46, gamma=26, alpha=0.1, alpha.ratio=1,
+ report.prob=0.11, boost.prob=0.75, polar.prob=0, foi.mod=0.5,
+ popsize=5e+6, noise.sigma=0, tau=0.01,
+ S.0=0.127128689912424, E.0=0.00126497004491763, I.0=0.00216092385991776,
+ R1.0=0.639879739889535, R2.0=0.229565676293206),
+ full.small=c(
+ birthrate=0.02, deathrate=0.02, mean.beta=150, ampl.beta=0.15,
+ imports=10, sigma=46, gamma=26, alpha=0.1, alpha.ratio=1,
+ report.prob=0.1, boost.prob=0.75, polar.prob=0.1, foi.mod=0.5,
+ popsize=5e+5, noise.sigma=0.01, tau=0.01,
+ S.0=0.132553922599906, E.0=0.0010539075727066, I.0=0.00166100642162314,
+ R1.0=0.641737544956371, R2.0=0.222993618449393),
+ full.big=c(
+ birthrate=0.02, deathrate=0.02, mean.beta=150, ampl.beta=0.15,
+ imports=10, sigma=46, gamma=26, alpha=0.1, alpha.ratio=1,
+ report.prob=0.1, boost.prob=0.75, polar.prob=0.1, foi.mod=0.5,
+ popsize=5e+6, noise.sigma=0.01, tau=0.01,
+ S.0=0.130980596244438, E.0=0.00115096693013597, I.0=0.0018994251960431,
+ R1.0=0.643957103848235, R2.0=0.222011907781148),
+ stop("unrecognized dataset ",sQuote(which),call.=FALSE)
+ ),
+ rprocess = euler.sim(
+ step.fun="pertussis_sveirr_EM",
+ delta.t=1/52/7, # Euler stepsize
+ PACKAGE="pompExamples"
+ ),
+ skeleton.type="vectorfield",
+ skeleton="pertussis_sveirr_skel",
+ PACKAGE="pompExamples",
+ paramnames=c(
+ "birthrate","deathrate","mean.beta","ampl.beta",
+ "imports","sigma","gamma","alpha","alpha.ratio",
+ "report.prob","boost.prob","polar.prob",
+ "foi.mod","noise.sigma","popsize","tau",
+ "S.0","E.0","I.0","R1.0","R2.0"
+ ),
+ statenames=c("S","E","I","R1","R2","cases","W","err","simpop"),
+ zeronames=c("cases","err"),
+ ivps=c("S.0","E.0","I.0","R1.0","R2.0"),
+ comp.names=c("S","E","I","R1","R2"),
+ rmeasure = "negbin_rmeasure",
+ dmeasure = "negbin_dmeasure",
+ parameter.inv.transform="pertussis_par_untrans",
+ parameter.transform="pertussis_par_trans",
+ initializer = function (params, t0, statenames, comp.names, ivps, ...) {
+ states <- numeric(length(statenames))
+ names(states) <- statenames
+ ## translate fractions into initial conditions
+ frac <- params[ivps]
+ states[comp.names] <- round(params['popsize']*frac/sum(frac))
+ states["simpop"] <- params["popsize"]
+ states
+ }
+ ),
+ seed=switch(
+ which,
+ SEIR.small=1831650124L,
+ SEIR.big=908022490L,
+ SEIRS.small=1111340406L,
+ SEIRS.big=1751228386L,
+ SEIRR.small=350421545L,
+ SEIRR.big=748454784L,
+ full.small=581894515L,
+ full.big=301057392L,
+ stop("unrecognized dataset ",sQuote(which),call.=FALSE)
+ )
+ )
+ }
+}
Modified: pkg/pompExamples/inst/NEWS
===================================================================
--- pkg/pompExamples/inst/NEWS 2013-01-11 19:03:31 UTC (rev 813)
+++ pkg/pompExamples/inst/NEWS 2013-01-13 13:30:49 UTC (rev 814)
@@ -1,4 +1,8 @@
NEWS
+0.22-1
+ o The pomp objects are no longer accessed using 'data'.
+ Instead, 'pertussis.sim' and 'budmoth.sim' are functions that return the pomp objects.
+
0.21-1
o Changed t0 for the budmoth model from 0 to -1.
The data are different as well.
Modified: pkg/pompExamples/man/budmoth.Rd
===================================================================
--- pkg/pompExamples/man/budmoth.Rd 2013-01-11 19:03:31 UTC (rev 813)
+++ pkg/pompExamples/man/budmoth.Rd 2013-01-13 13:30:49 UTC (rev 814)
@@ -1,22 +1,26 @@
\name{budmoth}
\alias{budmoth.sim}
-\docType{data}
\title{Larch budmoth model POMPs with real and simulated data.}
\description{
- \code{budmoth.sim} is a list of \code{pomp} objects containing the larch budmoth model, budmoth density, parasitism rate, and food quality (needle-length) simulated data.
- Four parameter regimes are represented.
+ \code{budmoth.sim} constructs a \code{pomp} object containing the larch budmoth model and simulated budmoth density, parasitism rate, and food quality (needle-length) data.
+ Four datasets, representing four distinct parameter regimes, are avaiable.
}
\usage{
-data(budmoth.sim)
+budmoth.sim(which)
}
+\arguments{
+ \item{which}{
+ dataset to load given as a name or literal character string.
+ Evoked without an argument, \code{budmoth.sim} lists all available datasets.
+ }
+}
\examples{
-data(budmoth.sim)
-names(budmoth.sim)
+budmoth.sim() ## print a list of all available datasets
## three regimes, high and low noise regimes for parasitism and tritrophic
-plot(budmoth.sim[["food"]])
-plot(budmoth.sim[["para1"]])
-plot(budmoth.sim[["para2"]])
-plot(budmoth.sim[["tri"]])
+plot(budmoth.sim(food))
+plot(budmoth.sim(para1))
+plot(budmoth.sim(para2))
+plot(budmoth.sim(tri))
}
-\seealso{the vignettes}
+\seealso{the dQuote{budmoth-model} vignette}
\keyword{datasets}
Modified: pkg/pompExamples/man/pertussis.Rd
===================================================================
--- pkg/pompExamples/man/pertussis.Rd 2013-01-11 19:03:31 UTC (rev 813)
+++ pkg/pompExamples/man/pertussis.Rd 2013-01-13 13:30:49 UTC (rev 814)
@@ -1,21 +1,24 @@
\name{pertussis}
\alias{pertussis.sim}
-\docType{data}
\title{Pertussis models with simulated data.}
\description{
- \code{pertussis.sim} is a list of \code{pomp} objects:
- each encodes one of the pertussis models and contains simulated data.
+ \code{pertussis.sim} constructs a \code{pomp} object containing one of the pertussis models and simulated incidence data.
}
\usage{
-data(pertussis.sim)
+pertussis.sim(which)
}
+\arguments{
+ \item{which}{
+ dataset to load given as a name or literal character string.
+ Evoked without an argument, \code{pertussis.sim} lists all available datasets.
+ }
+}
\examples{
-data(pertussis.sim)
-names(pertussis.sim)
-
-plot(pertussis.sim[["SEIR.small"]])
-plot(pertussis.sim[["SEIRS.big"]])
-plot(pertussis.sim[["SEIRR.small"]])
+pertussis.sim()
+plot(pertussis.sim(SEIR.small))
+plot(pertussis.sim(SEIRS.big))
+plot(pertussis.sim(SEIRR.small))
+plot(pertussis.sim(full.big))
}
-\seealso{the vignettes}
+\seealso{the dQuote{pertussis-model} vignette}
\keyword{datasets}
Modified: pkg/pompExamples/tests/budmoth.R
===================================================================
--- pkg/pompExamples/tests/budmoth.R 2013-01-11 19:03:31 UTC (rev 813)
+++ pkg/pompExamples/tests/budmoth.R 2013-01-13 13:30:49 UTC (rev 814)
@@ -1,17 +1,19 @@
require(pompExamples)
-data(budmoth.sim)
+all <- c("food","para1","para2","tri")
-names(budmoth.sim)
-x <- lapply(budmoth.sim,as,"data.frame")
+sapply(all,function(n)eval(bquote(budmoth.sim(.(n))))) -> bm
+names(bm)
+x <- lapply(bm,as,"data.frame")
+
print(lapply(x,tail))
-y <- simulate(budmoth.sim$food,seed=3434996L,as.data.frame=TRUE)
+y <- simulate(budmoth.sim(food),seed=3434996L,as.data.frame=TRUE)
tail(y)
-z <- trajectory(budmoth.sim$tri,as.data.frame=TRUE)
+z <- trajectory(budmoth.sim(tri),as.data.frame=TRUE)
tail(z)
-pf <- pfilter(budmoth.sim$food,seed=34348885L,Np=1000)
+pf <- pfilter(budmoth.sim(food),seed=34348885L,Np=1000)
logLik(pf)
Modified: pkg/pompExamples/tests/budmoth.Rout.save
===================================================================
--- pkg/pompExamples/tests/budmoth.Rout.save 2013-01-11 19:03:31 UTC (rev 813)
+++ pkg/pompExamples/tests/budmoth.Rout.save 2013-01-13 13:30:49 UTC (rev 814)
@@ -1,5 +1,5 @@
-R version 2.14.2 (2012-02-29)
+R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-linux-gnu (64-bit)
@@ -23,28 +23,30 @@
Loading required package: subplex
Loading required package: deSolve
>
-> data(budmoth.sim)
+> all <- c("food","para1","para2","tri")
>
-> names(budmoth.sim)
-[1] "tri" "para1" "food" "para2"
-> x <- lapply(budmoth.sim,as,"data.frame")
+> sapply(all,function(n)eval(bquote(budmoth.sim(.(n))))) -> bm
>
+> names(bm)
+[1] "food" "para1" "para2" "tri"
+> x <- lapply(bm,as,"data.frame")
+>
> print(lapply(x,tail))
-$tri
+$food
time Qobs Nobs Sobs Q N S
-56 55 33.36428 5.9679756 3.009644e-06 0.9830185 17.7373560 2.949154e-06
-57 56 29.34035 68.3377210 7.255485e-05 0.8608373 80.4984855 8.681423e-05
-58 57 21.75336 1.2686716 8.552474e-03 0.6310693 0.7176525 9.751675e-03
-59 58 29.15848 0.3500365 1.047452e-02 0.8025905 0.3723865 1.256469e-02
-60 59 32.98213 1.6128525 6.718390e-03 0.8941107 1.0912736 7.410424e-03
-61 60 33.51183 13.0954506 1.173616e-02 0.9378258 7.6959192 1.147648e-02
- Alpha Lambda A
-56 0.4833379 21.87970 1.738908
-57 0.5010435 22.22843 1.640583
-58 0.5190030 22.45659 1.404177
-59 0.5165432 22.20828 1.810056
-60 0.5183574 22.01444 1.591250
-61 0.4828455 22.35338 1.429722
+56 55 28.90370 20.4806075 2.870679e-05 0.8502012 19.7570009 0.0000257055
+57 56 25.04577 17.2422706 4.340917e-04 0.6794912 14.4355348 0.0005305078
+58 57 21.68427 2.1931554 7.856243e-03 0.6288764 1.9239332 0.0079677582
+59 58 27.25406 0.3103408 1.348361e-02 0.7899212 0.1973302 0.0141600465
+60 59 32.38417 0.1216027 3.499077e-03 0.8898312 0.1200628 0.0032722573
+61 60 32.80127 0.1632687 3.512564e-04 0.9346055 0.1908352 0.0003751743
+ Alpha Lambda A
+56 0.5413771 4.932935 0.8677828
+57 0.5082460 4.851747 1.0428964
+58 0.4871643 4.388592 1.0444642
+59 0.4316777 4.394991 0.9302542
+60 0.5012195 5.063878 1.1726507
+61 0.5703106 4.810311 0.9525775
$para1
time Qobs Nobs Sobs Q N S
@@ -62,22 +64,6 @@
60 0.5157446 21.68244 1.551352
61 0.4964261 21.93769 1.721177
-$food
- time Qobs Nobs Sobs Q N S
-56 55 28.90370 20.4806075 2.870679e-05 0.8502012 19.7570009 0.0000257055
-57 56 25.04577 17.2422706 4.340917e-04 0.6794912 14.4355348 0.0005305078
-58 57 21.68427 2.1931554 7.856243e-03 0.6288764 1.9239332 0.0079677582
-59 58 27.25406 0.3103408 1.348361e-02 0.7899212 0.1973302 0.0141600465
-60 59 32.38417 0.1216027 3.499077e-03 0.8898312 0.1200628 0.0032722573
-61 60 32.80127 0.1632687 3.512564e-04 0.9346055 0.1908352 0.0003751743
- Alpha Lambda A
-56 0.5413771 4.932935 0.8677828
-57 0.5082460 4.851747 1.0428964
-58 0.4871643 4.388592 1.0444642
-59 0.4316777 4.394991 0.9302542
-60 0.5012195 5.063878 1.1726507
-61 0.5703106 4.810311 0.9525775
-
$para2
time Qobs Nobs Sobs Q N S
56 55 33.87640 15.747819 2.248998e-06 0.9749960 6.921067 2.575005e-06
@@ -94,8 +80,24 @@
60 0.5059549 7.138668 6.7705675
61 0.5120129 10.895242 3.1116150
+$tri
+ time Qobs Nobs Sobs Q N S
+56 55 33.36428 5.9679756 3.009644e-06 0.9830185 17.7373560 2.949154e-06
+57 56 29.34035 68.3377210 7.255485e-05 0.8608373 80.4984855 8.681423e-05
+58 57 21.75336 1.2686716 8.552474e-03 0.6310693 0.7176525 9.751675e-03
+59 58 29.15848 0.3500365 1.047452e-02 0.8025905 0.3723865 1.256469e-02
+60 59 32.98213 1.6128525 6.718390e-03 0.8941107 1.0912736 7.410424e-03
+61 60 33.51183 13.0954506 1.173616e-02 0.9378258 7.6959192 1.147648e-02
+ Alpha Lambda A
+56 0.4833379 21.87970 1.738908
+57 0.5010435 22.22843 1.640583
+58 0.5190030 22.45659 1.404177
+59 0.5165432 22.20828 1.810056
+60 0.5183574 22.01444 1.591250
+61 0.4828455 22.35338 1.429722
+
>
-> y <- simulate(budmoth.sim$food,seed=3434996L,as.data.frame=TRUE)
+> y <- simulate(budmoth.sim(food),seed=3434996L,as.data.frame=TRUE)
> tail(y)
time Qobs Nobs Sobs Q N S
56 55 24.75707 1.2571930 0.0536837100 0.6960924 0.6909030 0.0655207892
@@ -112,7 +114,7 @@
60 0.4994891 5.038971 1.0048781 1
61 0.4818151 4.973050 0.9217527 1
>
-> z <- trajectory(budmoth.sim$tri,as.data.frame=TRUE)
+> z <- trajectory(budmoth.sim(tri),as.data.frame=TRUE)
> tail(z)
Q N S Alpha Lambda A time traj
56 0.9795835 16.9946885 0.0001199655 0.5 22 1.7 55 1
@@ -122,7 +124,10 @@
60 0.8998638 0.5133800 0.1500321097 0.5 22 1.7 59 1
61 0.9448503 3.3848665 0.1205147616 0.5 22 1.7 60 1
>
-> pf <- pfilter(budmoth.sim$food,seed=34348885L,Np=1000)
+> pf <- pfilter(budmoth.sim(food),seed=34348885L,Np=1000)
> logLik(pf)
[1] 360.1747
>
+> proc.time()
+ user system elapsed
+ 1.284 0.056 1.323
Modified: pkg/pompExamples/tests/pertussis.R
===================================================================
--- pkg/pompExamples/tests/pertussis.R 2013-01-11 19:03:31 UTC (rev 813)
+++ pkg/pompExamples/tests/pertussis.R 2013-01-13 13:30:49 UTC (rev 814)
@@ -1,19 +1,21 @@
require(pompExamples)
-data(pertussis.sim)
+all <- c("SEIR.small","SEIR.big","SEIRS.small","SEIRS.big","SEIRR.small","SEIRR.big","full.small","full.big")
-names(pertussis.sim)
-x <- lapply(pertussis.sim,as.data.frame)
+sapply(all,function(n)eval(bquote(pertussis.sim(.(n))))) -> pt
+names(pt)
+x <- lapply(pt,as.data.frame)
+
print(lapply(x,tail))
-x <- simulate(pertussis.sim$full.big,seed=395885L,as.data.frame=TRUE)
+x <- simulate(pertussis.sim(full.big),seed=395885L,as.data.frame=TRUE)
tail(x)
-y <- trajectory(pertussis.sim$SEIRS.small,as.data.frame=TRUE)
+y <- trajectory(pertussis.sim(SEIRS.small),as.data.frame=TRUE)
tail(y)
-system.time(pf <- pfilter(pertussis.sim$full.small,seed=3445886L,Np=1000))
+system.time(pf <- pfilter(pertussis.sim(full.small),seed=3445886L,Np=1000))
logLik(pf)
pttest <- function (po, digits = 15) {
@@ -23,7 +25,7 @@
)
}
-stopifnot(all(sapply(pertussis.sim,pttest)))
+stopifnot(all(sapply(pt,pttest)))
pttest <- function (po, digits = 15) {
identical(
@@ -32,5 +34,4 @@
)
}
-stopifnot(all(sapply(pertussis.sim,pttest)))
-
+stopifnot(all(sapply(pt,pttest)))
Modified: pkg/pompExamples/tests/pertussis.Rout.save
===================================================================
--- pkg/pompExamples/tests/pertussis.Rout.save 2013-01-11 19:03:31 UTC (rev 813)
+++ pkg/pompExamples/tests/pertussis.Rout.save 2013-01-13 13:30:49 UTC (rev 814)
@@ -1,5 +1,5 @@
-R version 2.15.0 (2012-03-30)
+R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-linux-gnu (64-bit)
@@ -23,12 +23,14 @@
Loading required package: subplex
Loading required package: deSolve
>
-> data(pertussis.sim)
+> all <- c("SEIR.small","SEIR.big","SEIRS.small","SEIRS.big","SEIRR.small","SEIRR.big","full.small","full.big")
>
-> names(pertussis.sim)
+> sapply(all,function(n)eval(bquote(pertussis.sim(.(n))))) -> pt
+>
+> names(pt)
[1] "SEIR.small" "SEIR.big" "SEIRS.small" "SEIRS.big" "SEIRR.small"
[6] "SEIRR.big" "full.small" "full.big"
-> x <- lapply(pertussis.sim,as.data.frame)
+> x <- lapply(pt,as.data.frame)
>
> print(lapply(x,tail))
$SEIR.small
@@ -111,7 +113,7 @@
1041 5002948
>
-> x <- simulate(pertussis.sim$full.big,seed=395885L,as.data.frame=TRUE)
+> x <- simulate(pertussis.sim(full.big),seed=395885L,as.data.frame=TRUE)
> tail(x)
time reports S E I R1 R2 cases W err
1036 19.90385 432 662434 5677 9458 3201351 1118149 4831 -0.4801610 7
@@ -128,7 +130,7 @@
1040 4997185 1
1041 4997213 1
>
-> y <- trajectory(pertussis.sim$SEIRS.small,as.data.frame=TRUE)
+> y <- trajectory(pertussis.sim(SEIRS.small),as.data.frame=TRUE)
> tail(y)
S E I R1 R2 cases W err simpop
1036 81409.73 558.4599 942.3100 227353.0 189736.5 487.2963 0 0 5e+05
@@ -145,9 +147,9 @@
1040 19.98077 1
1041 20.00000 1
>
-> system.time(pf <- pfilter(pertussis.sim$full.small,seed=3445886L,Np=1000))
+> system.time(pf <- pfilter(pertussis.sim(full.small),seed=3445886L,Np=1000))
user system elapsed
- 22.637 0.004 22.720
+ 65.132 0.216 65.495
> logLik(pf)
[1] -3829.33
>
@@ -158,7 +160,7 @@
+ )
+ }
>
-> stopifnot(all(sapply(pertussis.sim,pttest)))
+> stopifnot(all(sapply(pt,pttest)))
>
> pttest <- function (po, digits = 15) {
+ identical(
@@ -167,9 +169,8 @@
+ )
+ }
>
-> stopifnot(all(sapply(pertussis.sim,pttest)))
+> stopifnot(all(sapply(pt,pttest)))
>
->
> proc.time()
user system elapsed
- 23.089 0.020 23.200
+ 66.960 0.264 67.348
More information about the pomp-commits
mailing list