[Pomp-commits] r268 - in pkg: . data inst inst/data-R inst/doc man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 14 18:34:54 CEST 2010
Author: kingaa
Date: 2010-06-14 18:34:53 +0200 (Mon, 14 Jun 2010)
New Revision: 268
Added:
pkg/data/dacca.rda
pkg/inst/data-R/dacca.R
pkg/man/dacca.Rd
pkg/src/cholmodel.c
Modified:
pkg/DESCRIPTION
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/ou2.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/ChangeLog
pkg/inst/NEWS
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
Log:
The cholera transmission model from King, A. A., Ionides, E. L., Pascual, M., and Bouma, M. J. (2008) Inapparent infections and cholera dynamics, Nature 454:877--880 together with data from the Dacca district of historic Bengal are now included.
Data provided courtesy of Menno J. Bouma.
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-05-31 19:08:09 UTC (rev 267)
+++ pkg/DESCRIPTION 2010-06-14 16:34:53 UTC (rev 268)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.29-4
-Date: 2010-05-31
+Version: 0.29-5
+Date: 2010-06-14
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing
Maintainer: Aaron A. King <kingaa at umich.edu>
URL: http://pomp.r-forge.r-project.org
Added: pkg/data/dacca.rda
===================================================================
(Binary files differ)
Property changes on: pkg/data/dacca.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
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/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2010-05-31 19:08:09 UTC (rev 267)
+++ pkg/inst/ChangeLog 2010-06-14 16:34:53 UTC (rev 268)
@@ -1,3 +1,23 @@
+2010-05-31 kingaa
+
+ * [r267] DESCRIPTION, NAMESPACE, R/nlf.R, R/trajectory-pomp.R: -
+ move deSolve, subplex from Depends to Import
+ * [r266] DESCRIPTION, NAMESPACE, R/aaa.R, R/dmeasure-pomp.R,
+ R/dprocess-pomp.R, R/init.state-pomp.R, R/mif-class.R,
+ R/mif-methods.R, R/mif.R, R/particles-mif.R, R/pfilter.R,
+ R/pomp-fun.R, R/pomp-methods.R, R/pomp.R, R/rmeasure-pomp.R,
+ R/rprocess-pomp.R, R/skeleton-pomp.R, R/trajectory-pomp.R: - use
+ the Collate option in package DESCRIPTION
+
+2010-05-26 kingaa
+
+ * [r264] DESCRIPTION, inst/ChangeLog, inst/doc/ecology.bst,
+ inst/doc/fullnat.bst, inst/doc/intro_to_pomp.Rnw,
+ inst/doc/intro_to_pomp.pdf, src/dmeasure.c, src/pomp_fun.c,
+ src/rmeasure.c: - slight improvements to dmeasure & rmeasure
+ methods
+ - change bibliography style in vignettes
+
2010-05-23 kingaa
* [r263] man/mif-class.Rd, man/mif-methods.Rd, man/mif.Rd: - trying
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2010-05-31 19:08:09 UTC (rev 267)
+++ pkg/inst/NEWS 2010-06-14 16:34:53 UTC (rev 268)
@@ -1,3 +1,7 @@
+NEW IN VERSION 0.29-5:
+ o The cholera transmission model from King, A. A., Ionides, E. L., Pascual, M., and Bouma, M. J. (2008) Inapparent infections and cholera dynamics, Nature 454:877--880 together with data from the Dacca district of historic Bengal are now included.
+ See 'help(dacca)' for details.
+
NEW IN VERSION 0.29-1:
This is a major revision.
Added: pkg/inst/data-R/dacca.R
===================================================================
--- pkg/inst/data-R/dacca.R (rev 0)
+++ pkg/inst/data-R/dacca.R 2010-06-14 16:34:53 UTC (rev 268)
@@ -0,0 +1,194 @@
+## pomp object encoding the "SIRS model with seasonal reservoir" of
+## King, A. A., Ionides, E. L., Pascual, M., & Bouma, M. J.
+## Inapparent infections and cholera dynamics.
+## Nature 454:877-880 (2008)
+## Data are cholera deaths and decadal census figures from the Dacca district of Bengal province, 1891-1941.
+##
+## Data courtesy of Menno J. Bouma, London School of Tropical Medicine & Hygiene
+##
+## Native codes are in the package source.
+
+require(pomp)
+
+mle <- c(
+ gamma=20.8, eps=19.1, rho=0,
+ delta=0.02, deltaI=0.06, clin=1,
+ alpha=1, beta.trend=-0.00498,
+ log.beta1=0.747, log.beta2=6.38, log.beta3=-3.436, log.beta4=4.23, log.beta5=3.33, log.beta6=4.549,
+ omega1=0.184, omega2=0.07857, omega3=0.0584, omega4=0.0091719, omega5=0.00020821, omega6=0.01236,
+ sd.beta=3.1267, tau=0.2296,
+ S.0=6.3457, I.0=3.8615, Rs.0=0, R1.0=0.0086226, R2.0=0.009941, R3.0=1.18243e-06,
+ nbasis=6, nrstage=3
+ )
+
+census <- data.frame(
+ year = c(1891L, 1901L, 1911L, 1921L, 1931L, 1941L),
+ census = c(2420656L, 2649522L, 2960402L, 3125967L, 3432577L, 4222142L)
+ )
+
+cholera <- data.frame(
+ time=seq(from=1891+1/12,to=1941,by=1/12),
+ cholera.deaths = c(
+ 2641L, 939L, 905L, 1219L, 368L, 78L, 29L, 12L, 30L, 44L, 270L, 1149L,
+ 633L, 501L, 855L, 1271L, 666L, 101L, 62L, 23L, 20L, 28L, 461L,
+ 892L, 751L, 170L, 253L, 906L, 700L, 98L, 57L, 72L, 471L, 4217L,
+ 5168L, 4747L, 2380L, 852L, 1166L, 2122L, 576L, 60L, 53L, 62L,
+ 241L, 403L, 551L, 739L, 862L, 348L, 490L, 5596L, 1180L, 142L,
+ 41L, 28L, 39L, 748L, 3934L, 3562L, 587L, 311L, 1639L, 1903L,
+ 601L, 110L, 32L, 19L, 82L, 420L, 1014L, 1073L, 416L, 168L, 909L,
+ 1355L, 447L, 59L, 13L, 21L, 43L, 109L, 338L, 470L, 489L, 394L,
+ 483L, 842L, 356L, 29L, 17L, 16L, 57L, 110L, 488L, 1727L, 1253L,
+ 359L, 245L, 549L, 215L, 9L, 7L, 31L, 236L, 279L, 819L, 1728L,
+ 1942L, 1251L, 3521L, 3412L, 290L, 46L, 35L, 14L, 79L, 852L, 2951L,
+ 2656L, 607L, 172L, 325L, 2191L, 584L, 58L, 38L, 8L, 22L, 50L,
+ 380L, 2059L, 938L, 389L, 767L, 1882L, 286L, 94L, 61L, 10L, 106L,
+ 281L, 357L, 1388L, 810L, 306L, 381L, 1308L, 702L, 87L, 9L, 14L,
+ 36L, 46L, 553L, 1302L, 618L, 147L, 414L, 768L, 373L, 39L, 10L,
+ 36L, 151L, 1130L, 3437L, 4041L, 1415L, 207L, 92L, 128L, 147L,
+ 32L, 7L, 59L, 426L, 2644L, 2891L, 4249L, 2291L, 797L, 680L, 1036L,
+ 404L, 41L, 19L, 12L, 10L, 121L, 931L, 2158L, 1886L, 803L, 397L,
+ 613L, 132L, 48L, 17L, 22L, 26L, 34L, 344L, 657L, 117L, 75L, 443L,
+ 972L, 646L, 107L, 18L, 6L, 9L, 5L, 12L, 142L, 133L, 189L, 1715L,
+ 3115L, 1412L, 182L, 50L, 37L, 77L, 475L, 1730L, 1489L, 620L,
+ 190L, 571L, 1558L, 440L, 27L, 7L, 14L, 93L, 1462L, 2467L, 1703L,
+ 1262L, 458L, 453L, 717L, 232L, 26L, 16L, 18L, 9L, 78L, 353L,
+ 897L, 777L, 404L, 799L, 2067L, 613L, 98L, 19L, 26L, 47L, 171L,
+ 767L, 1896L, 887L, 325L, 816L, 1653L, 355L, 85L, 54L, 88L, 609L,
+ 882L, 1363L, 2178L, 580L, 396L, 1493L, 2154L, 683L, 78L, 19L,
+ 10L, 27L, 88L, 1178L, 1862L, 611L, 478L, 2697L, 3395L, 520L,
+ 67L, 41L, 36L, 209L, 559L, 971L, 2144L, 1099L, 494L, 586L, 508L,
+ 269L, 27L, 19L, 21L, 12L, 22L, 333L, 676L, 487L, 262L, 535L,
+ 979L, 170L, 25L, 9L, 19L, 13L, 45L, 229L, 673L, 432L, 107L, 373L,
+ 1126L, 339L, 19L, 11L, 3L, 15L, 101L, 539L, 709L, 200L, 208L,
+ 926L, 1783L, 831L, 103L, 37L, 17L, 33L, 179L, 426L, 795L, 481L,
+ 491L, 773L, 936L, 325L, 101L, 22L, 25L, 24L, 88L, 633L, 513L,
+ 298L, 93L, 687L, 1750L, 356L, 33L, 2L, 18L, 70L, 648L, 2471L,
+ 1270L, 616L, 193L, 706L, 1372L, 668L, 107L, 58L, 21L, 23L, 93L,
+ 318L, 867L, 332L, 118L, 437L, 2233L, 491L, 27L, 7L, 21L, 96L,
+ 360L, 783L, 1492L, 550L, 176L, 633L, 922L, 267L, 91L, 42L, 4L,
+ 10L, 7L, 43L, 377L, 563L, 284L, 298L, 625L, 131L, 35L, 12L, 8L,
+ 9L, 83L, 502L, 551L, 256L, 198L, 664L, 1701L, 425L, 76L, 17L,
+ 9L, 16L, 5L, 141L, 806L, 1603L, 587L, 530L, 771L, 511L, 97L,
+ 35L, 39L, 156L, 1097L, 1233L, 1418L, 1125L, 420L, 1592L, 4169L,
+ 1535L, 371L, 139L, 55L, 85L, 538L, 1676L, 1435L, 804L, 370L,
+ 477L, 394L, 306L, 132L, 84L, 87L, 53L, 391L, 1541L, 1859L, 894L,
+ 326L, 853L, 1891L, 1009L, 131L, 77L, 63L, 66L, 33L, 178L, 1003L,
+ 1051L, 488L, 911L, 1806L, 837L, 280L, 132L, 76L, 381L, 1328L,
+ 2639L, 2164L, 1082L, 326L, 254L, 258L, 119L, 106L, 93L, 29L,
+ 17L, 17L, 17L, 46L, 79L, 135L, 1290L, 2240L, 561L, 116L, 24L,
+ 15L, 33L, 18L, 16L, 38L, 26L, 45L, 151L, 168L, 57L, 32L, 29L,
+ 27L, 20L, 106L, 1522L, 2013L, 434L, 205L, 528L, 634L, 195L, 45L,
+ 33L, 19L, 20L, 46L, 107L, 725L, 572L, 183L, 2199L, 4018L, 428L,
+ 67L, 31L, 8L, 44L, 484L, 1324L, 2054L, 467L, 216L, 673L, 887L,
+ 353L, 73L, 46L, 15L, 20L, 27L, 25L, 38L, 158L, 312L, 1226L, 1021L,
+ 222L, 90L, 31L, 93L, 368L, 657L, 2208L, 2178L, 702L, 157L, 317L,
+ 146L, 63L, 27L, 22L, 23L, 28L, 225L, 483L, 319L, 120L, 59L, 274L,
+ 282L, 155L, 31L, 16L, 15L, 12L, 14L, 14L, 42L
+ )
+ )
+
+covar.dt <- 0.01
+nbasis <- as.integer(mle["nbasis"])
+nrstage <- as.integer(mle["nrstage"])
+
+t0 <- with(cholera,2*time[1]-time[2])
+tcovar <- seq(from=t0,to=max(cholera$time)+2/12,by=covar.dt)
+covartable <- data.frame(
+ time=tcovar,
+ seas=periodic.bspline.basis(tcovar-1/12,nbasis=nbasis,degree=3,period=1),
+ pop=predict(smooth.spline(x=census$year,y=census$census),x=tcovar)$y,
+ dpopdt=predict(smooth.spline(x=census$year,y=census$census),x=tcovar,deriv=1)$y,
+ trend=tcovar-mean(tcovar)
+ )
+
+dacca <- list()
+
+pomp(
+ data=cholera,
+ times='time',
+ t0=t0,
+ nrstage = nrstage,
+ rprocess = euler.sim(
+ step.fun = "_cholmodel_one",
+ PACKAGE="pomp",
+ delta.t=1/240
+ ),
+ dmeasure = "_cholmodel_norm_dmeasure",
+ rmeasure="_cholmodel_norm_rmeasure",
+ PACKAGE = "pomp",
+ covar=covartable,
+ tcovar='time',
+ zeronames = c("M","count"),
+ obsnames = "cholera.deaths",
+ statenames = c("S","I","Rs","R1","M","W","count"),
+ paramnames = c("tau","gamma","eps","delta","deltaI",
+ "omega1","sd.beta","beta.trend","log.beta1",
+ "alpha","rho","clin","nbasis","nrstage"),
+ covarnames = c("pop","dpopdt","seas.1","trend"),
+ ## log.trans=c( # parameters to log transform
+ ## "gamma","eps","rho","delta","deltaI","alpha",
+ ## "tau","sd.beta",
+ ## paste("omega",seq(length=nbasis),sep=''),
+ ## "S.0","I.0","Rs.0","R1.0","R2.0","R3.0"
+ ## ),
+ ## logit.trans="clin", # parameters to logit transform
+ ## transform=function (params, log.trans, logit.trans, ...) {
+ ## logit <- function(p){log(p/(1-p))} # (0,1) -> (-inf,inf)
+ ## params[logit.trans] <- logit(params[logit.trans])
+ ## params[log.trans] <- log(params[log.trans])
+ ## params
+ ## },
+ ## inv.transform=function (params, log.trans, logit.trans, ...) {
+ ## expit <- function(r){1/(1+exp(-r))} # (-inf,inf) -> (0,1)
+ ## params[logit.trans] <- expit(params[logit.trans])
+ ## params[log.trans] <- exp(params[log.trans])
+ ## params
+ ## },
+ initializer = function (params, t0, covars, all.state.names, comp.names, nrstage, ...) {
+ all.state.names <- c("S","I","Rs","R1","R2","R3","M","W","count")
+ comp.names <- c("S","I","Rs",paste("R",1:nrstage,sep=''))
+ comp.ic.names <- paste(comp.names,"0",sep='.')
+ states <- numeric(length(all.state.names))
+ names(states) <- all.state.names
+ ## translate fractions into initial conditions
+ frac <- exp(params[comp.ic.names])
+ states[comp.names] <- round(covars['pop']*frac/sum(frac))
+ states
+ }
+ ) -> dacca$po
+
+## Parameter transformations
+## Parameters are fit in the transformed space.
+## Positive parameters are log transformed, probabilities are logit transformed.
+## Initial conditions are parameterized in terms of log-transformed fractions in each compartment.
+dacca$transform <- function (params, dir = c("forward","inverse")) {
+ dir <- match.arg(dir)
+ r <- length(dim(params))
+ nm <- if (r>0) rownames(params) else names(params)
+ log.trans=c( # parameters to log transform
+ "gamma","eps","rho","delta","deltaI","alpha",
+ "tau","sd.beta",
+ grep("omega.$",nm,val=TRUE),
+ "S.0","I.0","Rs.0","R1.0","R2.0","R3.0"
+ )
+ logit.trans="clin"
+ logit <- function(p){log(p/(1-p))} # (0,1) -> (-inf,inf)
+ expit <- function(r){1/(1+exp(-r))} # (-inf,inf) -> (0,1)
+ par.trans.vec <- function (x) {
+ x[logit.trans] <- logit(x[logit.trans])
+ x[log.trans] <- log(x[log.trans])
+ x
+ }
+ par.untrans.vec <- function (x) {
+ x[logit.trans] <- expit(x[logit.trans])
+ x[log.trans] <- exp(x[log.trans])
+ x
+ }
+ switch(
+ dir,
+ forward=if (r > 1) apply(params,2:r,par.trans.vec) else par.trans.vec(params),
+ inverse=if (r > 1) apply(params,2:r,par.untrans.vec) else par.untrans.vec(params)
+ )
+}
+
+coef(dacca$po) <- dacca$transform(mle,dir="forward")
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Added: pkg/man/dacca.Rd
===================================================================
--- pkg/man/dacca.Rd (rev 0)
+++ pkg/man/dacca.Rd 2010-06-14 16:34:53 UTC (rev 268)
@@ -0,0 +1,43 @@
+\name{dacca}
+\alias{dacca}
+\docType{data}
+\title{Model of cholera transmission for historic Bengal.}
+\description{
+ \code{dacca} contains census and cholera mortality data from the Dacca district of the former British province of Bengal over the years 1891 to 1940 together with a stochastic differential equation transmission model.
+ The model is that of King et al. (2008).
+ It also has the MLE for the SIRS model with seasonal reservoir.
+
+ Data are provided courtesy of Dr. Menno J. Bouma, London School of Tropical Medicine and Hygiene.
+}
+\usage{
+data(dacca)
+}
+\details{
+ \code{dacca} is a list consisting of
+ \describe{
+ \item{po}{
+ A \code{pomp} object containing the model, data, and MLE parameters.
+ }
+ \item{transform}{
+ A function for transforming the parameters.
+ }
+ }
+}
+\examples{
+data(dacca)
+plot(dacca$po)
+theta <- dacca$transform(coef(dacca$po),dir="inverse") #MLEs
+plot(simulate(dacca$po))
+theta["eps"] <- 1
+plot(simulate(dacca$po,params=dacca$transform(theta,dir="forward")))
+}
+\references{
+ King, A. A., Ionides, E. L., Pascual, M., and Bouma, M. J.
+ Inapparent infections and cholera dynamics.
+ Nature 454:877-880 (2008)
+}
+\seealso{
+ \code{\link{euler.sir}},
+ \code{\link{pomp}}
+}
+\keyword{datasets}
Added: pkg/src/cholmodel.c
===================================================================
--- pkg/src/cholmodel.c (rev 0)
+++ pkg/src/cholmodel.c 2010-06-14 16:34:53 UTC (rev 268)
@@ -0,0 +1,204 @@
+// -*- C++ -*-
+
+#include <R.h>
+#include <Rmath.h>
+#include "pomp.h"
+
+double expit (double x) {return 1.0/(1.0+exp(-x));}
+
+#define TAU (p[parindex[0]])
+#define GAMMA (p[parindex[1]])
+#define EPS (p[parindex[2]])
+#define DELTA (p[parindex[3]])
+#define DELTA_I (p[parindex[4]])
+#define OMEGA (p[parindex[5]])
+#define SD_BETA (p[parindex[6]])
+#define BETATREND (p[parindex[7]])
+#define LOGBETA (p[parindex[8]])
+#define ALPHA (p[parindex[9]])
+#define RHO (p[parindex[10]])
+#define CLIN (p[parindex[11]])
+#define NBASIS (p[parindex[12]])
+#define NRSTAGE (p[parindex[13]])
+
+#define SUSCEP (x[stateindex[0]])
+#define INFECT (x[stateindex[1]])
+#define RSHORT (x[stateindex[2]])
+#define RLONG (x[stateindex[3]])
+#define DEATHS (x[stateindex[4]])
+#define NOISE (x[stateindex[5]])
+#define COUNT (x[stateindex[6]])
+
+#define POP (covar[covindex[0]])
+#define DPOPDT (covar[covindex[1]])
+#define SEASBASIS (covar[covindex[2]])
+#define TREND (covar[covindex[3]])
+
+#define DATADEATHS (y[obsindex[0]])
+
+void _cholmodel_norm_rmeasure (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t)
+{
+ double v, tol = 1.0e-18;
+ v = DEATHS*exp(TAU);
+ if ((COUNT > 0) || (!(R_FINITE(v)))) {
+ DATADEATHS = R_NaReal;
+ } else {
+ DATADEATHS = rnorm(DEATHS,v+tol);
+ }
+}
+
+void _cholmodel_norm_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)
+{
+ double v, tol = 1.0e-18;
+ v = DEATHS*exp(TAU);
+ if ((COUNT>0.0) || (!(R_FINITE(v)))) {
+ *lik = tol;
+ } else {
+ *lik = dnorm(DATADEATHS,DEATHS,v+tol,0)+tol;
+ }
+ if (give_log) *lik = log(*lik);
+}
+
+#undef DATADEATHS
+
+// two-path SIRS cholera model using SDEs
+// exponent (alpha) on I/n
+// only "severe" infections are infectious
+// truncation is not used
+// instead, particles with negative states are killed
+void _cholmodel_one (double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar,
+ double t, double dt)
+{ // implementation of the SIRS cholera model
+ int nrstage = (int) NRSTAGE;
+ int nbasis = (int) NBASIS;
+ double births;
+ double infections;
+ double sdeaths;
+ double ideaths;
+ double rsdeaths;
+ double rldeaths[nrstage];
+ double disease;
+ double wanings;
+ double passages[nrstage+1];
+ double effI;
+ double beta;
+ double omega;
+ double gamma;
+ double rho;
+ double clin;
+ double eps;
+ double delta;
+ double deltaI;
+ double sd_beta;
+ double alpha;
+ double dw;
+ const double *pt;
+ int j;
+
+ if (!(R_FINITE(SUSCEP))) return;
+ if (!(R_FINITE(INFECT))) return;
+ if (!(R_FINITE(RSHORT))) return;
+ for (pt = &RLONG, j = 0; j < nrstage; j++) {
+ if (!(R_FINITE(pt[j]))) return;
+ }
+
+ gamma = exp(GAMMA);
+ eps = exp(EPS)*NRSTAGE;
+ rho = exp(RHO);
+ delta = exp(DELTA);
+ deltaI = exp(DELTA_I);
+ clin = expit(CLIN);
+ sd_beta = exp(SD_BETA);
+ alpha = exp(ALPHA);
+
+ if (!(R_FINITE(gamma))) return;
+ if (!(R_FINITE(eps))) return;
+ if (!(R_FINITE(rho))) return;
+ if (!(R_FINITE(delta))) return;
+ if (!(R_FINITE(deltaI))) return;
+ if (!(R_FINITE(clin))) return;
+ if (!(R_FINITE(sd_beta))) return;
+ if (!(R_FINITE(alpha))) return;
+ if (!(R_FINITE(BETATREND))) return;
+
+ beta = exp(dot_product(nbasis,&SEASBASIS,&LOGBETA)+BETATREND*TREND);
+ omega = exp(dot_product(nbasis,&SEASBASIS,&OMEGA));
+
+ if (!(R_FINITE(beta))) return;
+ if (!(R_FINITE(omega))) return;
+
+ dw = rnorm(0,sqrt(dt)); // white noise
+
+ effI = pow(INFECT/POP,alpha);
+ births = DPOPDT + delta*POP; // births
+
+ passages[0] = gamma*INFECT; // recovery
+ ideaths = delta*INFECT; // natural I deaths
+ disease = deltaI*INFECT; // disease death
+ rsdeaths = delta*RSHORT; // natural Rs deaths
+ wanings = rho*RSHORT; // loss of immunity
+
+ for (pt = &RLONG, j = 0; j < nrstage; j++) {
+ rldeaths[j] = delta*pt[j]; // natural R deaths
+ passages[j+1] = eps*pt[j]; // passage to the next immunity class
+ }
+
+ infections = (omega+(beta+sd_beta*dw/dt)*effI)*SUSCEP; // infection
+ sdeaths = delta*SUSCEP; // natural S deaths
+ if (infections > 0.0) {
+ if ((infections+sdeaths)*dt > SUSCEP) { // too many leaving S class
+ COUNT += 1.0e10;
+ return;
+ }
+ } else {
+ if ((-clin*infections+disease+ideaths+passages[0])*dt > INFECT) { // too many leaving I class
+ COUNT += 1.0e5;
+ return;
+ }
+ if ((-(1-clin)*infections+rsdeaths+wanings)*dt > RSHORT) { // too many leaving Rs class
+ COUNT += 1;
+ return;
+ }
+ }
+
+ SUSCEP += (births - infections - sdeaths + passages[nrstage] + wanings)*dt;
+ INFECT += (clin*infections - disease - ideaths - passages[0])*dt;
+ RSHORT += ((1-clin)*infections - rsdeaths - wanings)*dt;
+ for (j = 0; j < nrstage; j++)
+ (&RLONG)[j] += (passages[j] - passages[j+1] - rldeaths[j])*dt;
+ DEATHS += disease*dt; // cumulative deaths due to disease
+ NOISE += dw;
+
+}
+
+#undef GAMMA
+#undef EPS
+#undef DELTA
+#undef DELTA_I
+#undef OMEGA
+#undef SD_BETA
+#undef BETATREND
+#undef LOGBETA
+#undef ALPHA
+#undef RHO
+#undef CLIN
+#undef NBASIS
+#undef NRSTAGE
+
+#undef SUSCEP
+#undef INFECT
+#undef RSHORT
+#undef RLONG
+#undef DEATHS
+#undef NOISE
+#undef COUNT
+
+#undef POP
+#undef DPOPDT
+#undef SEASBASIS
Property changes on: pkg/src/cholmodel.c
___________________________________________________________________
Added: svn:eol-style
+ native
More information about the pomp-commits
mailing list