[Pomp-commits] r75 - in pkg: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 3 22:22:04 CET 2009
Author: kingaa
Date: 2009-02-03 22:22:04 +0100 (Tue, 03 Feb 2009)
New Revision: 75
Added:
pkg/src/euler_sir.c
Removed:
pkg/src/sir.c
Modified:
pkg/DESCRIPTION
pkg/R/mif.R
pkg/man/mif.Rd
pkg/man/pomp-methods.Rd
pkg/man/trajectory-pomp.Rd
Log:
changes to documentation and default arguments so as to bring these into agreement. (discrepancies revealed by Kurt Hornik's extensive codoc testing).
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2009-02-01 14:07:02 UTC (rev 74)
+++ pkg/DESCRIPTION 2009-02-03 21:22:04 UTC (rev 75)
@@ -1,11 +1,11 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.22-4
-Date: 2009-01-31
+Version: 0.22-5
+Date: 2009-02-03
Author: Aaron A. King, Edward L. Ionides, Carles Martinez Breto, Steve Ellner, Bruce Kendall
Maintainer: Aaron A. King <kingaa at umich.edu>
Description: Inference methods for partially-observed Markov processes
-Depends: R(>= 2.6.1), stats, methods, graphics, deSolve, subplex
+Depends: R(>= 2.7.1), stats, methods, graphics, deSolve, subplex
License: GPL (>= 2)
LazyLoad: true
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2009-02-01 14:07:02 UTC (rev 74)
+++ pkg/R/mif.R 2009-02-03 21:22:04 UTC (rev 75)
@@ -20,14 +20,18 @@
"pomp",
function (object, Nmif = 1,
start,
- pars = stop("mif error: ",sQuote("pars")," must be specified",call.=FALSE),
- ivps = character(0),
+ pars, ivps = character(0),
particles,
- rw.sd = stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE),
- alg.pars = stop("mif error: ",sQuote("alg.pars")," must be specified",call.=FALSE),
+ rw.sd, alg.pars,
weighted = TRUE, tol = 1e-17, warn = TRUE, max.fail = 0,
verbose = FALSE, .ndone = 0)
{
+ if (missing(pars))
+ stop("mif error: ",sQuote("pars")," must be specified",call.=FALSE)
+ if (missing(rw.sd))
+ stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
+ if (missing(alg.pars))
+ stop("mif error: ",sQuote("alg.pars")," must be specified",call.=FALSE)
if (missing(particles)) { # use default: normal distribution
particles <- function (Np, center, sd, ...) {
matrix(
@@ -152,12 +156,16 @@
setMethod(
"mif",
"mif",
- function (object, Nmif = object at Nmif, start = coef(object),
- pars = object at pars, ivps = object at ivps, rw.sd = object at random.walk.sd,
- alg.pars = object at alg.pars,
+ function (object, Nmif, start, pars, ivps, rw.sd, alg.pars,
weighted = TRUE, tol = 1e-17, warn = TRUE, max.fail = 0,
verbose = FALSE, .ndone = 0)
{
+ if (missing(Nmif)) Nmif <- object at Nmif
+ if (missing(start)) start <- coef(object)
+ if (missing(pars)) pars <- object at pars
+ if (missing(ivps)) ivps <- object at ivps
+ if (missing(rw.sd)) rw.sd <- object at random.walk.sd
+ if (missing(alg.pars)) alg.pars <- object at alg.pars
theta <- start
start.names <- names(start)
if (is.null(start.names))
Modified: pkg/man/mif.Rd
===================================================================
--- pkg/man/mif.Rd 2009-02-01 14:07:02 UTC (rev 74)
+++ pkg/man/mif.Rd 2009-02-03 21:22:04 UTC (rev 75)
@@ -14,10 +14,12 @@
}
\usage{
mif(object, \dots)
-\S4method{mif}{pomp}(object, Nmif, start, pars, ivps = character(0),
+\S4method{mif}{pomp}(object, Nmif = 1, start, pars, ivps = character(0),
particles, rw.sd, alg.pars, weighted = TRUE,
- tol = 1e-17, warn = TRUE, max.fail = 0, verbose = FALSE)
-\S4method{mif}{mif}(object, \dots)
+ tol = 1e-17, warn = TRUE, max.fail = 0, verbose = FALSE, .ndone = 0)
+\S4method{mif}{mif}(object, Nmif, start, pars, ivps, rw.sd, alg.pars,
+ weighted = TRUE, tol = 1e-17, warn = TRUE, max.fail = 0,
+ verbose = FALSE, .ndone = 0)
\S4method{continue}{mif}(object, Nmif, \dots)
}
\arguments{
@@ -92,6 +94,10 @@
\item{verbose}{
logical; if TRUE, print progress reports.
}
+ \item{.ndone}{
+ for internal use by \code{continue}.
+ do not meddle with this!
+ }
\item{\dots}{
Additional arguments that can be used to override the defaults.
}
Modified: pkg/man/pomp-methods.Rd
===================================================================
--- pkg/man/pomp-methods.Rd 2009-02-01 14:07:02 UTC (rev 74)
+++ pkg/man/pomp-methods.Rd 2009-02-03 21:22:04 UTC (rev 75)
@@ -30,18 +30,18 @@
coef(object, pars, \dots) <- value
\S4method{data.array}{pomp}(object, vars, \dots)
\S4method{states}{pomp}(object, vars, \dots)
-\S4method{time}{pomp}(object, t0 = FALSE)
+\S4method{time}{pomp}(x, t0 = FALSE, \dots)
\S4method{show}{pomp}(object)
\S4method{as}{pomp}(object, class)
-\S4method{coerce}{pomp,data.frame}(from, to = "data.frame")
+\S4method{coerce}{pomp,data.frame}(from, to = "data.frame", strict = TRUE)
\S4method{print}{pomp}(x, \dots)
-\S4method{plot}{pomp}(x, y = NULL, variables, panel = lines,
+\S4method{plot}{pomp}(x, y, variables, panel = lines,
nc = NULL, yax.flip = FALSE,
mar = c(0, 5.1, 0, if (yax.flip) 5.1 else 2.1),
oma = c(6, 0, 5, 0), axes = TRUE, \dots)
}
\arguments{
- \item{object}{The \code{pomp} object.}
+ \item{object, x}{The \code{pomp} object.}
\item{pars}{
optional character;
names of parameters to be retrieved or set.
@@ -65,7 +65,9 @@
\item{from, to}{
the classes betwen which coercion should be performed.
}
- \item{x}{the \code{pomp} object.}
+ \item{strict}{
+ ignored.
+ }
\item{y}{ignored.}
\item{variables}{
optional character;
Modified: pkg/man/trajectory-pomp.Rd
===================================================================
--- pkg/man/trajectory-pomp.Rd 2009-02-01 14:07:02 UTC (rev 74)
+++ pkg/man/trajectory-pomp.Rd 2009-02-03 21:22:04 UTC (rev 75)
@@ -12,13 +12,14 @@
}
\usage{
trajectory(object, params, times, \dots)
-\S4method{trajectory}{pomp}(object, params, times = time(object), \dots)
+\S4method{trajectory}{pomp}(object, params, times, \dots)
}
\arguments{
\item{object}{an object of class \code{pomp}.}
\item{times}{
a numeric vector containing the times at which a trajectory is desired.
The first of these will be the initial time.
+ By default, \code{times=time(object,t0=TRUE)}.
}
\item{params}{
a rank-2 array of parameters.
Copied: pkg/src/euler_sir.c (from rev 74, pkg/src/sir.c)
===================================================================
--- pkg/src/euler_sir.c (rev 0)
+++ pkg/src/euler_sir.c 2009-02-03 21:22:04 UTC (rev 75)
@@ -0,0 +1,265 @@
+// 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 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 BIRTHS (x[stateindex[5]]) // births
+#define dW (x[stateindex[6]]) // white noise process
+
+#define SEASBASIS (covar[covindex[0]]) // first column of seasonality basis in lookup table
+
+// 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; // transition numbers
+ double gamma, mu, iota, beta_sd, beta_var, popsize;
+ double beta;
+
+ // 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;
+
+ beta = exp(dot_product(covdim,&SEASBASIS,&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;
+
+ GetRNGstate(); // initialize R's random number generator
+
+ 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 = &BIRTHS;
+ 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]);
+
+ PutRNGstate(); // finished with R's random number generator
+
+ // 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
+ }
+
+}
+
+#undef BIRTHS
+#undef W
+#undef dW
+
+#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;
+
+ // untransform the parameters
+ gamma = exp(LOGGAMMA);
+ mu = exp(LOGMU);
+ iota = exp(LOGIOTA);
+ popsize = exp(LOGPOPSIZE);
+
+ beta = exp(dot_product(covdim,&SEASBASIS,&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
+
+#define SUSC (x1[stateindex[0]]) // number of susceptibles
+#define INFD (x1[stateindex[1]]) // number of infectives
+#define RCVD (x1[stateindex[2]]) // number of recovereds
+#define CASE (x1[stateindex[3]]) // number of cases (accumulated per reporting period)
+#define W (x1[stateindex[4]]) // integrated white noise
+#define BIRTHS (x2[stateindex[5]]) // births
+#define dW (x2[stateindex[6]]) // white noise process
+
+// 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_density (double *f, double *x1, double *x2, double t1, double t2, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar)
+{
+ int nrate = 6;
+ double rate[nrate]; // transition rates
+ double *trans; // transition numbers
+ double gamma, mu, iota, beta_sd, popsize;
+ double beta;
+ double dt = t2-t1;
+
+ // untransform the parameters
+ gamma = exp(LOGGAMMA);
+ mu = exp(LOGMU);
+ iota = exp(LOGIOTA);
+ beta_sd = exp(LOGBETA_SD);
+ popsize = exp(LOGPOPSIZE);
+
+ beta = exp(dot_product(covdim,&SEASBASIS,&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))) {
+ *f = R_NaN;
+ return;
+ }
+
+ // compute the transition rates
+ trans = &BIRTHS;
+ if (beta_sd > 0.0) { // environmental noise is ON
+ double beta_var = beta_sd*beta_sd;
+ *f = dgamma(dW,dt/beta_var,beta_var,1);
+ } else { // environmental noise is OFF
+ *f = 0; // THIS ASSUMES THAT dw = dt !!!
+ }
+ 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
+ *f += dpois(trans[0],rate[0]*dt,1); // births are Poisson
+ *f += deulermultinom(2,SUSC,&rate[1],dt,&trans[1],1);
+ *f += deulermultinom(2,INFD,&rate[3],dt,&trans[3],1);
+ *f += deulermultinom(1,RCVD,&rate[5],dt,&trans[5],1);
+}
+
+#undef GAMMA
+#undef MU
+#undef IOTA
+#undef BETA
+#undef BETA_SD
+#undef POPSIZE
+
+#undef SUSC
+#undef INFD
+#undef RCVD
+#undef CASE
+#undef W
+#undef BIRTHS
+#undef DW
+
+#undef SEASBASIS
Deleted: pkg/src/sir.c
===================================================================
--- pkg/src/sir.c 2009-02-01 14:07:02 UTC (rev 74)
+++ pkg/src/sir.c 2009-02-03 21:22:04 UTC (rev 75)
@@ -1,265 +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 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 BIRTHS (x[stateindex[5]]) // births
-#define dW (x[stateindex[6]]) // white noise process
-
-#define SEASBASIS (covar[covindex[0]]) // first column of seasonality basis in lookup table
-
-// 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; // transition numbers
- double gamma, mu, iota, beta_sd, beta_var, popsize;
- double beta;
-
- // 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;
-
- beta = exp(dot_product(covdim,&SEASBASIS,&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;
-
- GetRNGstate(); // initialize R's random number generator
-
- 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 = &BIRTHS;
- 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]);
-
- PutRNGstate(); // finished with R's random number generator
-
- // 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
- }
-
-}
-
-#undef BIRTHS
-#undef W
-#undef dW
-
-#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;
-
- // untransform the parameters
- gamma = exp(LOGGAMMA);
- mu = exp(LOGMU);
- iota = exp(LOGIOTA);
- popsize = exp(LOGPOPSIZE);
-
- beta = exp(dot_product(covdim,&SEASBASIS,&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
-
-#define SUSC (x1[stateindex[0]]) // number of susceptibles
-#define INFD (x1[stateindex[1]]) // number of infectives
-#define RCVD (x1[stateindex[2]]) // number of recovereds
-#define CASE (x1[stateindex[3]]) // number of cases (accumulated per reporting period)
-#define W (x1[stateindex[4]]) // integrated white noise
-#define BIRTHS (x2[stateindex[5]]) // births
-#define dW (x2[stateindex[6]]) // white noise process
-
-// 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_density (double *f, double *x1, double *x2, double t1, double t2, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
- int covdim, const double *covar)
-{
- int nrate = 6;
- double rate[nrate]; // transition rates
- double *trans; // transition numbers
- double gamma, mu, iota, beta_sd, popsize;
- double beta;
- double dt = t2-t1;
-
- // untransform the parameters
- gamma = exp(LOGGAMMA);
- mu = exp(LOGMU);
- iota = exp(LOGIOTA);
- beta_sd = exp(LOGBETA_SD);
- popsize = exp(LOGPOPSIZE);
-
- beta = exp(dot_product(covdim,&SEASBASIS,&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))) {
- *f = R_NaN;
- return;
- }
-
- // compute the transition rates
- trans = &BIRTHS;
- if (beta_sd > 0.0) { // environmental noise is ON
- double beta_var = beta_sd*beta_sd;
- *f = dgamma(dW,dt/beta_var,beta_var,1);
- } else { // environmental noise is OFF
- *f = 0; // THIS ASSUMES THAT dw = dt !!!
- }
- 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
- *f += dpois(trans[0],rate[0]*dt,1); // births are Poisson
- *f += deulermultinom(2,SUSC,&rate[1],dt,&trans[1],1);
- *f += deulermultinom(2,INFD,&rate[3],dt,&trans[3],1);
- *f += deulermultinom(1,RCVD,&rate[5],dt,&trans[5],1);
-}
-
-#undef GAMMA
-#undef MU
-#undef IOTA
-#undef BETA
-#undef BETA_SD
-#undef POPSIZE
-
-#undef SUSC
-#undef INFD
-#undef RCVD
-#undef CASE
-#undef W
-#undef BIRTHS
-#undef DW
-
-#undef SEASBASIS
More information about the pomp-commits
mailing list