[Pomp-commits] r1023 - in pkg/pompExamples: . R inst/examples man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 17 20:12:21 CET 2014
Author: kingaa
Date: 2014-12-17 20:12:21 +0100 (Wed, 17 Dec 2014)
New Revision: 1023
Added:
pkg/pompExamples/inst/examples/budmoth.R
Removed:
pkg/pompExamples/R/budmoth.R
Modified:
pkg/pompExamples/DESCRIPTION
pkg/pompExamples/NAMESPACE
pkg/pompExamples/man/budmoth.Rd
pkg/pompExamples/tests/budmoth.R
pkg/pompExamples/tests/budmoth.Rout.save
pkg/pompExamples/tests/examples.R
Log:
- implement budmoth example for 'pompExample'
Modified: pkg/pompExamples/DESCRIPTION
===================================================================
--- pkg/pompExamples/DESCRIPTION 2014-12-17 19:12:15 UTC (rev 1022)
+++ pkg/pompExamples/DESCRIPTION 2014-12-17 19:12:21 UTC (rev 1023)
@@ -21,4 +21,4 @@
License: GPL (>= 2)
LazyData: false
BuildVignettes: true
-Collate: aaa.R budmoth.R pertussis.R
+Collate: aaa.R pertussis.R
Modified: pkg/pompExamples/NAMESPACE
===================================================================
--- pkg/pompExamples/NAMESPACE 2014-12-17 19:12:15 UTC (rev 1022)
+++ pkg/pompExamples/NAMESPACE 2014-12-17 19:12:21 UTC (rev 1023)
@@ -1,3 +1,3 @@
useDynLib(pompExamples)
import(pomp)
-export(pertussis.sim,budmoth.sim)
+export(pertussis.sim)
Deleted: pkg/pompExamples/R/budmoth.R
===================================================================
--- pkg/pompExamples/R/budmoth.R 2014-12-17 19:12:15 UTC (rev 1022)
+++ pkg/pompExamples/R/budmoth.R 2014-12-17 19:12:21 UTC (rev 1023)
@@ -1,106 +0,0 @@
-budmoth.sim <- function (which) {
- if (missing(which)) {
- datasets <- c("food","para1","para2","tri")
- cat("available datasets:",sQuote(datasets),"\n")
- invisible(datasets)
- } else {
- which <- 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(
- which,
- 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(which),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(
- which,
- tri=1691699385L,
- food=1054866677L,
- para1=1116757478L,
- para2=1361101458L,
- stop("unrecognized dataset ",sQuote(which),call.=FALSE)
- )
- )
- }
-}
Added: pkg/pompExamples/inst/examples/budmoth.R
===================================================================
--- pkg/pompExamples/inst/examples/budmoth.R (rev 0)
+++ pkg/pompExamples/inst/examples/budmoth.R 2014-12-17 19:12:21 UTC (rev 1023)
@@ -0,0 +1,113 @@
+require(pomp)
+
+budmoth.example <- function (which = c("food","para1","para2","tri")) {
+ which <- match.arg(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(
+ which,
+ 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(which),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(
+ which,
+ tri=1691699385L,
+ food=1054866677L,
+ para1=1116757478L,
+ para2=1361101458L,
+ stop("unrecognized dataset ",sQuote(which),call.=FALSE)
+ )
+ )
+}
+
+if (exists("which",where=environment(),inherits=FALSE)) {
+ budmoth.example(which) -> budmoth
+ c("budmoth")
+} else {
+ budmoth.example(which="food") -> food
+ budmoth.example(which="para1") -> para1
+ budmoth.example(which="para2") -> para2
+ budmoth.example(which="tri") -> tri
+ c("food","para1","para2","tri")
+}
Modified: pkg/pompExamples/man/budmoth.Rd
===================================================================
--- pkg/pompExamples/man/budmoth.Rd 2014-12-17 19:12:15 UTC (rev 1022)
+++ pkg/pompExamples/man/budmoth.Rd 2014-12-17 19:12:21 UTC (rev 1023)
@@ -1,26 +1,17 @@
\name{budmoth}
+\docType{data}
\alias{budmoth.sim}
\title{Larch budmoth model POMPs with real and simulated data.}
\description{
- \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.
+ \code{pompExample(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{
-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{
-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"))
+bm <- pompExample(budmoth,envir=NULL)
+plot(bm$food)
+plot(bm$para1)
+plot(bm$para2)
+plot(bm$tri)
}
-\seealso{the \dQuote{budmoth-model} vignette}
\keyword{datasets}
Modified: pkg/pompExamples/tests/budmoth.R
===================================================================
--- pkg/pompExamples/tests/budmoth.R 2014-12-17 19:12:15 UTC (rev 1022)
+++ pkg/pompExamples/tests/budmoth.R 2014-12-17 19:12:21 UTC (rev 1023)
@@ -2,18 +2,18 @@
all <- c("food","para1","para2","tri")
-sapply(all,function(n)eval(bquote(budmoth.sim(.(n))))) -> bm
+bm <- pompExample(budmoth,envir=NULL)
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(bm$food,seed=3434996L,as.data.frame=TRUE)
tail(y)
-z <- trajectory(budmoth.sim(tri),as.data.frame=TRUE)
+z <- trajectory(bm$tri,as.data.frame=TRUE)
tail(z)
-pf <- pfilter(budmoth.sim(food),seed=34348885L,Np=1000)
+pf <- pfilter(bm$para1,seed=34348885L,Np=1000)
logLik(pf)
Modified: pkg/pompExamples/tests/budmoth.Rout.save
===================================================================
--- pkg/pompExamples/tests/budmoth.Rout.save 2014-12-17 19:12:15 UTC (rev 1022)
+++ pkg/pompExamples/tests/budmoth.Rout.save 2014-12-17 19:12:21 UTC (rev 1023)
@@ -1,6 +1,6 @@
-R version 3.0.1 (2013-05-16) -- "Good Sport"
-Copyright (C) 2013 The R Foundation for Statistical Computing
+R version 3.1.2 (2014-10-31) -- "Pumpkin Helmet"
+Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -17,13 +17,12 @@
> library(pompExamples)
Loading required package: pomp
-Loading required package: mvtnorm
Loading required package: subplex
-Loading required package: deSolve
+Loading required package: nloptr
>
> all <- c("food","para1","para2","tri")
>
-> sapply(all,function(n)eval(bquote(budmoth.sim(.(n))))) -> bm
+> bm <- pompExample(budmoth,envir=NULL)
>
> names(bm)
[1] "food" "para1" "para2" "tri"
@@ -95,7 +94,7 @@
61 0.4828455 22.35338 1.429722
>
-> y <- simulate(budmoth.sim(food),seed=3434996L,as.data.frame=TRUE)
+> y <- simulate(bm$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 +111,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(bm$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,10 +121,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(bm$para1,seed=34348885L,Np=1000)
> logLik(pf)
-[1] 360.1747
+[1] 10.68836
>
> proc.time()
user system elapsed
- 0.604 0.064 0.689
+ 0.649 0.033 0.670
Modified: pkg/pompExamples/tests/examples.R
===================================================================
--- pkg/pompExamples/tests/examples.R 2014-12-17 19:12:15 UTC (rev 1022)
+++ pkg/pompExamples/tests/examples.R 2014-12-17 19:12:21 UTC (rev 1023)
@@ -23,6 +23,6 @@
pf <- pfilter(simulate(po$parus),Np=100,max.fail=Inf)
po <- pompExample(bbp,envir=NULL)
-pf <- pfilter(simulate(pf$bbp),Np=100,max.fail=Inf)
+pf <- pfilter(simulate(po$bbp),Np=100,max.fail=Inf)
## dev.off()
More information about the pomp-commits
mailing list