[Pomp-commits] r485 - in pkg: data inst/examples src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 19 21:27:02 CEST 2011
Author: kingaa
Date: 2011-05-19 21:27:01 +0200 (Thu, 19 May 2011)
New Revision: 485
Modified:
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/inst/examples/sir.c
pkg/src/sir.c
Log:
- change the measurement model for the SIR example
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/examples/sir.c
===================================================================
--- pkg/inst/examples/sir.c 2011-05-19 19:22:36 UTC (rev 484)
+++ pkg/inst/examples/sir.c 2011-05-19 19:27:01 UTC (rev 485)
@@ -12,6 +12,12 @@
return log(x/(1-x));
}
+// some macros to make the code easier to read.
+// note how variables and parameters use the indices.
+// the integer vectors parindex, stateindex, obsindex
+// are computed by pomp by matching the names of input vectors
+// against parnames, statenames, and obsnames, respectively.
+
#define LOGGAMMA (p[parindex[0]]) // recovery rate
#define LOGMU (p[parindex[1]]) // baseline birth and death rate
#define LOGIOTA (p[parindex[2]]) // import rate
@@ -23,32 +29,36 @@
#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 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 REPORTS (y[obsindex[0]])
+#define REPORTS (y[obsindex[0]]) // the data
+// the measurement model.
+// note that dmeasure and rmeasure are mirrors of each other.
+
void binomial_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 = dbinom(REPORTS,CASE,exp(LOGRHO),give_log);
+ *lik = dbinom(REPORTS,nearbyint(CASE),exp(LOGRHO),give_log);
}
void binomial_rmeasure (double *y, double *x, double *p,
int *obsindex, int *stateindex, int *parindex, int *covindex,
int ncovars, double *covars, double t) {
- REPORTS = rbinom(CASE,exp(LOGRHO));
+ REPORTS = rbinom(nearbyint(CASE),exp(LOGRHO));
}
#undef REPORTS
-// SIR model with Euler multinomial step
-// forced transmission (basis functions passed as covariates)
-// constant population size as a parameter
-// environmental stochasticity on transmission
+// the process model:
+// an SIR model with Euler-multinomial step,
+// transmission is seasonal, implemented using B-spline basis functions passed as covariates.
+// the population size is constant and specified by a parameter.
+// there is environmental stochasticity in transmission (dW).
void sir_euler_simulator (double *x, const double *p,
const int *stateindex, const int *parindex, const int *covindex,
int covdim, const double *covar,
@@ -58,8 +68,8 @@
double rate[nrate]; // transition rates
double trans[nrate]; // transition numbers
double gamma, mu, iota, beta_sd, beta_var, popsize;
- double beta;
- double dW;
+ double beta; // transmission rate
+ double dW; // white noise increment
int nseas = (int) NBASIS; // number of seasonal basis functions
int deg = (int) DEGREE; // degree of seasonal basis functions
double seasonality[nseas];
@@ -75,15 +85,20 @@
popsize = exp(LOGPOPSIZE);
beta_var = beta_sd*beta_sd;
+ // to evaluate the basis functions and compute the transmission rate, use some of
+ // pomp's built-in C-level facilities:
eval_basis = (void (*)(double,double,int,int,double*)) R_GetCCallable("pomp","periodic_bspline_basis_eval");
+ dot = (double (*)(int,const double *,const double *)) R_GetCCallable("pomp","dot_product");
+ // pomp's C-level eulermultinomial simulator
reulmult = (void (*)(int,double,double*,double,double*)) R_GetCCallable("pomp","reulermultinom");
- dot = (double (*)(int,const double *,const double *)) R_GetCCallable("pomp","dot_product");
+ // compute transmission rate from seasonality
if (nseas <= 0) return;
(*eval_basis)(t,PERIOD,deg,nseas,&seasonality[0]); // evaluate the periodic B-spline basis
- beta = exp((*dot)(nseas,&seasonality[0],&LOGBETA));
+ beta = exp((*dot)(nseas,&seasonality[0],&LOGBETA)); // compute the transmission rate
// test to make sure the parameters and state variable values are sane
+ // if not, return without doing any work
if (!(R_FINITE(beta)) ||
!(R_FINITE(gamma)) ||
!(R_FINITE(mu)) ||
@@ -97,6 +112,7 @@
!(R_FINITE(W)))
return;
+ // compute the environmental stochasticity
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;
@@ -114,9 +130,9 @@
// compute the transition numbers
trans[0] = rpois(rate[0]*dt); // births are Poisson
- (*reulmult)(2,SUSC,&rate[1],dt,&trans[1]);
- (*reulmult)(2,INFD,&rate[3],dt,&trans[3]);
- (*reulmult)(1,RCVD,&rate[5],dt,&trans[5]);
+ (*reulmult)(2,SUSC,&rate[1],dt,&trans[1]); // euler-multinomial exits from S class
+ (*reulmult)(2,INFD,&rate[3],dt,&trans[3]); // euler-multinomial exits from I class
+ (*reulmult)(1,RCVD,&rate[5],dt,&trans[5]); // euler-multinomial exits from R class
// balance the equations
SUSC += trans[0]-trans[1]-trans[2];
@@ -124,7 +140,7 @@
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
+ W += (dW-dt)/beta_sd; // increment mean zero, variance = dt
}
}
@@ -155,9 +171,12 @@
iota = exp(LOGIOTA);
popsize = exp(LOGPOPSIZE);
+ // to evaluate the basis functions and compute the transmission rate, use some of
+ // pomp's built-in C-level facilities:
eval_basis = (void (*)(double,double,int,int,double*)) R_GetCCallable("pomp","periodic_bspline_basis_eval");
dot = (double (*)(int,const double *,const double *)) R_GetCCallable("pomp","dot_product");
+ // compute transmission rate from seasonality
if (nseas <= 0) return;
(*eval_basis)(t,PERIOD,deg,nseas,&seasonality[0]); // evaluate the periodic B-spline basis
beta = exp((*dot)(nseas,&seasonality[0],&LOGBETA));
@@ -178,11 +197,14 @@
term[4] = rate[4]*INFD;
term[5] = rate[5]*RCVD;
- // balance the equations
+ // assemble the differential 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
+ // note that we cannot do the same accumulation trick using 'zeronames'
+ // that we do in the stochastic simulator
+ // instead, we'll try to get roughly the number of new infections in a week using:
+ DCDT = DIDT*gamma/52; // this assumes that the reporting rate is weekly!
}
Modified: pkg/src/sir.c
===================================================================
--- pkg/src/sir.c 2011-05-19 19:22:36 UTC (rev 484)
+++ pkg/src/sir.c 2011-05-19 19:27:01 UTC (rev 485)
@@ -46,13 +46,35 @@
void _sir_binom_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 = dbinom(REPORTS,nearbyint(CASE),exp(LOGRHO),give_log);
+ double size, prob;
+ double mean, sd;
+ double f;
+ size = CASE;
+ prob = exp(LOGRHO);
+ mean = size*prob;
+ sd = sqrt(size*prob*(1-prob));
+ if (REPORTS > 0) {
+ f = pnorm(REPORTS+0.5,mean,sd,1,0)-pnorm(REPORTS-0.5,mean,sd,1,0);
+ } else {
+ f = pnorm(REPORTS+0.5,mean,sd,1,0);
+ }
+ *lik = (give_log) ? log(f) : f;
+ // *lik = dbinom(REPORTS,nearbyint(CASE),exp(LOGRHO),give_log);
}
void _sir_binom_rmeasure (double *y, double *x, double *p,
int *obsindex, int *stateindex, int *parindex, int *covindex,
int ncovars, double *covars, double t) {
- REPORTS = rbinom(nearbyint(CASE),exp(LOGRHO));
+ double size, prob;
+ double mean, sd;
+ double rep;
+ size = CASE;
+ prob = exp(LOGRHO);
+ mean = size*prob;
+ sd = sqrt(size*prob*(1-prob));
+ rep = nearbyint(rnorm(mean,sd));
+ REPORTS = (rep > 0) ? rep : 0;
+ // REPORTS = rbinom(nearbyint(CASE),exp(LOGRHO));
}
#undef REPORTS
More information about the pomp-commits
mailing list