[Pomp-commits] r223 - in pkg: R data inst inst/data-R inst/doc inst/examples src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 10 15:50:37 CEST 2010
Author: kingaa
Date: 2010-05-10 15:50:37 +0200 (Mon, 10 May 2010)
New Revision: 223
Added:
pkg/data/gillespie.sir.rda
pkg/inst/data-R/gillespie.sir.R
pkg/inst/examples/sir.c
pkg/src/SSA.f
pkg/src/SSA_wrapper.c
pkg/src/sir.c
Removed:
pkg/inst/examples/euler_sir.c
pkg/src/euler_sir.c
pkg/src/euler_sir_density.c
Modified:
pkg/R/euler.R
pkg/data/euler.sir.rda
pkg/data/ou2.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/NEWS
pkg/inst/data-R/euler.sir.R
pkg/inst/data-R/ou2.R
pkg/inst/data-R/rw2.R
pkg/inst/data-R/verhulst.R
pkg/inst/doc/advanced_topics_in_pomp.Rnw
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/examples/euler_sir.R
pkg/inst/examples/logistic.R
pkg/inst/examples/rw2.R
pkg/src/euler.c
pkg/src/pomp.h
Log:
- encoding models is now much easier through the use of "plugin" modules.
- old plugins are now deprecated.
- several example pomps are now loadable via the data() command.
- new gillespie simulator (most of the work due to Helen Wearing)
Modified: pkg/R/euler.R
===================================================================
--- pkg/R/euler.R 2010-05-10 13:45:57 UTC (rev 222)
+++ pkg/R/euler.R 2010-05-10 13:50:37 UTC (rev 223)
@@ -6,6 +6,7 @@
zeronames = character(0),
tcovar, covar, PACKAGE)
{
+ .Deprecated(new="onestep.sim",package="pomp")
if (is.character(step.fun)) {
efun <- try(
getNativeSymbolInfo(step.fun,PACKAGE)$address,
@@ -48,6 +49,7 @@
zeronames = character(0),
tcovar, covar, PACKAGE)
{
+ .Deprecated(new="euler.sim",package="pomp")
if (is.character(step.fun)) {
efun <- try(
getNativeSymbolInfo(step.fun,PACKAGE)$address,
@@ -90,6 +92,7 @@
tcovar, covar, log = FALSE,
PACKAGE)
{
+ .Deprecated(new="onestep.dens",package="pomp")
if (is.character(dens.fun)) {
efun <- try(
getNativeSymbolInfo(dens.fun,PACKAGE)$address,
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Added: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Property changes on: pkg/data/gillespie.sir.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/verhulst.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2010-05-10 13:45:57 UTC (rev 222)
+++ pkg/inst/NEWS 2010-05-10 13:50:37 UTC (rev 223)
@@ -1,3 +1,11 @@
+NEW IN VERSION 0.29-1:
+
+ o big improvements to the documentation, including worked examples in vignettes.
+
+ o encoding models is now much easier through the use of "plugin" modules.
+
+ o several example pomps are now loadable via the data() command.
+
NEW IN VERSION 0.22-1:
o the nonlinear forecasting method of Kendall, Ellner, et al. is now implemented in the package. This is an "indirect inference" method using quasi-likelihood as an objective function.
Modified: pkg/inst/data-R/euler.sir.R
===================================================================
--- pkg/inst/data-R/euler.sir.R 2010-05-10 13:45:57 UTC (rev 222)
+++ pkg/inst/data-R/euler.sir.R 2010-05-10 13:50:37 UTC (rev 223)
@@ -5,7 +5,6 @@
times=seq(1/52,4,by=1/52),
data=rbind(measles=numeric(52*4)),
t0=0,
- delta.t=1/52/20,
statenames=c("S","I","R","cases","W"),
paramnames=c(
"gamma","mu","iota",
@@ -14,8 +13,11 @@
),
zeronames=c("cases"),
comp.names=c("S","I","R"),
- step.fun="sir_euler_simulator",
- rprocess=euler.simulate,
+ rprocess=euler.sim(
+ step.fun="sir_euler_simulator",
+ delta.t=1/52/20,
+ PACKAGE="pomp"
+ ),
skeleton.vectorfield="sir_ODE",
rmeasure="binom_rmeasure",
dmeasure="binom_dmeasure",
Added: pkg/inst/data-R/gillespie.sir.R
===================================================================
--- pkg/inst/data-R/gillespie.sir.R (rev 0)
+++ pkg/inst/data-R/gillespie.sir.R 2010-05-10 13:50:37 UTC (rev 223)
@@ -0,0 +1,69 @@
+library(pomp)
+
+params <- c(
+ nu=log(1/70),
+ mu=log(1/70),
+ beta1=log(330),
+ beta2=log(410),
+ beta3=log(490),
+ gamma=log(24),
+ iota=log(0.1),
+ S.0=log(0.05),
+ I.0=log(1e-4),
+ R.0=log(0.95),
+ N.0=1000000,
+ cases.0=0,
+ nbasis=3,
+ degree=3,
+ period=1
+ )
+
+simulate(
+ pomp(
+ data=data.frame(
+ time=seq(from=0,to=10,by=1/52),
+ reports=NA
+ ),
+ times="time",
+ t0=0,
+ rprocess=gillespie.sim(
+ rate.fun="sir_rates",
+ PACKAGE="pomp",
+ v=cbind(
+ birth=c(1,0,0,1,0),
+ sdeath=c(-1,0,0,-1,0),
+ infection=c(-1,1,0,0,0),
+ ideath=c(0,-1,0,-1,0),
+ recovery=c(0,-1,1,0,1),
+ rdeath=c(0,0,-1,-1,0)
+ ),
+ d=cbind(
+ birth=c(0,0,0,1,0),
+ sdeath=c(1,0,0,0,0),
+ infection=c(1,1,0,1,0),
+ ideath=c(0,1,0,0,0),
+ recovery=c(0,1,0,0,0),
+ rdeath=c(0,0,1,0,0)
+ )
+ ),
+ statenames=c("S","I","R","N","cases"),
+ paramnames=c("gamma","mu","iota","beta1","nu","nbasis","degree","period"),
+ zeronames=c("cases"),
+ measurement.model=reports~binom(size=cases,prob=0.1),
+ initializer=function(params, t0, ...){
+ comp.names <- c("S","I","R")
+ icnames <- paste(comp.names,"0",sep=".")
+ snames <- c("S","I","R","N","cases")
+ fracs <- exp(params[icnames])
+ x0 <- numeric(length(snames))
+ names(x0) <- snames
+ x0["N"] <- params["N.0"]
+ x0[comp.names] <- round(params['N.0']*fracs/sum(fracs))
+ x0
+ }
+ ),
+ params=params,
+ nsim=1,
+ seed=1165270654L
+ ) -> gillespie.sir
+
Modified: pkg/inst/data-R/ou2.R
===================================================================
--- pkg/inst/data-R/ou2.R 2010-05-10 13:45:57 UTC (rev 222)
+++ pkg/inst/data-R/ou2.R 2010-05-10 13:50:37 UTC (rev 223)
@@ -75,9 +75,10 @@
statenames = c("x1","x2")
),
params=c(
- alpha.1=0.9,alpha.2=0,alpha.3=0,alpha.4=0.99,
- sigma.1=1,sigma.2=0,sigma.3=2,
- tau=1,x1.0=50,x2.0=-50
+ alpha.1=0.8, alpha.2=-0.5, alpha.3=0.3, alpha.4=0.9,
+ sigma.1=3, sigma.2=-0.5, sigma.3=2,
+ tau=1,
+ x1.0=-3, x2.0=4
),
nsim=1,
seed=377456545L
Modified: pkg/inst/data-R/rw2.R
===================================================================
--- pkg/inst/data-R/rw2.R 2010-05-10 13:45:57 UTC (rev 222)
+++ pkg/inst/data-R/rw2.R 2010-05-10 13:50:37 UTC (rev 223)
@@ -2,25 +2,27 @@
simulate(
pomp(
- rprocess = onestep.simulate,
- dprocess = onestep.density,
- step.fun = function(x, t, params, delta.t, ...) {
- c(
- x1=rnorm(n=1,mean=x['x1'],sd=params['s1']*delta.t),
- x2=rnorm(n=1,mean=x['x2'],sd=params['s2']*delta.t)
- )
- },
- dens.fun = function (x1, t1, x2, t2, params, ...) {
- sum(
- dnorm(
- x=x2[c('x1','x2')],
- mean=x1[c('x1','x2')],
- sd=params[c('s1','s2')]*(t2-t1),
- log=TRUE
- ),
- na.rm=TRUE
+ rprocess = onestep.sim(
+ step.fun = function(x, t, params, delta.t, ...) {
+ c(
+ x1=rnorm(n=1,mean=x['x1'],sd=params['s1']*delta.t),
+ x2=rnorm(n=1,mean=x['x2'],sd=params['s2']*delta.t)
)
- },
+ }
+ ),
+ dprocess = onestep.dens(
+ dens.fun = function (x1, t1, x2, t2, params, ...) {
+ sum(
+ dnorm(
+ x=x2[c('x1','x2')],
+ mean=x1[c('x1','x2')],
+ sd=params[c('s1','s2')]*(t2-t1),
+ log=TRUE
+ ),
+ na.rm=TRUE
+ )
+ }
+ ),
measurement.model=list(
y1 ~ norm(mean=x1,sd=tau),
y2 ~ norm(mean=x2,sd=tau)
Modified: pkg/inst/data-R/verhulst.R
===================================================================
--- pkg/inst/data-R/verhulst.R 2010-05-10 13:45:57 UTC (rev 222)
+++ pkg/inst/data-R/verhulst.R 2010-05-10 13:50:37 UTC (rev 223)
@@ -5,38 +5,40 @@
data=rbind(obs=rep(0,1000)),
times=seq(0.1,by=0.1,length=1000),
t0=0,
- rprocess=euler.simulate,
- step.fun=function(x,t,params,delta.t,...){
- with(
- as.list(c(x,params)),
- rnorm(
- n=1,
- mean=n+r*n*(1-n/K)*delta.t,
- sd=sigma*n*sqrt(delta.t)
- )
- )
- },
- dprocess=onestep.density,
- dens.fun=function(x1,x2,t1,t2,params,log,...){
- delta.t <- t2-t1
- with(
- as.list(c(x1,params)),
- dnorm(
- x=x2['n'],
- mean=n+r*n*(1-n/K)*delta.t,
- sd=sigma*n*sqrt(delta.t),
- log=log
- )
- )
- },
+ rprocess=euler.sim(
+ step.fun=function(x,t,params,delta.t,...){
+ with(
+ as.list(c(x,params)),
+ rnorm(
+ n=1,
+ mean=n+r*n*(1-n/K)*delta.t,
+ sd=sigma*n*sqrt(delta.t)
+ )
+ )
+ },
+ delta.t=0.01
+ ),
+ dprocess=onestep.dens(
+ dens.fun=function(x1,x2,t1,t2,params,log,...){
+ delta.t <- t2-t1
+ with(
+ as.list(c(x1,params)),
+ dnorm(
+ x=x2['n'],
+ mean=n+r*n*(1-n/K)*delta.t,
+ sd=sigma*n*sqrt(delta.t),
+ log=log
+ )
+ )
+ }
+ ),
measurement.model=obs~lnorm(meanlog=log(n),sdlog=log(1+tau)),
skeleton.vectorfield=function(x,t,params,...){
with(
as.list(c(x,params)),
r*n*(1-n/K)
)
- },
- delta.t=0.01
+ }
),
params=c(
n.0=10000,
Modified: pkg/inst/doc/advanced_topics_in_pomp.Rnw
===================================================================
--- pkg/inst/doc/advanced_topics_in_pomp.Rnw 2010-05-10 13:45:57 UTC (rev 222)
+++ pkg/inst/doc/advanced_topics_in_pomp.Rnw 2010-05-10 13:50:37 UTC (rev 223)
@@ -242,13 +242,14 @@
times=seq(1/52,4,by=1/52),
data=rbind(measles=numeric(52*4)),
t0=0,
- delta.t=1/52/20,
statenames=c("S","I","R","cases","W"),
paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho","nbasis","degree","period"),
zeronames=c("cases"),
comp.names=c("S","I","R"),
- step.fun="sir_euler_simulator",
- rprocess=euler.simulate,
+ rprocess=euler.sim(
+ step.fun="sir_euler_simulator",
+ delta.t=1/52/20
+ ),
skeleton.vectorfield="sir_ODE",
rmeasure="binom_rmeasure",
dmeasure="binom_dmeasure",
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/examples/euler_sir.R
===================================================================
--- pkg/inst/examples/euler_sir.R 2010-05-10 13:45:57 UTC (rev 222)
+++ pkg/inst/examples/euler_sir.R 2010-05-10 13:50:37 UTC (rev 223)
@@ -22,33 +22,35 @@
t0=0,
tcovar=tbasis,
covar=basis,
- delta.t=1/52/20,
zeronames=c("cases"),
- step.fun=function(t,x,params,covars,delta.t,...) {
- params <- exp(params)
- with(
- as.list(c(x,params)),
- {
- beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
- beta.var <- beta.sd^2
- dW <- rgamma(n=1,shape=delta.t/beta.var,scale=beta.var)
- foi <- (iota+beta*I*dW/delta.t)/pop
- trans <- c(
- rpois(n=1,lambda=mu*pop*delta.t),
- reulermultinom(n=1,size=S,rate=c(foi,mu),dt=delta.t),
- reulermultinom(n=1,size=I,rate=c(gamma,mu),dt=delta.t),
- reulermultinom(n=1,size=R,rate=c(mu),dt=delta.t)
- )
- c(
- S=S+trans[1]-trans[2]-trans[3],
- I=I+trans[2]-trans[4]-trans[5],
- R=R+trans[4]-trans[6],
- cases=cases+trans[4],
- W=if (beta.sd>0) W+(dW-delta.t)/beta.sd else W
- )
- }
- )
- },
+ rprocess=euler.sim(
+ step.fun=function(t,x,params,covars,delta.t,...) {
+ params <- exp(params)
+ with(
+ as.list(c(x,params)),
+ {
+ beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
+ beta.var <- beta.sd^2
+ dW <- rgamma(n=1,shape=delta.t/beta.var,scale=beta.var)
+ foi <- (iota+beta*I*dW/delta.t)/pop
+ trans <- c(
+ rpois(n=1,lambda=mu*pop*delta.t),
+ reulermultinom(n=1,size=S,rate=c(foi,mu),dt=delta.t),
+ reulermultinom(n=1,size=I,rate=c(gamma,mu),dt=delta.t),
+ reulermultinom(n=1,size=R,rate=c(mu),dt=delta.t)
+ )
+ c(
+ S=S+trans[1]-trans[2]-trans[3],
+ I=I+trans[2]-trans[4]-trans[5],
+ R=R+trans[4]-trans[6],
+ cases=cases+trans[4],
+ W=if (beta.sd>0) W+(dW-delta.t)/beta.sd else W
+ )
+ }
+ )
+ },
+ delta.t=1/52/20
+ ),
skeleton.vectorfield=function(x,t,params,covars,...) {
params <- exp(params)
with(
@@ -74,7 +76,6 @@
}
)
},
- rprocess=euler.simulate,
measurement.model=measles~binom(size=cases,prob=exp(rho)),
initializer=function(params,t0,...){
p <- exp(params)
@@ -94,11 +95,12 @@
)
## alternatively, one can define the computationally intensive bits using native routines:
-## the C codes "sir_euler_simulator" and "sir_euler_density" are included in the "examples" directory (file "euler_sir.c")
+## the C codes "sir_euler_simulator", "sir_ODE", "binom_rmeasure", and "binom_dmeasure"
+## are included in the "examples" directory (file "sir.c")
if (Sys.info()['sysname']=='Linux') {
- model <- "euler_sir"
+ model <- "sir"
pkg <- "pomp"
modelfile <- paste(model,".c",sep="")
headerfile <- system.file("include/pomp.h",package=pkg)
@@ -118,16 +120,17 @@
times=seq(1/52,4,by=1/52),
data=rbind(measles=numeric(52*4)),
t0=0,
- delta.t=1/52/20,
statenames=c("S","I","R","cases","W"),
paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho","nbasis","degree","period"),
zeronames=c("cases"),
- step.fun="sir_euler_simulator",
- rprocess=euler.simulate,
+ rprocess=euler.sim(
+ step.fun="sir_euler_simulator",
+ delta.t=1/52/20
+ ),
skeleton.vectorfield="sir_ODE",
rmeasure="binom_rmeasure",
dmeasure="binom_dmeasure",
- PACKAGE="euler_sir", ## name of the shared-object library
+ PACKAGE="sir", ## name of the shared-object library
initializer=function(params,t0,...){
p <- exp(params)
with(
Deleted: pkg/inst/examples/euler_sir.c
===================================================================
--- pkg/inst/examples/euler_sir.c 2010-05-10 13:45:57 UTC (rev 222)
+++ pkg/inst/examples/euler_sir.c 2010-05-10 13:50:37 UTC (rev 223)
@@ -1,208 +0,0 @@
-// dear emacs, please treat this as -*- C++ -*-
-
-#include <Rmath.h>
-
-#include "pomp.h"
-
-static double expit (double x) {
- return 1.0/(1.0 + exp(-x));
-}
-
-static double logit (double x) {
- return log(x/(1-x));
-}
-
-static double term_time (double t, double b0, double b1)
-{
- static double correction = 0.52328767123287671233;
- double day = 365.0 * (t - floor(t));
- int tt = (day >= 7.0 && day <= 100.0)
- || (day >= 115.0 && day <= 200.0)
- || (day >= 252.0 && day <= 300.0)
- || (day >= 308.0 && day <= 356.0);
- return b0 * (1.0 + b1 * (-1.0 + 2.0 * ((double) tt))) / (1.0 + correction * b1);
-}
-
-#define LOGGAMMA (p[parindex[0]]) // recovery rate
-#define LOGMU (p[parindex[1]]) // baseline birth and death rate
-#define LOGIOTA (p[parindex[2]]) // import rate
-#define LOGBETA (p[parindex[3]]) // transmission rate
-#define LOGBETA_SD (p[parindex[4]]) // environmental stochasticity SD in transmission rate
-#define LOGPOPSIZE (p[parindex[5]]) // population size
-#define LOGRHO (p[parindex[6]]) // reporting probability
-#define NBASIS (p[parindex[7]]) // number of periodic B-spline basis functions
-#define DEGREE (p[parindex[8]]) // degree of periodic B-spline basis functions
-#define PERIOD (p[parindex[9]]) // period of B-spline basis functions
-
-#define SUSC (x[stateindex[0]]) // number of susceptibles
-#define INFD (x[stateindex[1]]) // number of infectives
-#define RCVD (x[stateindex[2]]) // number of recovereds
-#define CASE (x[stateindex[3]]) // number of cases (accumulated per reporting period)
-#define W (x[stateindex[4]]) // integrated white noise
-
-#define MEASLES (y[0])
-
-void binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
- int *stateindex, int *parindex, int *covindex,
- int ncovars, double *covars, double t) {
- *lik = dbinom(MEASLES,CASE,exp(LOGRHO),give_log);
-}
-
-void binom_rmeasure (double *y, double *x, double *p,
- int *stateindex, int *parindex, int *covindex,
- int ncovars, double *covars, double t) {
- MEASLES = rbinom(CASE,exp(LOGRHO));
-}
-
-#undef MEASLES
-
-// SIR model with Euler multinomial step
-// forced transmission (basis functions passed as covariates)
-// constant population size as a parameter
-// environmental stochasticity on transmission
-void sir_euler_simulator (double *x, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
- int covdim, const double *covar,
- double t, double dt)
-{
- int nrate = 6;
- double rate[nrate]; // transition rates
- double trans[nrate]; // transition numbers
- double gamma, mu, iota, beta_sd, beta_var, popsize;
- double beta;
- double dW;
- int nseas = (int) NBASIS; // number of seasonal basis functions
- int deg = (int) DEGREE; // degree of seasonal basis functions
- double seasonality[nseas];
-
- // untransform the parameters
- gamma = exp(LOGGAMMA);
- mu = exp(LOGMU);
- iota = exp(LOGIOTA);
- beta_sd = exp(LOGBETA_SD);
- popsize = exp(LOGPOPSIZE);
- beta_var = beta_sd*beta_sd;
-
- if (nseas <= 0) return;
- periodic_bspline_basis_eval(t,PERIOD,deg,nseas,&seasonality[0]);
- beta = exp(dot_product(nseas,&seasonality[0],&LOGBETA));
-
- // test to make sure the parameters and state variable values are sane
- if (!(R_FINITE(beta)) ||
- !(R_FINITE(gamma)) ||
- !(R_FINITE(mu)) ||
- !(R_FINITE(beta_sd)) ||
- !(R_FINITE(iota)) ||
- !(R_FINITE(popsize)) ||
- !(R_FINITE(SUSC)) ||
- !(R_FINITE(INFD)) ||
- !(R_FINITE(RCVD)) ||
- !(R_FINITE(CASE)) ||
- !(R_FINITE(W)))
- return;
-
- if (beta_sd > 0.0) { // environmental noise is ON
- dW = rgamma(dt/beta_var,beta_var); // gamma noise, mean=dt, variance=(beta_sd^2 dt)
- if (!(R_FINITE(dW))) return;
- } else { // environmental noise is OFF
- dW = dt;
- }
-
- // compute the transition rates
- rate[0] = mu*popsize; // birth into susceptible class
- rate[1] = (iota+beta*INFD*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
-
- // compute the transition numbers
- trans[0] = rpois(rate[0]*dt); // births are Poisson
- reulermultinom(2,SUSC,&rate[1],dt,&trans[1]);
- reulermultinom(2,INFD,&rate[3],dt,&trans[3]);
- reulermultinom(1,RCVD,&rate[5],dt,&trans[5]);
-
- // balance the equations
- SUSC += trans[0]-trans[1]-trans[2];
- INFD += trans[1]-trans[3]-trans[4];
- RCVD += trans[3]-trans[5];
- CASE += trans[3]; // cases are cumulative recoveries
- if (beta_sd > 0.0) {
- W += (dW-dt)/beta_sd; // mean zero, variance = dt
- }
-
-}
-
-#define DSDT (f[stateindex[0]])
-#define DIDT (f[stateindex[1]])
-#define DRDT (f[stateindex[2]])
-#define DCDT (f[stateindex[3]])
-
-void sir_ODE (double *f, double *x, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
- int covdim, const double *covar, double t)
-{
- int nrate = 6;
- double rate[nrate]; // transition rates
- double term[nrate]; // terms in the equations
- double gamma, mu, iota, popsize;
- double beta;
- int nseas = (int) NBASIS; // number of seasonal basis functions
- int deg = (int) DEGREE; // degree of seasonal basis functions
- double seasonality[nseas];
-
- // untransform the parameters
- gamma = exp(LOGGAMMA);
- mu = exp(LOGMU);
- iota = exp(LOGIOTA);
- popsize = exp(LOGPOPSIZE);
-
- if (nseas <= 0) return;
- periodic_bspline_basis_eval(t,PERIOD,deg,nseas,&seasonality[0]);
- beta = exp(dot_product(nseas,&seasonality[0],&LOGBETA));
-
- // compute the transition rates
- rate[0] = mu*popsize; // birth into susceptible class
- rate[1] = (iota+beta*INFD)/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 the several terms
- term[0] = rate[0];
- term[1] = rate[1]*SUSC;
- term[2] = rate[2]*SUSC;
- term[3] = rate[3]*INFD;
- term[4] = rate[4]*INFD;
- term[5] = rate[5]*RCVD;
-
- // balance the equations
- DSDT = term[0]-term[1]-term[2];
- DIDT = term[1]-term[3]-term[4];
- DRDT = term[3]-term[5];
- DCDT = term[3]; // cases are cumulative recoveries
-
-}
-
-#undef DSDT
-#undef DIDT
-#undef DRDT
-#undef DCDT
-
-#undef SUSC
-#undef INFD
-#undef RCVD
-#undef CASE
-#undef W
-
-#undef LOGGAMMA
-#undef LOGMU
-#undef LOGIOTA
-#undef LOGBETA
-#undef LOGBETA_SD
-#undef LOGPOPSIZE
-#undef LOGRHO
-#undef NBASIS
-#undef DEGREE
-#undef PERIOD
Modified: pkg/inst/examples/logistic.R
===================================================================
--- pkg/inst/examples/logistic.R 2010-05-10 13:45:57 UTC (rev 222)
+++ pkg/inst/examples/logistic.R 2010-05-10 13:50:37 UTC (rev 223)
@@ -1,41 +1,45 @@
require(pomp)
+## a stochastic version of the Verhulst-Pearl logistic model
+
po <- pomp(
data=rbind(obs=rep(0,1000)),
times=seq(0.1,by=0.1,length=1000),
t0=0,
- rprocess=euler.simulate,
- step.fun=function(x,t,params,delta.t,...){
- with(
- as.list(c(x,params)),
- rnorm(
- n=1,
- mean=n+r*n*(1-n/K)*delta.t,
- sd=sigma*n*sqrt(delta.t)
- )
- )
- },
- dprocess=onestep.density,
- dens.fun=function(x1,x2,t1,t2,params,log,...){
- delta.t <- t2-t1
- with(
- as.list(c(x1,params)),
- dnorm(
- x=x2['n'],
- mean=n+r*n*(1-n/K)*delta.t,
- sd=sigma*n*sqrt(delta.t),
- log=log
- )
- )
- },
+ rprocess=euler.sim(
+ step.fun=function(x,t,params,delta.t,...){
+ with(
+ as.list(c(x,params)),
+ rnorm(
+ n=1,
+ mean=n+r*n*(1-n/K)*delta.t,
+ sd=sigma*n*sqrt(delta.t)
+ )
+ )
+ },
+ delta.t=0.01
+ ),
+ dprocess=onestep.dens(
+ dens.fun=function(x1,x2,t1,t2,params,log,...){
+ delta.t <- t2-t1
+ with(
+ as.list(c(x1,params)),
+ dnorm(
+ x=x2['n'],
+ mean=n+r*n*(1-n/K)*delta.t,
+ sd=sigma*n*sqrt(delta.t),
+ log=log
+ )
+ )
+ }
+ ),
measurement.model=obs~lnorm(meanlog=log(n),sdlog=log(1+tau)),
skeleton.vectorfield=function(x,t,params,...){
with(
as.list(c(x,params)),
r*n*(1-n/K)
)
- },
- delta.t=0.01
+ }
)
params <- c(n.0=10000,K=10000,r=0.9,sigma=0.4,tau=0.1)
Modified: pkg/inst/examples/rw2.R
===================================================================
--- pkg/inst/examples/rw2.R 2010-05-10 13:45:57 UTC (rev 222)
+++ pkg/inst/examples/rw2.R 2010-05-10 13:50:37 UTC (rev 223)
@@ -1,27 +1,30 @@
require(pomp)
-## using the "onestep" plugins
+## a simple two-dimensional random walk
rw2 <- pomp(
- rprocess = onestep.simulate,
- dprocess = onestep.density,
- step.fun = function(x, t, params, delta.t, ...) {
- c(
- x1=rnorm(n=1,mean=x['x1'],sd=params['s1']*delta.t),
- x2=rnorm(n=1,mean=x['x2'],sd=params['s2']*delta.t)
- )
- },
- dens.fun = function (x1, t1, x2, t2, params, ...) {
- sum(
- dnorm(
- x=x2[c('x1','x2')],
- mean=x1[c('x1','x2')],
- sd=params[c('s1','s2')]*(t2-t1),
- log=TRUE
- ),
- na.rm=TRUE
+ rprocess = euler.sim(
+ step.fun = function(x, t, params, delta.t, ...) {
+ c(
+ x1=rnorm(n=1,mean=x['x1'],sd=params['s1']*delta.t),
+ x2=rnorm(n=1,mean=x['x2'],sd=params['s2']*delta.t)
)
- },
+ },
+ delta.t=1
+ ),
+ dprocess = onestep.dens(
+ dens.fun = function (x1, t1, x2, t2, params, ...) {
+ sum(
+ dnorm(
+ x=x2[c('x1','x2')],
+ mean=x1[c('x1','x2')],
+ sd=params[c('s1','s2')]*(t2-t1),
+ log=TRUE
+ ),
+ na.rm=TRUE
+ )
+ }
+ ),
measurement.model=list(
y1 ~ norm(mean=x1,sd=tau),
y2 ~ norm(mean=x2,sd=tau)
@@ -33,78 +36,3 @@
),
t0=0
)
-
-## writing rprocess and dprocess from scratch
-
-rw.rprocess <- function (params, xstart, times, ...) {
- ## this function simulates two independent random walks with intensities s1, s2
- nvars <- nrow(xstart)
- nreps <- ncol(params)
- ntimes <- length(times)
- dt <- diff(times)
- x <- array(0,dim=c(nvars,nreps,ntimes))
- rownames(x) <- rownames(xstart)
- noise.sds <- params[c('s1','s2'),]
- x[,,1] <- xstart
- for (j in 2:ntimes) {
- ## we are mimicking a continuous-time process, so the increments have SD ~ sqrt(dt)
- ## note that we do not have to assume that 'times' are equally spaced
- x[c("x1","x2"),,j] <- rnorm(
- n=2*nreps,
- mean=x[c("x1","x2"),,j-1],
- sd=noise.sds*dt[j-1]
- )
- }
- x
-}
-
-rw.dprocess <- function (x, times, params, log = FALSE, ...) {
- ## given a sequence of consecutive states in 'x', this function computes the p.d.f.
- nreps <- ncol(params)
- ntimes <- length(times)
- dt <- diff(times)
- d <- array(0,dim=c(2,nreps,ntimes-1))
- noise.sds <- params[c('s1','s2'),]
- for (j in 2:ntimes)
- d[,,j-1] <- dnorm(x[,,j]-x[,,j-1],mean=0,sd=noise.sds*dt[j-1],log=TRUE)
- d <- apply(d,c(2,3),sum)
- if (log) d else exp(d)
-}
-
-bvnorm.rmeasure <- function (t, x, params, ...) {
- ## noisy observations of the two walks with common noise SD 'tau'
- c(
- y1=rnorm(n=1,mean=x['x1'],sd=params['tau']),
- y2=rnorm(n=1,mean=x['x2'],sd=params['tau'])
- )
-}
-
-bvnorm.dmeasure <- function (y, x, t, params, log = FALSE, ...) {
- f <- sum(
- dnorm(
- x=y[c("y1","y2")],
- mean=x[c("x1","x2")],
- sd=params["tau"],
- log=TRUE
- ),
- na.rm=TRUE
- )
- if (log) f else exp(f)
-}
-
-rw2 <- pomp(
- rprocess = rw.rprocess,
- dprocess = rw.dprocess,
- measurement.model=list(
- y1 ~ norm(mean=x1,sd=tau),
- y2 ~ norm(mean=x2,sd=tau)
- ),
- times=1:100,
- data=rbind(
- y1=rep(0,100),
- y2=rep(0,100)
- ),
- t0=0,
- useless=23
- )
-
Copied: pkg/inst/examples/sir.c (from rev 219, pkg/inst/examples/euler_sir.c)
===================================================================
--- pkg/inst/examples/sir.c (rev 0)
+++ pkg/inst/examples/sir.c 2010-05-10 13:50:37 UTC (rev 223)
@@ -0,0 +1,208 @@
+// dear emacs, please treat this as -*- C++ -*-
+
+#include <Rmath.h>
+
+#include "pomp.h"
+
+static double expit (double x) {
+ return 1.0/(1.0 + exp(-x));
+}
+
+static double logit (double x) {
+ return log(x/(1-x));
+}
+
+static double term_time (double t, double b0, double b1)
+{
+ static double correction = 0.52328767123287671233;
+ double day = 365.0 * (t - floor(t));
+ int tt = (day >= 7.0 && day <= 100.0)
+ || (day >= 115.0 && day <= 200.0)
+ || (day >= 252.0 && day <= 300.0)
+ || (day >= 308.0 && day <= 356.0);
+ return b0 * (1.0 + b1 * (-1.0 + 2.0 * ((double) tt))) / (1.0 + correction * b1);
+}
+
+#define LOGGAMMA (p[parindex[0]]) // recovery rate
+#define LOGMU (p[parindex[1]]) // baseline birth and death rate
+#define LOGIOTA (p[parindex[2]]) // import rate
+#define LOGBETA (p[parindex[3]]) // transmission rate
+#define LOGBETA_SD (p[parindex[4]]) // environmental stochasticity SD in transmission rate
+#define LOGPOPSIZE (p[parindex[5]]) // population size
+#define LOGRHO (p[parindex[6]]) // reporting probability
+#define NBASIS (p[parindex[7]]) // number of periodic B-spline basis functions
+#define DEGREE (p[parindex[8]]) // degree of periodic B-spline basis functions
+#define PERIOD (p[parindex[9]]) // period of B-spline basis functions
+
+#define SUSC (x[stateindex[0]]) // number of susceptibles
+#define INFD (x[stateindex[1]]) // number of infectives
+#define RCVD (x[stateindex[2]]) // number of recovereds
+#define CASE (x[stateindex[3]]) // number of cases (accumulated per reporting period)
+#define W (x[stateindex[4]]) // integrated white noise
+
+#define MEASLES (y[0])
+
+void binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ *lik = dbinom(MEASLES,CASE,exp(LOGRHO),give_log);
+}
+
+void binom_rmeasure (double *y, double *x, double *p,
+ int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ MEASLES = rbinom(CASE,exp(LOGRHO));
+}
+
+#undef MEASLES
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 223
More information about the pomp-commits
mailing list