[Pomp-commits] r1019 - in pkg/pompExamples: . R inst inst/examples src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 17 20:12:01 CET 2014
Author: kingaa
Date: 2014-12-17 20:12:00 +0100 (Wed, 17 Dec 2014)
New Revision: 1019
Added:
pkg/pompExamples/R/aaa.R
pkg/pompExamples/inst/examples/
pkg/pompExamples/inst/examples/bbp.R
pkg/pompExamples/inst/examples/parus.R
pkg/pompExamples/src/parus.c
pkg/pompExamples/tests/examples.R
Modified:
pkg/pompExamples/DESCRIPTION
pkg/pompExamples/inst/NEWS
pkg/pompExamples/inst/NEWS.Rd
pkg/pompExamples/tests/pertussis.Rout.save
Log:
- revamp the 'pompExample' facility to accept a path of directories in which to search
- it is now possible to pass variables to the 'pompExample' scripts
- move the 'parus' example to the pompExamples package
- new "bbp" example for the Bombay plague outbreak of 1905-06.
- pompExample can now return the objects as a list instead of creating them in a specified environment, if desired
Modified: pkg/pompExamples/DESCRIPTION
===================================================================
--- pkg/pompExamples/DESCRIPTION 2014-12-17 15:45:48 UTC (rev 1018)
+++ pkg/pompExamples/DESCRIPTION 2014-12-17 19:12:00 UTC (rev 1019)
@@ -1,8 +1,8 @@
Package: pompExamples
Type: Package
Title: Additional pomp examples
-Version: 0.23-3
-Date: 2014-07-15
+Version: 0.24-1
+Date: 2014-12-16
Maintainer: Aaron A. King <kingaa at umich.edu>
Authors at R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),
email="kingaa at umich.edu"),
@@ -21,4 +21,4 @@
License: GPL (>= 2)
LazyData: false
BuildVignettes: true
-Collate: budmoth.R pertussis.R
+Collate: aaa.R budmoth.R pertussis.R
Added: pkg/pompExamples/R/aaa.R
===================================================================
--- pkg/pompExamples/R/aaa.R (rev 0)
+++ pkg/pompExamples/R/aaa.R 2014-12-17 19:12:00 UTC (rev 1019)
@@ -0,0 +1,12 @@
+.onAttach <- function (...) {
+ exampleDir <- getOption("pomp.examples")
+ newDir <- system.file("examples",package="pompExamples")
+ options(pomp.examples=c(exampleDir,newDir,recursive=TRUE))
+}
+
+.onDetach <- function (...) {
+ exampleDir <- getOption("pomp.examples")
+ newDir <- system.file("examples",package="pompExamples")
+ exampleDir <- exampleDir[exampleDir!=newDir]
+ options(pomp.examples=exampleDir)
+}
Modified: pkg/pompExamples/inst/NEWS
===================================================================
--- pkg/pompExamples/inst/NEWS 2014-12-17 15:45:48 UTC (rev 1018)
+++ pkg/pompExamples/inst/NEWS 2014-12-17 19:12:00 UTC (rev 1019)
@@ -1,5 +1,12 @@
_N_e_w_s _f_o_r _P_a_c_k_a_g_e '_p_o_m_p_E_x_a_m_p_l_e_s'
+_C_h_a_n_g_e_s _i_n '_p_o_m_p_E_x_a_m_p_l_e_s' _v_e_r_s_i_o_n _0._2_4-_1:
+
+ • Examples included in the package can now be accessed using
+ ‘pomp’'s ‘pompExample’ function.
+
+ • New _Parus major_ example.
+
_C_h_a_n_g_e_s _i_n '_p_o_m_p_E_x_a_m_p_l_e_s' _v_e_r_s_i_o_n _0._2_3-_2:
• Update for use with ‘pomp’ version 0.50.
Modified: pkg/pompExamples/inst/NEWS.Rd
===================================================================
--- pkg/pompExamples/inst/NEWS.Rd 2014-12-17 15:45:48 UTC (rev 1018)
+++ pkg/pompExamples/inst/NEWS.Rd 2014-12-17 19:12:00 UTC (rev 1019)
@@ -1,5 +1,11 @@
\name{NEWS}
\title{News for Package 'pompExamples'}
+\section{Changes in \pkg{pompExamples} version 0.24-1}{
+ \itemize{
+ \item Examples included in the package can now be accessed using \pkg{pomp}'s \code{pompExample} function.
+ \item New \emph{Parus major} example.
+ }
+}
\section{Changes in \pkg{pompExamples} version 0.23-2}{
\itemize{
\item Update for use with \pkg{pomp} version 0.50.
Added: pkg/pompExamples/inst/examples/bbp.R
===================================================================
--- pkg/pompExamples/inst/examples/bbp.R (rev 0)
+++ pkg/pompExamples/inst/examples/bbp.R 2014-12-17 19:12:00 UTC (rev 1019)
@@ -0,0 +1,89 @@
+require(pomp)
+
+dat <- read.csv(text='
+## Deaths due to plague during an outbreak on the Island of Bombay
+## over the period 17 Dec 1905 to 21 July 1906.
+## From Kermack, W. O. & McKendrick, A. G. (1927)
+## A Contribution to the Mathematical Theory of Epidemics
+## Proceedings of the Royal Society of London, Series A, 115: 700--721.
+## "date" is date of end of the week (Saturday)
+"week","date","deaths"
+1,1905-12-23,4
+2,1905-12-30,10
+3,1906-01-06,15
+4,1906-01-13,18
+5,1906-01-20,21
+6,1906-01-27,31
+7,1906-02-03,51
+8,1906-02-10,53
+9,1906-02-17,97
+10,1906-02-24,125
+11,1906-03-03,183
+12,1906-03-10,292
+13,1906-03-17,390
+14,1906-03-24,448
+15,1906-03-31,641
+16,1906-04-07,771
+17,1906-04-14,701
+18,1906-04-21,696
+19,1906-04-28,867
+20,1906-05-05,925
+21,1906-05-12,801
+22,1906-05-19,580
+23,1906-05-26,409
+24,1906-06-02,351
+25,1906-06-09,210
+26,1906-06-16,113
+27,1906-06-23,65
+28,1906-06-30,52
+29,1906-07-07,51
+30,1906-07-14,39
+31,1906-07-21,33
+',comment.char="#")
+
+pomp(data=subset(dat,select=c(week,deaths)),
+ times="week",
+ t0=0,
+ params=c(
+ beta=2,delta=1.5,y0=0.0004,theta=54,
+ sigma=0.02,
+ mu=0,gamma=0.2,ratio=10000
+ ),
+ rprocess=euler.sim(
+ step.fun=Csnippet("
+ double X = exp(x);
+ double Y = exp(y);
+ double dx, dy, dn, dW, ito;
+ dx = (mu*(1.0/X-1)+(delta-beta)*Y)*dt;
+ dy = (beta*X+delta*(Y-1)-gamma-mu)*dt;
+ dn = -delta*Y*dt;
+ dW = rnorm(0,sigma*sqrt(dt));
+ ito = 0.5*sigma*sigma*dt;
+ x += dx - beta*Y*(dW-beta*Y*ito);
+ y += dy + beta*X*(dW+beta*X*ito);
+ n += dn;
+ "
+ ),
+ delta.t=1/24),
+ paramnames=c("beta","delta","mu","gamma","sigma","theta","ratio"),
+ statenames=c("x","y","n"),
+ measurement.model=deaths~nbinom(mu=ratio*exp(y),size=theta),
+ logvar=c("beta","delta","ratio","sigma","theta"),
+ logitvar=c("y0"),
+ parameter.inv.transform=function (params, logvar, logitvar, ...) {
+ params[logvar] <- log(params[logvar])
+ params[logitvar] <- qlogis(params[logitvar])
+ params
+ },
+ parameter.transform=function (params, logvar, logitvar, ...) {
+ params[logvar] <- exp(params[logvar])
+ params[logitvar] <- plogis(params[logitvar])
+ params
+ },
+ initializer=function(params, t0, ...) {
+ y0 <- unname(params["y0"])
+ c(x=log(1-y0),y=log(y0),n=log(1))
+ }
+ ) -> bbp
+
+c("bbp")
Added: pkg/pompExamples/inst/examples/parus.R
===================================================================
--- pkg/pompExamples/inst/examples/parus.R (rev 0)
+++ pkg/pompExamples/inst/examples/parus.R 2014-12-17 19:12:00 UTC (rev 1019)
@@ -0,0 +1,90 @@
+require(pomp)
+
+proc.avail <- c("Gompertz","Ricker")
+meas.avail <- c("lognormal","Poisson","negbin")
+
+if (!exists("proc",where=environment()))
+ stop("choose a process model: proc = ",sQuote(proc.avail))
+if (!exists("meas",where=environment()))
+ stop("choose a measurement model: meas = ",sQuote(meas.avail))
+
+proc <- proc.avail[pmatch(proc,proc.avail)]
+meas <- meas.avail[pmatch(meas,meas.avail)]
+
+dat <- 'year,pop
+1960,148
+1961,258
+1962,185
+1963,170
+1964,267
+1965,239
+1966,196
+1967,132
+1968,167
+1969,186
+1970,128
+1971,227
+1972,174
+1973,177
+1974,137
+1975,172
+1976,119
+1977,226
+1978,166
+1979,161
+1980,199
+1981,306
+1982,206
+1983,350
+1984,214
+1985,175
+1986,211
+'
+
+dat <- read.csv(text=dat)
+
+pomp(
+ data=dat,
+ times="year",
+ t0=1960,
+ params=c(K=190,r=2.7,sigma=0.2,theta=0.05,N.0=148),
+ rprocess=discrete.time.sim(
+ step.fun=switch(proc,
+ Gompertz="_parus_gompertz_simulator",
+ Ricker="_parus_ricker_simulator",
+ stop("unrecognized value of ",sQuote("proc"))
+ ),
+ delta.t=1
+ ),
+ skeleton=switch(proc,
+ Gompertz="_parus_gompertz_skeleton",
+ Ricker="_parus_ricker_skeleton",
+ stop("unrecognized value of ",sQuote("proc"))
+ ),
+ skeleton.type="map",
+ skelmap.delta.t=1,
+ rmeasure=switch(meas,
+ lognormal="_parus_lognormal_rmeasure",
+ Poisson="_parus_poisson_rmeasure",
+ negbin="_parus_nbinom_rmeasure",
+ stop("unrecognized value of ",sQuote("meas"))
+ ),
+ dmeasure=switch(meas,
+ lognormal="_parus_lognormal_dmeasure",
+ Poisson="_parus_poisson_dmeasure",
+ negbin="_parus_nbinom_dmeasure",
+ stop("unrecognized value of ",sQuote("meas"))
+ ),
+ paramnames=c("r","K","sigma","theta"),
+ statenames=c("N"),
+ obsnames=c("pop"),
+ parameter.transform=function(params,...){
+ exp(params)
+ },
+ parameter.inv.transform=function(params,...){
+ log(params)
+ },
+ PACKAGE="pompExamples"
+ ) -> parus
+
+c("parus")
Added: pkg/pompExamples/src/parus.c
===================================================================
--- pkg/pompExamples/src/parus.c (rev 0)
+++ pkg/pompExamples/src/parus.c 2014-12-17 19:12:00 UTC (rev 1019)
@@ -0,0 +1,87 @@
+// dear emacs, please treat this as -*- C++ -*-
+
+#include <Rmath.h>
+
+#include "pomp.h"
+
+#define R (p[parindex[0]]) // growth rate
+#define K (p[parindex[1]]) // carrying capacity
+#define SIGMA (p[parindex[2]]) // process noise level
+#define THETA (p[parindex[3]]) // measurement noise level
+
+#define POP (y[obsindex[0]])
+#define N (x[stateindex[0]])
+#define NPRIME (f[stateindex[0]])
+
+void _parus_lognormal_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ *lik = dlnorm(POP,log(N),THETA,give_log);
+}
+
+void _parus_lognormal_rmeasure (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ POP = rlnorm(log(N),THETA);
+}
+
+void _parus_poisson_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ *lik = dpois(POP,N+1.0e-10,give_log);
+}
+
+void _parus_poisson_rmeasure (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ POP = rpois(N+1.0e-10);
+}
+
+void _parus_nbinom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ *lik = dnbinom_mu(POP,1.0/THETA,N+1.0e-10,give_log);
+}
+
+void _parus_nbinom_rmeasure (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ POP = rnbinom_mu(1.0/THETA,N+1.0e-10);
+}
+
+void _parus_gompertz_simulator (double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar,
+ double t, double dt)
+{
+ double S = exp(-R*dt);
+ double eps = (SIGMA > 0.0) ? exp(rnorm(0,SIGMA)) : 1.0;
+ N = pow(K,(1-S))*pow(N,S)*eps;
+}
+
+// the deterministic skeleton
+void _parus_gompertz_skeleton (double *f, double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar, double t)
+{
+ double dt = 1.0;
+ double S = exp(-R*dt);
+ NPRIME = pow(K,(1-S))*pow(N,S);
+}
+
+// Ricker model with log-normal process noise
+void _parus_ricker_simulator (double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar,
+ double t, double dt)
+{
+ double e = (SIGMA > 0.0) ? rnorm(0,SIGMA) : 0.0;
+ N = exp(log(N)+R*(1-N/K)+e);
+}
+
+void _parus_ricker_skeleton (double *f, double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar, double t)
+{
+ NPRIME = exp(log(N)+R*(1-N/K));
+}
Added: pkg/pompExamples/tests/examples.R
===================================================================
--- pkg/pompExamples/tests/examples.R (rev 0)
+++ pkg/pompExamples/tests/examples.R 2014-12-17 19:12:00 UTC (rev 1019)
@@ -0,0 +1,28 @@
+library(pompExamples)
+
+## pdf.options(useDingbats=FALSE)
+## pdf(file="examples.pdf")
+
+set.seed(47575684L)
+
+po <- pompExample(parus,proc="Ricker",meas="lognormal",envir=NULL)
+pf <- pfilter(simulate(po$parus),Np=100,max.fail=Inf)
+tj <- trajectory(po$parus)
+
+po <- pompExample(parus,proc="Ricker",meas="negbin",envir=NULL)
+pf <- pfilter(simulate(po$parus),Np=100,max.fail=Inf)
+
+po <- pompExample(parus,proc="Ricker",meas="Poisson",envir=NULL)
+pf <- pfilter(simulate(po$parus),Np=100,max.fail=Inf)
+
+po <- pompExample(parus,proc="Gompertz",meas="Poisson",envir=NULL)
+pf <- pfilter(simulate(po[[1]]),Np=100,max.fail=Inf)
+tj <- trajectory(po[[1]])
+
+po <- pompExample(parus,proc="Gompertz",meas="lognormal",envir=NULL)
+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)
+
+## dev.off()
Modified: pkg/pompExamples/tests/pertussis.Rout.save
===================================================================
--- pkg/pompExamples/tests/pertussis.Rout.save 2014-12-17 15:45:48 UTC (rev 1018)
+++ pkg/pompExamples/tests/pertussis.Rout.save 2014-12-17 19:12:00 UTC (rev 1019)
@@ -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,9 +17,8 @@
> library(pompExamples)
Loading required package: pomp
-Loading required package: mvtnorm
Loading required package: subplex
-Loading required package: deSolve
+Loading required package: nloptr
>
> all <- c("SEIR.small","SEIR.big","SEIRS.small","SEIRS.big","SEIRR.small","SEIRR.big","full.small","full.big")
>
@@ -147,7 +146,7 @@
>
> system.time(pf <- pfilter(pertussis.sim(full.small),seed=3445886L,Np=1000))
user system elapsed
- 19.717 0.000 19.778
+ 18.793 0.000 18.864
> logLik(pf)
[1] -3829.33
>
@@ -171,4 +170,4 @@
>
> proc.time()
user system elapsed
- 20.497 0.064 20.642
+ 19.533 0.040 19.671
More information about the pomp-commits
mailing list