[Pomp-commits] r667 - in pkg/pompExamples: . data inst/data-R inst/doc tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 18 12:01:33 CEST 2012
Author: kingaa
Date: 2012-04-18 12:01:33 +0200 (Wed, 18 Apr 2012)
New Revision: 667
Added:
pkg/pompExamples/inst/data-R/budmoth-params.rda
pkg/pompExamples/inst/data-R/pertussis-params.rda
pkg/pompExamples/inst/doc/pertussis-model-true-loglik.rda
Removed:
pkg/pompExamples/inst/data-R/budmoth-simdata.rda
pkg/pompExamples/inst/data-R/pertussis-simdata.rda
Modified:
pkg/pompExamples/.Rinstignore
pkg/pompExamples/DESCRIPTION
pkg/pompExamples/data/budmoth.sim.rda
pkg/pompExamples/data/pertussis.sim.rda
pkg/pompExamples/inst/data-R/budmoth.sim.R
pkg/pompExamples/inst/data-R/pertussis.sim.R
pkg/pompExamples/inst/doc/budmoth-model-slices.rda
pkg/pompExamples/inst/doc/budmoth-model-true-loglik.rda
pkg/pompExamples/inst/doc/budmoth-model.Rnw
pkg/pompExamples/inst/doc/budmoth-model.pdf
pkg/pompExamples/inst/doc/pertussis-model.Rnw
pkg/pompExamples/inst/doc/pertussis-model.pdf
pkg/pompExamples/tests/budmoth.Rout.save
pkg/pompExamples/tests/pertussis.R
pkg/pompExamples/tests/pertussis.Rout.save
Log:
- new datasets
- revert dependence to R>2.14.2
- add log likelihood calculation to pertussis vignette
- true parameters are included with pomps
Modified: pkg/pompExamples/.Rinstignore
===================================================================
--- pkg/pompExamples/.Rinstignore 2012-04-17 21:37:50 UTC (rev 666)
+++ pkg/pompExamples/.Rinstignore 2012-04-18 10:01:33 UTC (rev 667)
@@ -1,3 +1,4 @@
inst/doc/Makefile
inst/data-R/Makefile
inst/doc/(.+?)\.bst$
+inst/doc/(.+?)\.rda$
Modified: pkg/pompExamples/DESCRIPTION
===================================================================
--- pkg/pompExamples/DESCRIPTION 2012-04-17 21:37:50 UTC (rev 666)
+++ pkg/pompExamples/DESCRIPTION 2012-04-18 10:01:33 UTC (rev 667)
@@ -1,12 +1,12 @@
Package: pompExamples
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.19-1
-Date: 2012-04-17
+Version: 0.20-1
+Date: 2012-04-18
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
-Depends: R(>= 2.15.0), stats, methods, graphics, pomp(>= 0.41-4)
+Depends: R(>= 2.14.2), stats, methods, graphics, pomp(>= 0.41-3)
License: GPL (>= 2)
LazyLoad: true
LazyData: false
Modified: pkg/pompExamples/data/budmoth.sim.rda
===================================================================
(Binary files differ)
Modified: pkg/pompExamples/data/pertussis.sim.rda
===================================================================
(Binary files differ)
Added: pkg/pompExamples/inst/data-R/budmoth-params.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pompExamples/inst/data-R/budmoth-params.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Deleted: pkg/pompExamples/inst/data-R/budmoth-simdata.rda
===================================================================
(Binary files differ)
Modified: pkg/pompExamples/inst/data-R/budmoth.sim.R
===================================================================
--- pkg/pompExamples/inst/data-R/budmoth.sim.R 2012-04-17 21:37:50 UTC (rev 666)
+++ pkg/pompExamples/inst/data-R/budmoth.sim.R 2012-04-18 10:01:33 UTC (rev 667)
@@ -1,14 +1,18 @@
require(pompExamples)
-load("budmoth-simdata.rda")
+load("budmoth-params.rda")
datasets <- rownames(params)
+seeds <- params["seed"]
params <- subset(params,select=-seed)
budmoth.sim <- list()
for (d in datasets) {
po <- pomp(
- data=subset(simdata,dataset==d,select=c(time,Qobs,Nobs,Sobs)),
+ data=data.frame(
+ time=seq(from=0,to=60,by=1),
+ Qobs=NA,Nobs=NA,Sobs=NA
+ ),
time="time",
t0=0,
rprocess=euler.sim(
@@ -58,7 +62,7 @@
}
)
coef(po) <- unlist(params[d,])
- budmoth.sim[[d]] <- po
+ budmoth.sim[[d]] <- simulate(po,seed=seeds[d,])
}
save(budmoth.sim,file="budmoth.sim.rda",compress="xz")
Added: pkg/pompExamples/inst/data-R/pertussis-params.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pompExamples/inst/data-R/pertussis-params.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Deleted: pkg/pompExamples/inst/data-R/pertussis-simdata.rda
===================================================================
(Binary files differ)
Modified: pkg/pompExamples/inst/data-R/pertussis.sim.R
===================================================================
--- pkg/pompExamples/inst/data-R/pertussis.sim.R 2012-04-17 21:37:50 UTC (rev 666)
+++ pkg/pompExamples/inst/data-R/pertussis.sim.R 2012-04-18 10:01:33 UTC (rev 667)
@@ -5,15 +5,16 @@
require(pompExamples)
-load("pertussis-simdata.rda")
+load("pertussis-params.rda")
datasets <- rownames(params)
-params <- subset(params,select=-c(model,pop,seed))
+seeds <- params["seed"]
+params <- subset(params,select=-seed)
pertussis.sim <- list()
for (d in datasets) {
po <- pomp(
- data=subset(simdata,dataset==d,select=c(time,reports)),
+ data=data.frame(time=seq(from=0,to=20,by=1/52),reports=NA),
times="time",
t0=-1/52,
rprocess = euler.sim(
@@ -65,7 +66,7 @@
}
)
coef(po) <- unlist(params[d,])
- pertussis.sim[[d]] <- po
+ pertussis.sim[[d]] <- simulate(po,seed=seeds[d,])
}
save(pertussis.sim,file="pertussis.sim.rda",compress="xz")
Modified: pkg/pompExamples/inst/doc/budmoth-model-slices.rda
===================================================================
(Binary files differ)
Modified: pkg/pompExamples/inst/doc/budmoth-model-true-loglik.rda
===================================================================
(Binary files differ)
Modified: pkg/pompExamples/inst/doc/budmoth-model.Rnw
===================================================================
--- pkg/pompExamples/inst/doc/budmoth-model.Rnw 2012-04-17 21:37:50 UTC (rev 666)
+++ pkg/pompExamples/inst/doc/budmoth-model.Rnw 2012-04-18 10:01:33 UTC (rev 667)
@@ -165,6 +165,7 @@
x <- ldply(budmoth.sim,as.data.frame)
x <- rename(x,c(.id="dataset"))
x <- melt(x,id.var=c("dataset","time"))
+x <- subset(x,variable%in%c("Qobs","Nobs","Sobs"))
pl <- ggplot(data=x,mapping=aes(x=time,y=value))+
geom_line()+facet_grid(variable~dataset,scale="free")
print(pl)
@@ -271,7 +272,7 @@
xtable(
loglik.truth,
caption=paste(
- "Estimated log likelihood at the true parameters for the first ",length(time(budmoth.sim[[1]]))," years of simulated data. ",
+ "Estimated log likelihood at the true parameters for the simulated budmoth data. ",
"To obtain these,",nrep,"particle filtering runs, each with",Np,"particles, were used. ",
"The column labeled",dQuote("se"),"gives the standard error of the Monte Carlo likelihood calculation. ",
"The computation took ",signif(etime,2),"~CPU~",units(etime),"on inexpensive processors. ",
@@ -308,8 +309,7 @@
estnames <- c("gam","lambda","g","delta","a","w")
par.range <- t(apply(true.pars[estnames,],1,function(x)c(0.5*min(x),1.5*max(x))))
-##ncpus <- as.integer(Sys.getenv("PBS_NP"))
-ncpus <- 6
+ncpus <- as.integer(Sys.getenv("PBS_NP"))
slices <- lapply(
budmoth.sim,
Modified: pkg/pompExamples/inst/doc/budmoth-model.pdf
===================================================================
(Binary files differ)
Added: pkg/pompExamples/inst/doc/pertussis-model-true-loglik.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pompExamples/inst/doc/pertussis-model-true-loglik.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Modified: pkg/pompExamples/inst/doc/pertussis-model.Rnw
===================================================================
--- pkg/pompExamples/inst/doc/pertussis-model.Rnw 2012-04-17 21:37:50 UTC (rev 666)
+++ pkg/pompExamples/inst/doc/pertussis-model.Rnw 2012-04-18 10:01:33 UTC (rev 667)
@@ -374,85 +374,156 @@
)
@
-Plots of the data are shown in Figs.~\ref{fig:seir-plots}--\ref{fig:full-plots}.
+Plots of the data are shown in Fig.~\ref{fig:dataplots}.
-<<echo=F>>=
-tm <- time(pertussis.sim[[1]])
-ps <- sapply(pertussis.sim,coef,pars="popsize")
-cases <- t(apply(sapply(pertussis.sim,data.array),1,function(x)1e5*x/ps))
-@
-
-<<dataplot,echo=F,eval=F>>=
-op <- par(mfrow=c(2,1),mar=c(3,3,1,1),mgp=c(2,1,0))
-matplot(tm,cases[,ds],type='l',bty='l',lty=1,xlab="time",ylab=expression(cases/10^5))
-matplot(tm[tm<=20],cases[tm<=20,ds],type='l',bty='l',lty=1,xlab="time",ylab=expression(cases/10^5))
-par(op)
-@
-
\begin{figure}[h]
\begin{center}
-<<SEIR-plots,fig=T,echo=F>>=
-ds <- c("SEIR.small","SEIR.big")
-<<dataplot>>
-@
- \end{center}
- \caption{
- Plots of the simulated data, SEIR model with initial conditions chosen to give stationary dynamics.
- Population size \scinot{5}{5} (black); \scinot{5}{6} (red).
- Top: the full data sets.
- Bottom: the first 20~yr of simulated data.
- }
- \label{fig:seir-plots}
-\end{figure}
+<<dataplots,fig=T,echo=F,png=T,pdf=F,height=8,width=6>>=
+require(ggplot2)
+require(plyr)
+require(reshape)
+require(pompExamples)
+data(pertussis.sim)
+x <- ldply(pertussis.sim,as.data.frame)
+x <- rename(x,c(.id="dataset"))
+x <- melt(x,id.var=c("dataset","time"))
+x <- subset(x,variable=="reports")
+x <- transform(
+ x,
+ size=ifelse(grepl("small",dataset),"small popn","big popn"),
+ model=sub("(.+?)\\..*","\\1 model",dataset,perl=T)
+ )
+x <- rename(x,c(value="reports"))
+pl <- ggplot(data=x,mapping=aes(x=time,y=reports))+
+ geom_line()+facet_wrap(model~size,scale="free",nrow=4)
+print(pl)
-\begin{figure}[h]
- \begin{center}
-<<SEIRS-plots,fig=T,echo=F>>=
-ds <- c("SEIRS.small","SEIRS.big")
-<<dataplot>>
@
\end{center}
\caption{
- Plots of the simulated data, SEIRS model with initial conditions chosen to give stationary dynamics.
+ Plots of the simulated pertussis data.
Population size \scinot{5}{5} (black); \scinot{5}{6} (red).
Top: the full data sets.
Bottom: the first 20~yr of simulated data.
}
- \label{fig:seirs-plots}
+ \label{fig:dataplots}
\end{figure}
-\begin{figure}[h]
- \begin{center}
-<<SEIRR-plots,fig=T,echo=F>>=
-ds <- c("SEIRR.small","SEIRR.big")
-<<dataplot>>
-@
- \end{center}
- \caption{
- Plots of the simulated data, SEIRR model with initial conditions chosen to give stationary dynamics.
- Population size \scinot{5}{5} (black); \scinot{5}{6} (red).
- Top: the full data sets.
- Bottom: the first 20~yr of simulated data.
- }
- \label{fig:seirr-plots}
-\end{figure}
+We can get a benchmark for likelihood-based fitting methods by computing the true likelihood at the true parameter values.
+To do this, we run the \code{pfilter} particle filtering code.
+<<true-loglik-calc,echo=F,eval=F>>=
+require(Rmpi)
+require(mpifarm)
+require(pompExamples)
-\begin{figure}[h]
- \begin{center}
-<<full-plots,fig=T,echo=F>>=
-ds <- c("full.small","full.big")
-<<dataplot>>
-@
- \end{center}
- \caption{
- Plots of the simulated data, full model with initial conditions chosen to give stationary dynamics.
- Population size \scinot{5}{5} (black); \scinot{5}{6} (red).
- Top: the full data sets.
- Bottom: the first 20~yr of simulated data.
- }
- \label{fig:full-plots}
-\end{figure}
+data(pertussis.sim)
+ncpus <- as.integer(Sys.getenv("PBS_NP"))
+
+nrep <- 10 ### number of particle filters to run
+Np <- 10000 ### number of particles
+
+set.seed(5384959)
+
+mpi.spawn.Rslaves(nslaves=ncpus)
+
+mpi.farmer(
+ chunk=1,
+ seeds=as.integer(floor(runif(n=nrep,1,2^31))),
+ jobs={
+ require(plyr)
+ dlply(
+ expand.grid(
+ seed=seeds,
+ dataset=names(pertussis.sim)
+ ),
+ ~dataset+seed
+ )
+ },
+ common=list(
+ pomps=pertussis.sim,
+ Np=Np
+ ),
+ main={
+ require(pompExamples)
+ save.seed <- .Random.seed
+ set.seed(seed)
+ tic <- Sys.time()
+ pf <- try(pfilter(pomps[[dataset]],Np=Np,max.fail=100,warn=FALSE))
+ if (inherits(pf,"try-error")) {
+ loglik <- NA
+ nfail <- NA
+ } else {
+ loglik <- logLik(pf)
+ nfail <- pf$nfail
+ }
+ toc <- Sys.time()
+ .Random.seed <<- save.seed
+ data.frame(
+ seed=seed,
+ dataset=dataset,
+ loglik=loglik,
+ nfail=nfail,
+ etime=toc-tic
+ )
+ },
+ post={
+ require(plyr)
+ ldply(results)
+ }
+ ) -> results
+
+mpi.close.Rslaves()
+
+if (any(with(results,nfail!=0)))
+ warning("filtering failures occurred!")
+
+etime <- sum(results$etime)
+
+ll.est <- function (x) {
+ bl <- mean(x$loglik)
+ loglik <- bl+log(mean(exp(x$loglik-bl)))
+ se <- sd(exp(x$loglik-bl))/exp(loglik-bl)
+ data.frame(loglik=loglik,se=se)
+}
+
+loglik.truth <- ddply(results,~dataset,ll.est)
+
+save(loglik.truth,etime,Np,nrep,ncpus,file=binary.file,compress="xz")
+@
+<<true-loglik-eval,echo=F,results=hide>>=
+binary.file <- "pertussis-model-true-loglik.rda"
+if (file.exists(binary.file)) {
+ load(binary.file)
+} else {
+<<true-loglik-calc>>
+}
+@
+Table~\ref{tab:true-loglik} shows these likelihoods.
+
+<<true-loglik,results=tex,echo=F>>=
+print(
+ xtable(
+ loglik.truth,
+ caption=paste(
+ "Estimated log likelihood at the true parameters for the simulated pertussis data. ",
+ "To obtain these,",nrep,"particle filtering runs, each with",Np,"particles, were used. ",
+ "The column labeled",dQuote("se"),"gives the standard error of the Monte Carlo likelihood calculation. ",
+ "The computation took ",signif(etime,2),"~CPU~",units(etime),"on inexpensive processors. ",
+ sep=" "
+ ),
+ label="tab:true-loglik",
+ digits=c(0,0,1,2)
+ ),
+ type="latex",
+ floating=TRUE,
+ caption.placement="top",
+ include.rownames=FALSE,
+ sanitize.text.function=identity,
+ hline.after=c(-1,-1,0,nrow(loglik.truth))
+ )
+@
+
\newpage
\section{Design of the study}
Modified: pkg/pompExamples/inst/doc/pertussis-model.pdf
===================================================================
(Binary files differ)
Modified: pkg/pompExamples/tests/budmoth.Rout.save
===================================================================
--- pkg/pompExamples/tests/budmoth.Rout.save 2012-04-17 21:37:50 UTC (rev 666)
+++ pkg/pompExamples/tests/budmoth.Rout.save 2012-04-18 10:01:33 UTC (rev 667)
@@ -1,7 +1,8 @@
-R version 2.11.1 (2010-05-31)
-Copyright (C) 2010 The R Foundation for Statistical Computing
+R version 2.15.0 (2012-03-30)
+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.
@@ -15,51 +16,116 @@
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
-> require(pomp.devel)
-Loading required package: pomp.devel
+> require(pompExamples)
+Loading required package: pompExamples
Loading required package: pomp
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: deSolve
>
> data(budmoth.sim)
>
> names(budmoth.sim)
-[1] "tri.1" "para.1" "food.1" "para.2"
+[1] "tri" "para1" "food" "para2"
> x <- lapply(budmoth.sim,as,"data.frame")
>
> print(lapply(x,tail))
-$tri.1
- time Qobs Nobs Sobs
-996 995 30.91647 28.122473104 6.067451e-01
-997 996 25.61859 7.204479070 9.085794e-01
-998 997 27.65132 0.034050505 8.923623e-01
-999 998 32.67276 0.002703389 2.093088e-02
-1000 999 31.88108 0.006839399 5.231837e-05
-1001 1000 35.65827 0.078684854 1.684440e-06
+$tri
+ time Qobs Nobs Sobs Q N S
+56 55 34.90777 1.3909427 8.760997e-07 0.9860433 1.0105838 1.109173e-06
+57 56 33.36428 5.9679756 3.009644e-06 0.9830185 17.7373560 2.949154e-06
+58 57 29.34035 68.3377210 7.255485e-05 0.8608373 80.4984855 8.681423e-05
+59 58 21.75336 1.2686716 8.552474e-03 0.6310693 0.7176525 9.751675e-03
+60 59 29.15848 0.3500365 1.047452e-02 0.8025905 0.3723865 1.256469e-02
+61 60 32.98213 1.6128525 6.718390e-03 0.8941107 1.0912736 7.410424e-03
+ Alpha Lambda A
+56 0.5246694 21.95700 1.813852
+57 0.4833379 21.87970 1.738908
+58 0.5010435 22.22843 1.640583
+59 0.5190030 22.45659 1.404177
+60 0.5165432 22.20828 1.810056
+61 0.5183574 22.01444 1.591250
-$para.1
- time Qobs Nobs Sobs
-996 995 24.48653 9.6502436 0.8652499
-997 996 26.86343 1.7416087 0.8813915
-998 997 30.67217 0.3623039 0.7700699
-999 998 32.92757 1.7122017 0.4362900
-1000 999 34.28261 16.2801567 0.6302552
-1001 1000 30.61003 39.0414370 0.8830192
+$para1
+ time Qobs Nobs Sobs Q N S
+56 55 21.38207 0.68381238 0.907924733 0.6143380 0.55733084 0.997889584
+57 56 29.54776 0.02697955 0.430252439 0.8126219 0.02020849 0.573040310
+58 57 31.78574 0.07253569 0.020786180 0.9177534 0.17066444 0.022035629
+59 58 32.39154 2.81430772 0.005905631 0.9567875 3.49854944 0.005907819
+60 59 34.36420 63.29278180 0.031631312 0.9456880 56.35560023 0.033871059
+61 60 24.63131 17.42221294 0.790037174 0.7153921 12.65588781 0.871324054
+ Alpha Lambda A
+56 0.5031487 22.14524 1.579994
+57 0.4707325 21.78521 1.754212
+58 0.4377213 21.75771 1.930505
+59 0.5049253 22.14170 1.576733
+60 0.4996302 21.90574 1.675765
+61 0.5157446 21.68244 1.551352
-$food.1
- time Qobs Nobs Sobs
-996 995 28.70032 0.0024043355 4.457058e-01
-997 996 30.32688 0.0008529304 9.675072e-04
-998 997 33.63074 0.0015256083 1.852273e-06
-999 998 34.77058 0.0046225425 7.992639e-07
-1000 999 34.80885 0.0163346510 7.917979e-07
-1001 1000 34.54063 0.0927045595 8.916789e-07
+$food
+ time Qobs Nobs Sobs Q N S
+56 55 34.13985 16.01354997 3.224917e-06 0.9529645 7.4387892 3.827248e-06
+57 56 27.79192 22.28079514 2.070576e-05 0.8502012 19.7570009 2.570550e-05
+58 57 24.07011 31.94991132 4.869824e-04 0.6794912 14.4355348 5.305078e-04
+59 58 21.81314 1.77819709 7.707581e-03 0.6288764 1.9239332 7.967758e-03
+60 59 27.66444 0.24999819 1.516026e-02 0.7899212 0.1973302 1.416005e-02
+61 60 28.79342 0.09924733 2.128647e-03 0.8898312 0.1200628 3.272257e-03
+ Alpha Lambda A
+56 0.5153850 4.617116 1.0869780
+57 0.5413771 4.932935 0.8677828
+58 0.5082460 4.851747 1.0428964
+59 0.4871643 4.388592 1.0444642
+60 0.4316777 4.394991 0.9302542
+61 0.5012195 5.063878 1.1726507
-$para.2
- time Qobs Nobs Sobs
-996 995 34.51711 15.30719787 0.0004411851
-997 996 29.61491 8.38435616 0.0030955092
-998 997 29.11093 20.96921368 0.1184931629
-999 998 26.50715 32.89609224 0.8868716606
-1000 999 23.49514 0.46694271 0.9023910818
-1001 1000 27.95434 0.02132294 0.7515883032
+$para2
+ time Qobs Nobs Sobs Q N S
+56 55 33.68890 3.458362 1.054469e-06 0.9845993 1.799560 1.225629e-06
+57 56 33.87640 15.747819 2.248998e-06 0.9749960 6.921067 2.575005e-06
+58 57 33.94564 50.157117 1.166632e-05 0.9265344 33.878434 1.350848e-05
+59 58 26.54568 39.257279 7.764646e-04 0.7668865 21.145329 1.028206e-03
+60 59 25.06245 37.080107 1.942558e-02 0.7334130 32.430161 2.124212e-02
+61 60 24.65754 10.630755 8.364070e-01 0.6707489 14.812001 9.357010e-01
+ Alpha Lambda A
+56 0.5199566 23.368648 2.8482104
+57 0.5034550 4.475840 0.7140993
+58 0.4982577 8.622703 0.7018725
+59 0.5168533 9.734355 2.2460464
+60 0.4777904 9.364085 0.9906918
+61 0.5059549 7.138668 6.7705675
>
+> y <- simulate(budmoth.sim$food,seed=3434996L,as.data.frame=TRUE)
+> tail(y)
+ time Qobs Nobs Sobs Q N S
+56 55 21.60813 9.9416747 0.0084565260 0.6260519 6.6374470 0.0094829150
+57 56 23.69758 0.4068012 0.0555572562 0.6960924 0.6909030 0.0655207892
+58 57 27.36280 0.1056701 0.0461785098 0.8328430 0.1516087 0.0367590933
+59 58 31.64283 0.2311652 0.0054429206 0.9128312 0.1262158 0.0061363639
+60 59 33.07830 0.2439202 0.0007323085 0.9501525 0.2684836 0.0007689092
+61 60 34.70798 0.4126931 0.0001924636 0.9684718 0.8167865 0.0002084246
+ Alpha Lambda A sim
+56 0.4538116 5.723480 1.0150383 1
+57 0.4386467 5.049427 1.0766197 1
+58 0.4944825 4.972702 0.8273030 1
+59 0.4989262 4.612384 1.1043024 1
+60 0.5386598 5.130300 0.9918653 1
+61 0.4994891 5.038971 1.0048781 1
+>
+> z <- trajectory(budmoth.sim$tri,as.data.frame=TRUE)
+> tail(z)
+ Q N S Alpha Lambda A time traj
+56 1 0 9.99998e-07 0.5 22 1.7 55 1
+57 1 0 9.99998e-07 0.5 22 1.7 56 1
+58 1 0 9.99998e-07 0.5 22 1.7 57 1
+59 1 0 9.99998e-07 0.5 22 1.7 58 1
+60 1 0 9.99998e-07 0.5 22 1.7 59 1
+61 1 0 9.99998e-07 0.5 22 1.7 60 1
+>
+> pf <- pfilter(budmoth.sim$food,seed=34348885L,Np=1000)
+> logLik(pf)
+[1] 355.9053
+>
+> proc.time()
+ user system elapsed
+ 0.564 0.020 0.602
Modified: pkg/pompExamples/tests/pertussis.R
===================================================================
--- pkg/pompExamples/tests/pertussis.R 2012-04-17 21:37:50 UTC (rev 666)
+++ pkg/pompExamples/tests/pertussis.R 2012-04-18 10:01:33 UTC (rev 667)
@@ -13,8 +13,5 @@
y <- trajectory(pertussis.sim$SEIRS.small,as.data.frame=TRUE)
tail(y)
-x <- simulate(pertussis.sim$full.small,seed=395885L,times=100,states=TRUE)[1:6,,1]
-coef(pertussis.sim$full.small,paste0(names(x),".0")) <- x/sum(x)
-
pf <- pfilter(pertussis.sim$full.small,seed=3445886L,Np=1000)
logLik(pf)
Modified: pkg/pompExamples/tests/pertussis.Rout.save
===================================================================
--- pkg/pompExamples/tests/pertussis.Rout.save 2012-04-17 21:37:50 UTC (rev 666)
+++ pkg/pompExamples/tests/pertussis.Rout.save 2012-04-18 10:01:33 UTC (rev 667)
@@ -1,7 +1,8 @@
-R version 2.11.1 (2010-05-31)
-Copyright (C) 2010 The R Foundation for Statistical Computing
+R version 2.15.0 (2012-03-30)
+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.
@@ -15,88 +16,139 @@
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
-> require(pomp.devel)
-Loading required package: pomp.devel
+> require(pompExamples)
+Loading required package: pompExamples
Loading required package: pomp
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: deSolve
>
> data(pertussis.sim)
>
> names(pertussis.sim)
[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(pertussis.sim,as.data.frame)
>
> print(lapply(x,tail))
$SEIR.small
- time reports
-5196 99.90385 57
-5197 99.92308 71
-5198 99.94231 74
-5199 99.96154 60
-5200 99.98077 67
-5201 100.00000 51
+ time reports S E I R1 R2 V cases W err simpop
+1036 19.90385 70 26873 244 409 472429 0 0 225 0 7 499955
+1037 19.92308 69 26810 239 435 472458 0 0 221 0 7 499942
+1038 19.94231 54 26797 244 435 472497 0 0 208 0 7 499973
+1039 19.96154 46 26746 274 417 472542 0 0 205 0 7 499979
+1040 19.98077 53 26709 256 450 472538 0 0 228 0 7 499953
+1041 20.00000 71 26679 273 415 472560 0 0 183 0 7 499927
$SEIR.big
- time reports
-5196 99.90385 435
-5197 99.92308 460
-5198 99.94231 376
-5199 99.96154 478
-5200 99.98077 428
-5201 100.00000 463
+ time reports S E I R1 R2 V cases W err simpop
+1036 19.90385 520 255547 1703 2917 4739686 0 0 1470 0 7 4999853
+1037 19.92308 484 255784 1759 3003 4739211 0 0 1489 0 7 4999757
+1038 19.94231 406 256081 1856 2995 4738900 0 0 1444 0 7 4999832
+1039 19.96154 533 256320 1869 3129 4738606 0 0 1597 0 7 4999924
+1040 19.98077 425 256551 1879 3164 4738317 0 0 1545 0 7 4999911
+1041 20.00000 412 257106 1639 3151 4738010 0 0 1495 0 7 4999906
$SEIRS.small
- time reports
-5196 99.90385 48
-5197 99.92308 49
-5198 99.94231 56
-5199 99.96154 59
-5200 99.98077 68
-5201 100.00000 51
+ time reports S E I R1 R2 V cases W err simpop
+1036 19.90385 47 80392 539 920 228000 190629 0 467 0 7 500480
+1037 19.92308 41 80392 578 909 227952 190651 0 431 0 7 500482
+1038 19.94231 49 80395 591 914 227879 190690 0 464 0 7 500469
+1039 19.96154 67 80471 575 975 227795 190660 0 500 0 7 500476
+1040 19.98077 42 80516 553 993 227732 190669 0 513 0 7 500463
+1041 20.00000 61 80660 492 986 227653 190677 0 461 0 7 500468
$SEIRS.big
- time reports
-5196 99.90385 425
-5197 99.92308 416
-5198 99.94231 507
-5199 99.96154 514
-5200 99.98077 568
-5201 100.00000 466
+ time reports S E I R1 R2 V cases W err simpop
+1036 19.90385 489 773965 5409 9193 2299336 1912393 0 4484 0 7 5000296
+1037 19.92308 568 774652 5472 9318 2298462 1912401 0 4574 0 7 5000305
+1038 19.94231 397 775170 5539 9426 2297714 1912390 0 4671 0 7 5000239
+1039 19.96154 411 775614 5740 9491 2297059 1912330 0 4609 0 7 5000234
+1040 19.98077 442 776437 5571 9697 2296283 1912273 0 4748 0 7 5000261
+1041 20.00000 518 778090 4858 9441 2295594 1912196 0 4335 0 7 5000179
$SEIRR.small
- time reports
-5196 99.90385 24
-5197 99.92308 34
-5198 99.94231 37
-5199 99.96154 50
-5200 99.98077 26
-5201 100.00000 30
+ time reports S E I R1 R2 V cases W err simpop
+1036 19.90385 73 64537 864 1441 317687 115760 0 714 0 7 500289
+1037 19.92308 69 64260 896 1502 318088 115525 0 758 0 7 500271
+1038 19.94231 76 63954 981 1543 318494 115282 0 748 0 7 500254
+1039 19.96154 84 63664 965 1638 318953 115037 0 836 0 7 500257
+1040 19.98077 91 63373 992 1670 319485 114723 0 802 0 7 500243
+1041 20.00000 81 63266 861 1646 319863 114604 0 760 0 7 500240
$SEIRR.big
- time reports
-5196 99.90385 449
-5197 99.92308 549
-5198 99.94231 608
-5199 99.96154 503
-5200 99.98077 525
-5201 100.00000 432
+ time reports S E I R1 R2 V cases W err simpop
+1036 19.90385 449 636820 4839 8046 3202063 1148211 0 3956 0 7 4999979
+1037 19.92308 472 637284 4951 8179 3201137 1148414 0 4114 0 7 4999965
+1038 19.94231 438 637653 5019 8240 3200093 1148879 0 4143 0 7 4999884
+1039 19.96154 474 637980 5084 8458 3199122 1149316 0 4261 0 7 4999960
+1040 19.98077 499 638422 4929 8626 3198103 1149922 0 4307 0 7 5000002
+1041 20.00000 469 639527 4413 8258 3196671 1151115 0 3834 0 7 4999984
$full.small
- time reports
-5196 99.90385 44
-5197 99.92308 28
-5198 99.94231 42
-5199 99.96154 42
-5200 99.98077 61
-5201 100.00000 44
+ time reports S E I R1 R2 V cases W err simpop
+1036 19.90385 38 60812 492 892 330814 107242 0 437 -9.164189 7 500252
+1037 19.92308 35 60904 514 881 330659 107326 0 413 -9.153313 7 500284
+1038 19.94231 44 60988 464 873 330490 107453 0 420 -9.279501 7 500268
+1039 19.96154 46 61032 517 876 330295 107543 0 398 -9.092366 7 500263
+1040 19.98077 39 61122 504 845 330138 107631 0 403 -9.103926 7 500240
+1041 20.00000 43 61258 474 820 329898 107806 0 384 -8.973250 7 500256
$full.big
- time reports
-5196 99.90385 714
-5197 99.92308 712
-5198 99.94231 741
-5199 99.96154 787
-5200 99.98077 799
-5201 100.00000 550
+ time reports S E I R1 R2 V cases W err
+1036 19.90385 381 685681 5306 8066 3166871 1137704 0 3943 -5.397648 7
+1037 19.92308 526 685987 5712 8423 3165789 1137727 0 4434 -5.215067 7
+1038 19.94231 406 686713 5594 8728 3164496 1138137 0 4487 -5.376581 7
+1039 19.96154 377 687178 5673 8981 3163446 1138391 0 4463 -5.484270 7
+1040 19.98077 549 687709 5515 9310 3162464 1138650 0 4717 -5.563903 7
+1041 20.00000 419 688789 4914 9220 3161008 1139607 0 4358 -5.568812 7
+ simpop
+1036 5003628
+1037 5003638
+1038 5003668
+1039 5003669
+1040 5003648
+1041 5003538
>
+> x <- simulate(pertussis.sim$full.big,seed=395885L,as.data.frame=TRUE)
+> tail(x)
+ time reports S E I R1 R2 V cases W err
+1036 19.90385 529 634087 5541 9277 3260155 1090900 0 4646 0.9650237 7
+1037 19.92308 345 634343 5687 9453 3259271 1091140 0 4625 0.9681362 7
+1038 19.94231 459 634432 6021 9647 3258692 1091145 0 4844 1.0560729 7
+1039 19.96154 391 634660 6040 9905 3258019 1091364 0 4918 1.0452093 7
+1040 19.98077 565 635319 5595 10107 3257412 1091594 0 5033 0.8786933 7
+1041 20.00000 411 636665 4728 9753 3256002 1092869 0 4394 0.7840614 7
+ simpop sim
+1036 4999960 1
+1037 4999894 1
+1038 4999937 1
+1039 4999988 1
+1040 5000027 1
+1041 5000017 1
+>
+> y <- trajectory(pertussis.sim$SEIRS.small,as.data.frame=TRUE)
+> tail(y)
+ S E I R1 R2 V cases W err simpop
+1036 81409.73 558.4599 942.3100 227353.0 189736.5 0 0 0 0 5e+05
+1037 81420.36 573.0107 965.5155 227305.3 189735.8 0 0 0 0 5e+05
+1038 81418.14 587.5604 989.6877 227269.6 189735.0 0 0 0 0 5e+05
+1039 81402.73 602.3328 1014.5493 227246.2 189734.2 0 0 0 0 5e+05
+1040 81415.22 580.7639 1035.4780 227235.2 189733.4 0 0 0 0 5e+05
+1041 81532.40 510.2405 1002.6699 227222.2 189732.5 0 0 0 0 5e+05
+ time traj
+1036 19.90385 1
+1037 19.92308 1
+1038 19.94231 1
+1039 19.96154 1
+1040 19.98077 1
+1041 20.00000 1
+>
+> pf <- pfilter(pertussis.sim$full.small,seed=3445886L,Np=1000)
+> logLik(pf)
+[1] -3829.33
+>
+> proc.time()
+ user system elapsed
+ 22.637 0.056 22.773
More information about the pomp-commits
mailing list