[Pomp-commits] r703 - in pkg/pomp: . R demo
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 3 18:27:45 CEST 2012
Author: kingaa
Date: 2012-05-03 18:27:45 +0200 (Thu, 03 May 2012)
New Revision: 703
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/R/builder.R
pkg/pomp/R/trajectory-pomp.R
pkg/pomp/demo/sir.R
Log:
- fix bug in 'trajectory' for vectorfield case that shows up when 'zeronames' has length > 1
- re-work 'sir' demo so that it uses 'pompBuilder'
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2012-05-02 16:50:38 UTC (rev 702)
+++ pkg/pomp/DESCRIPTION 2012-05-03 16:27:45 UTC (rev 703)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.42-1
-Date: 2012-05-02
+Version: 0.42-2
+Date: 2012-05-03
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>
URL: http://pomp.r-forge.r-project.org
Modified: pkg/pomp/R/builder.R
===================================================================
--- pkg/pomp/R/builder.R 2012-05-02 16:50:38 UTC (rev 702)
+++ pkg/pomp/R/builder.R 2012-05-03 16:27:45 UTC (rev 703)
@@ -99,8 +99,12 @@
skeleton="\n}\n\n"
)
-pompCBuilder <- function (name, statenames, paramnames, obsnames, rmeasure, dmeasure, step.fn, skeleton) {
-
+pompCBuilder <- function (name, statenames, paramnames, obsnames, rmeasure, dmeasure, step.fn, skeleton)
+{
+ if (missing(name)) stop(sQuote("name")," must be supplied");
+ if (missing(statenames)) stop(sQuote("name")," must be supplied");
+ if (missing(paramnames)) stop(sQuote("name")," must be supplied");
+ if (missing(obsnames)) stop(sQuote("name")," must be supplied");
name <- cleanForC(name)
statenames <- cleanForC(statenames)
paramnames <- cleanForC(paramnames)
Modified: pkg/pomp/R/trajectory-pomp.R
===================================================================
--- pkg/pomp/R/trajectory-pomp.R 2012-05-02 16:50:38 UTC (rev 702)
+++ pkg/pomp/R/trajectory-pomp.R 2012-05-03 16:27:45 UTC (rev 703)
@@ -92,8 +92,8 @@
}
- if (length(znames)>0)
- x[znames,,-1] <- apply(x[znames,,,drop=FALSE],c(1,2),diff)
+ for (z in znames)
+ x[z,,-1] <- apply(x[z,,,drop=FALSE],c(1,2),diff)
if (as.data.frame) {
x <- lapply(
Modified: pkg/pomp/demo/sir.R
===================================================================
--- pkg/pomp/demo/sir.R 2012-05-02 16:50:38 UTC (rev 702)
+++ pkg/pomp/demo/sir.R 2012-05-03 16:27:45 UTC (rev 703)
@@ -1,112 +1,165 @@
require(pomp)
-## Coding up the SIR example using native routines results in much faster computations.
-## The C codes "sir_euler_simulator", "sir_ODE", "binom_rmeasure", and "binom_dmeasure"
-## are included in the "examples" directory (file "sir.c")
+dmeas <- "
+ lik = dbinom(reports,nearbyint(cases),rho,give_log);
+"
+rmeas <- "
+ reports = rbinom(nearbyint(cases),rho);
+"
+step.fn <- '
+ int nrate = 6;
+ double rate[nrate]; // transition rates
+ double trans[nrate]; // transition numbers
+ double beta; // transmission rate
+ double dW; // white noise increment
+ double period = 1.0; // period of the seasonality
+ int nbasis = 3; // number of seasonality basis functions
+ int deg = 3; // degree of the B-spline basis functions
+ double seasonality[nbasis];
+ int k;
-if (Sys.info()['sysname']=='Linux') { # only run this under linux
+ // compute transmission rate from seasonality
+ periodic_bspline_basis_eval(t,period,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
+ for (k = 0, beta = 0; k < nbasis; k++)
+ beta += (&beta1)[k]*seasonality[k];
- model <- "sir"
- ## name of the file holding the model codes:
- modelfile <- system.file("examples",paste(model,".c",sep=""),package="pomp")
- ## name of the shared-object library:
- solib <- paste(model,.Platform$dynlib.ext,sep="")
- ## environment variables needed to locate the pomp header file:
- cflags <- paste("PKG_CFLAGS=\"-I",system.file("include/",package="pomp"),"\"")
+ // compute the environmental stochasticity
+ dW = rgammawn(beta_sd,dt);
- ## compile the shared-object library containing the model codes:
- rv <- system2(
- command=R.home("bin/R"),
- args=c("CMD","SHLIB","-o",solib,modelfile),
- env=cflags
- )
- if (rv!=0)
- stop("cannot compile shared-object library ",sQuote(solib))
+ // compute the transition rates
+ rate[0] = mu*popsize; // birth into susceptible class
+ rate[1] = (iota+beta*I*dW/dt)/popsize; // force of infection
+ rate[2] = mu; // death from susceptible class
+ rate[3] = gamma; // recovery
+ rate[4] = mu; // death from infectious class
+ rate[5] = mu; // death from recovered class
- 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", # native routine for the simulation step
- delta.t=1/52/20 # Euler stepsize: 1/20 wk
- ),
- skeleton.type="vectorfield",
- skeleton="sir_ODE", # native routine for the skeleton
- rmeasure="binomial_rmeasure", # binomial measurement model
- dmeasure="binomial_dmeasure", # binomial measurement model
- PACKAGE="sir", ## name of the shared-object library
- ## the order of the observables assumed in the native routines:
- obsnames="reports",
- ## the order of the state variables assumed in the native routines:
- statenames=c("S","I","R","cases","W"),
- ## the order of the parameters assumed in the native routines:
- paramnames=c(
- "gamma","mu","iota","beta1","beta.sd","pop","rho",
- "nbasis","degree","period"
- ),
- ## reset cases to zero after each new observation:
- zeronames=c("cases"),
- logvar=c(
- "gamma","mu","iota",
- "beta1","beta2","beta3","beta.sd",
- "S.0","I.0","R.0"
- ),
- logitvar="rho",
- comp.names=c("S","I","R"),
- ic.names=c("S.0","I.0","R.0"),
- parameter.transform=function (params, logvar, logitvar, ic.names, ...) {
- params[logvar] <- exp(params[logvar])
- params[logitvar] <- plogis(params[logitvar])
- params[ic.names] <- params[ic.names]/sum(params[ic.names])
- params
- },
- parameter.inv.transform=function (params, logvar, logitvar, ic.names, ...) {
- params[ic.names] <- params[ic.names]/sum(params[ic.names])
- params[logvar] <- log(params[logvar])
- params[logitvar] <- qlogis(params[logitvar])
- params
- },
- initializer=function(params, t0, comp.names, ic.names, ...) {
- x0 <- numeric(5)
- names(x0) <- c("S","I","R","cases","W")
- fracs <- params[ic.names]
- x0[comp.names] <- round(params['pop']*fracs/sum(fracs))
- x0
- }
- )
+ // compute the transition numbers
+ trans[0] = rpois(rate[0]*dt); // births are Poisson
+ reulermultinom(2,S,&rate[1],dt,&trans[1]);
+ reulermultinom(2,I,&rate[3],dt,&trans[3]);
+ reulermultinom(1,R,&rate[5],dt,&trans[5]);
- coef(po) <- c(
- gamma=26,mu=0.02,iota=0.01,
- 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-26/1200,
- nbasis=3,degree=3,period=1
- )
+ // balance the equations
+ S += trans[0]-trans[1]-trans[2];
+ I += trans[1]-trans[3]-trans[4];
+ R += trans[3]-trans[5];
+ cases += trans[3]; // cases are cumulative recoveries
+ if (beta_sd > 0.0) W += (dW-dt)/beta_sd; // increment has mean = 0, variance = dt
+'
+skel <- '
+ int nrate = 6;
+ double rate[nrate]; // transition rates
+ double term[nrate]; // transition numbers
+ double beta; // transmission rate
+ double dW; // white noise increment
+ double period = 1.0; // period of the seasonality
+ int nbasis = 3;
+ int deg = 3;
+ double seasonality[nbasis];
+ int k;
+
+ // compute transmission rate from seasonality
+ if (nbasis <= 0) return;
+ periodic_bspline_basis_eval(t,period,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
+ for (k = 0, beta = 0; k < nbasis; k++)
+ beta += (&beta1)[k]*seasonality[k];
- dyn.load(solib) ## load the shared-object library
+ // compute the transition rates
+ rate[0] = mu*popsize; // birth into susceptible class
+ rate[1] = (iota+beta*I)/popsize; // force of infection
+ rate[2] = mu; // death from susceptible class
+ rate[3] = gamma; // recovery
+ rate[4] = mu; // death from infectious class
+ rate[5] = mu; // death from recovered class
- ## compute a trajectory of the deterministic skeleton
- tic <- Sys.time()
- X <- trajectory(po,hmax=1/52,as.data.frame=TRUE)
- toc <- Sys.time()
- print(toc-tic)
+ // compute the several terms
+ term[0] = rate[0];
+ term[1] = rate[1]*S;
+ term[2] = rate[2]*S;
+ term[3] = rate[3]*I;
+ term[4] = rate[4]*I;
+ term[5] = rate[5]*R;
- plot(cases~time,data=X,type='l')
+ // assemble the differential equations
+ DS = term[0]-term[1]-term[2];
+ DI = term[1]-term[3]-term[4];
+ DR = term[3]-term[5];
+ Dcases = term[3]; // accumulate the new I->R transitions
+ DW = 0;
+'
- ## simulate from the model
- tic <- Sys.time()
- x <- simulate(po,nsim=3,as.data.frame=TRUE)
- toc <- Sys.time()
- print(toc-tic)
-
- plot(cases~time,data=x,col=as.factor(x$sim),pch=16)
+pompBuilder(
+ name="SIR",
+ data=data.frame(
+ time=seq(from=1/52,to=4,by=1/52),
+ reports=NA
+ ),
+ times="time",
+ t0=0,
+ dmeasure=dmeas,
+ rmeasure=rmeas,
+ step.fn=step.fn,
+ step.fn.delta.t=1/52/20,
+ skeleton.type="vectorfield",
+ skeleton=skel,
+ statenames=c("S","I","R","cases","W"),
+ paramnames=c(
+ "gamma","mu","iota","beta1","beta.sd","popsize","rho"
+ ),
+ zeronames=c("cases","W"),
+ logvar=c(
+ "gamma","mu","iota",
+ "beta1","beta2","beta3","beta.sd",
+ "S.0","I.0","R.0"
+ ),
+ logitvar="rho",
+ comp.names=c("S","I","R"),
+ ic.names=c("S.0","I.0","R.0"),
+ parameter.transform=function (params, logvar, logitvar, ic.names, ...) {
+ params[logvar] <- exp(params[logvar])
+ params[logitvar] <- plogis(params[logitvar])
+ params[ic.names] <- params[ic.names]/sum(params[ic.names])
+ params
+ },
+ parameter.inv.transform=function (params, logvar, logitvar, ic.names, ...) {
+ params[ic.names] <- params[ic.names]/sum(params[ic.names])
+ params[logvar] <- log(params[logvar])
+ params[logitvar] <- qlogis(params[logitvar])
+ params
+ },
+ initializer=function(params, t0, comp.names, ic.names, ...) {
+ x0 <- numeric(5)
+ names(x0) <- c("S","I","R","cases","W")
+ fracs <- params[ic.names]
+ x0[comp.names] <- round(params['popsize']*fracs/sum(fracs))
+ x0
+ }
+ ) -> po
- dyn.unload(solib)
+coef(po) <- c(
+ gamma=26,mu=0.02,iota=0.01,
+ beta1=120,beta2=140,beta3=100,
+ beta.sd=1e-3,
+ popsize=5e5,
+ rho=0.6,
+ S.0=26/120,I.0=0.001,R.0=1-26/120
+ )
-}
+## compute a trajectory of the deterministic skeleton
+tic <- Sys.time()
+X <- trajectory(po,hmax=1/52,as.data.frame=TRUE)
+toc <- Sys.time()
+print(toc-tic)
+
+plot(cases~time,data=X,type='l')
+
+## simulate from the model
+tic <- Sys.time()
+x <- simulate(po,nsim=3,as.data.frame=TRUE)
+toc <- Sys.time()
+print(toc-tic)
+
+plot(cases~time,data=x,col=as.factor(x$sim),pch=16)
+
+dyn.unload("SIR.so")
More information about the pomp-commits
mailing list