[Pomp-commits] r900 - in pkg/pomp: R demo inst/examples man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 20 22:51:41 CET 2014
Author: kingaa
Date: 2014-03-20 22:51:41 +0100 (Thu, 20 Mar 2014)
New Revision: 900
Removed:
pkg/pomp/inst/examples/sir.c
Modified:
pkg/pomp/R/mif.R
pkg/pomp/demo/sir.R
pkg/pomp/inst/examples/bbs.R
pkg/pomp/inst/examples/euler.sir.R
pkg/pomp/man/mif.Rd
pkg/pomp/tests/bbs.Rout.save
pkg/pomp/tests/dacca.Rout.save
pkg/pomp/tests/ou2-mif.Rout.save
pkg/pomp/tests/pfilter.Rout.save
Log:
- removed inst/examples/sir.c
- edited demo/sir.R and pompBuilder section in Advanced Topics vignette to agree with one another
- small change to bbs example
- euler.sir example now uses negative binomial measurement error
- when 'mif' is called on a 'mif' object, pars are no longer copied but rather taken from 'rw.sd'
- clean up 'dacca' test
- update tests
Modified: pkg/pomp/R/mif.R
===================================================================
--- pkg/pomp/R/mif.R 2014-03-20 21:51:17 UTC (rev 899)
+++ pkg/pomp/R/mif.R 2014-03-20 21:51:41 UTC (rev 900)
@@ -491,7 +491,7 @@
signature=signature(object="mif"),
function (object, Nmif,
start,
- pars, ivps,
+ ivps,
particles, rw.sd,
Np, ic.lag, var.factor,
cooling.type, cooling.fraction,
@@ -502,7 +502,6 @@
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(particles)) particles <- object at particles
if (missing(rw.sd)) rw.sd <- object at random.walk.sd
@@ -520,7 +519,6 @@
object=as(object,"pomp"),
Nmif=Nmif,
start=start,
- pars=pars,
ivps=ivps,
particles=particles,
rw.sd=rw.sd,
Modified: pkg/pomp/demo/sir.R
===================================================================
--- pkg/pomp/demo/sir.R 2014-03-20 21:51:17 UTC (rev 899)
+++ pkg/pomp/demo/sir.R 2014-03-20 21:51:41 UTC (rev 900)
@@ -1,15 +1,18 @@
require(pomp)
## negative binomial measurement model
-dmeas <- "
- double prob = theta/(theta+rho*incid);
- lik = dnbinom(cases,theta,prob,give_log);
-"
-rmeas <- "
- double prob = theta/(theta+rho*incid);
- cases = rnbinom(theta,prob);
-"
+## E[cases|incid] = rho*incid
+## Var[cases|incid] = rho*incid*(1+rho*incid/theta)
+rmeas <- '
+ cases = rnbinom_mu(theta,rho*incid);
+'
+
+dmeas <- '
+ lik = dnbinom_mu(cases,theta,rho*incid,give_log);
+'
+
## SIR process model with extra-demographic stochasticity
+## and seasonal transmission
step.fn <- '
int nrate = 6;
double rate[nrate]; // transition rates
@@ -18,6 +21,7 @@
double dW; // white noise increment
int k;
+ // seasonality in transmission
beta = beta1*seas1+beta2*seas2+beta3*seas3;
// compute the environmental stochasticity
@@ -41,9 +45,10 @@
S += trans[0]-trans[1]-trans[2];
I += trans[1]-trans[3]-trans[4];
R += trans[3]-trans[5];
- incid += trans[3]; // cases are cumulative recoveries
+ incid += trans[3]; // incidence is 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
@@ -77,9 +82,12 @@
Dincid = term[3]; // accumulate the new I->R transitions
DW = 0;
'
+
## parameter transformations
+## note we use barycentric coordinates for the initial conditions
+## the success of this depends on S0, I0, R0 being in
+## adjacent memory locations, in that order
partrans <- "
- double sum;
Tgamma = exp(gamma);
Tmu = exp(mu);
Tiota = exp(iota);
@@ -89,16 +97,10 @@
Tbeta_sd = exp(beta_sd);
Trho = expit(rho);
Ttheta = exp(theta);
- TS_0 = exp(S_0);
- TI_0 = exp(I_0);
- TR_0 = exp(R_0);
- sum = TS_0+TI_0+TR_0;
- TS_0 /= sum;
- TI_0 /= sum;
- TR_0 /= sum;
+ from_log_barycentric(&TS_0,&S_0,3);
"
+
paruntrans <- "
- double sum;
Tgamma = log(gamma);
Tmu = log(mu);
Tiota = log(iota);
@@ -108,10 +110,7 @@
Tbeta_sd = log(beta_sd);
Trho = logit(rho);
Ttheta = log(theta);
- sum = S_0+I_0+R_0;
- TS_0 = log(S_0/sum);
- TI_0 = log(I_0/sum);
- TR_0 = log(R_0/sum);
+ to_log_barycentric(&TS_0,&S_0,3);
"
data(LondonYorke)
@@ -170,7 +169,7 @@
beta1=120,beta2=140,beta3=100,
beta.sd=0.01,
popsize=5e6,
- rho=0.2,theta=0.001,
+ rho=0.1,theta=1,
S.0=0.22,I.0=0.0018,R.0=0.78
)
@@ -189,18 +188,3 @@
print(toc-tic)
plot(incid~time,data=x,col=as.factor(x$sim),pch=16)
-
-coef(po) <- coef(
- traj.match(
- pomp(
- window(po,end=1930)
- ## window(po,end=1930),
- ## measurement.model=cases~norm(mean=rho*incid,sd=100)
- ),
- est=c("S.0","I.0","R.0"),
- transform=TRUE
- )
- )
-
-pf <- pfilter(po,Np=1000,max.fail=100)
-print(round(logLik(pf),1))
Modified: pkg/pomp/inst/examples/bbs.R
===================================================================
--- pkg/pomp/inst/examples/bbs.R 2014-03-20 21:51:17 UTC (rev 899)
+++ pkg/pomp/inst/examples/bbs.R 2014-03-20 21:51:41 UTC (rev 900)
@@ -37,13 +37,13 @@
),
skeleton.type="vectorfield",
skeleton="_sir_ODE",
- measurement.model=reports~nbinom(mu=rho*cases,size=1/sigma),
+ measurement.model=reports~nbinom(mu=rho*cases,size=1/sigma^2),
PACKAGE="pomp",
obsnames = c("reports"),
statenames=c("S","I","R","cases","W"),
paramnames=c(
"gamma","mu","iota",
- "beta","beta.sd","pop","rho",
+ "beta","beta.sd","pop","rho","sigma",
"S.0","I.0","R.0"
),
zeronames=c("cases"),
Modified: pkg/pomp/inst/examples/euler.sir.R
===================================================================
--- pkg/pomp/inst/examples/euler.sir.R 2014-03-20 21:51:17 UTC (rev 899)
+++ pkg/pomp/inst/examples/euler.sir.R 2014-03-20 21:51:41 UTC (rev 900)
@@ -220,7 +220,7 @@
beta1=400,beta2=480,beta3=320,
beta.sd=1e-3,
pop=2.1e6,
- rho=0.6,
+ rho=0.6,overdisp=1,
S.0=26/400,I.0=0.001,R.0=1-26/400
),
rprocess=euler.sim(
@@ -237,7 +237,7 @@
statenames=c("S","I","R","cases","W"),
paramnames=c(
"gamma","mu","iota",
- "beta1","beta.sd","pop","rho",
+ "beta1","beta.sd","pop","rho","overdisp",
"S.0","I.0","R.0"
),
zeronames=c("cases"),
Deleted: pkg/pomp/inst/examples/sir.c
===================================================================
--- pkg/pomp/inst/examples/sir.c 2014-03-20 21:51:17 UTC (rev 899)
+++ pkg/pomp/inst/examples/sir.c 2014-03-20 21:51:41 UTC (rev 900)
@@ -1,165 +0,0 @@
-// dear emacs, please treat this as -*- C++ -*-
-
-#include <R.h>
-#include <Rmath.h>
-#include <R_ext/Rdynload.h>
-#include <pomp.h>
-
-// SIR example as described in the "Advanced Topics in pomp" vignette.
-// for a demonstration of how to compile, link, and run this example,
-// do 'demo("sir",package="pomp")' at the R prompt.
-
-// 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 GAMMA (p[parindex[0]]) // recovery rate
-#define MU (p[parindex[1]]) // baseline birth and death rate
-#define IOTA (p[parindex[2]]) // import rate
-#define BETA (&p[parindex[3]]) // transmission rate
-#define BETA_SD (p[parindex[4]]) // environmental stochasticity SD in transmission rate
-#define POPSIZE (p[parindex[5]]) // population size
-#define RHO (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 REPORTS (y[obsindex[0]]) // the data
-
-#define DSDT (f[stateindex[0]])
-#define DIDT (f[stateindex[1]])
-#define DRDT (f[stateindex[2]])
-#define DCDT (f[stateindex[3]])
-#define DWDT (f[stateindex[4]])
-
-// 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,nearbyint(CASE),RHO,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(nearbyint(CASE),RHO);
-}
-
-// 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,
- double t, double dt)
-{
- int nrate = 6;
- double rate[nrate]; // transition rates
- double trans[nrate]; // transition numbers
- double beta; // transmission rate
- double dW; // white noise increment
- int nbasis = (int) NBASIS; // number of seasonal basis functions
- int deg = (int) DEGREE; // degree of seasonal basis functions
- int k;
- double seasonality[nbasis];
- void (*eval_basis)(double,double,int,int,double*);
- void (*reulmult)(int,double,double*,double,double*);
-
- // 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");
- // pomp's C-level eulermultinomial simulator
- reulmult = (void (*)(int,double,double*,double,double*)) R_GetCCallable("pomp","reulermultinom");
-
- // compute transmission rate from seasonality
- if (nbasis <= 0) return;
- eval_basis(t,PERIOD,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
- for (k = 0, beta = 0; k < nbasis; k++)
- beta += BETA[k]*seasonality[k];
-
- // compute the environmental stochasticity
- dW = rgammawn(BETA_SD,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
- (*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];
- 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; // increment has mean = 0, variance = dt
-
-}
-
-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 beta;
- int nbasis = (int) NBASIS; // number of seasonal basis functions
- int deg = (int) DEGREE; // degree of seasonal basis functions
- int k;
- double seasonality[nbasis];
- void (*eval_basis)(double,double,int,int,double*);
-
- // 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");
-
- // compute transmission rate from seasonality
- if (nbasis <= 0) return;
- eval_basis(t,PERIOD,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
- for (k = 0, beta = 0; k < nbasis; k++)
- beta += BETA[k]*seasonality[k];
-
- // 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;
-
- // 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]; // accumulate the new I->R transitions
- DWDT = 0;
-
-}
Modified: pkg/pomp/man/mif.Rd
===================================================================
--- pkg/pomp/man/mif.Rd 2014-03-20 21:51:17 UTC (rev 899)
+++ pkg/pomp/man/mif.Rd 2014-03-20 21:51:41 UTC (rev 900)
@@ -23,7 +23,7 @@
tol = 1e-17, max.fail = Inf,
verbose = getOption("verbose"), transform = FALSE, \dots)
\S4method{mif}{pfilterd.pomp}(object, Nmif = 1, Np, tol, \dots)
-\S4method{mif}{mif}(object, Nmif, start, pars, ivps,
+\S4method{mif}{mif}(object, Nmif, start, ivps,
particles, rw.sd, Np, ic.lag, var.factor,
cooling.type, cooling.fraction,
method, tol, transform, \dots)
Modified: pkg/pomp/tests/bbs.Rout.save
===================================================================
--- pkg/pomp/tests/bbs.Rout.save 2014-03-20 21:51:17 UTC (rev 899)
+++ pkg/pomp/tests/bbs.Rout.save 2014-03-20 21:51:41 UTC (rev 900)
@@ -18,6 +18,7 @@
> library(pomp)
Loading required package: mvtnorm
Loading required package: subplex
+Loading required package: nloptr
Loading required package: deSolve
>
> pompExample(bbs)
@@ -51,17 +52,17 @@
Walker's alias method used: results are different from R < 2.2.0
> signif(coef(fit1),3)
gamma mu iota beta beta.sd pop rho sigma
- 0.333 0.000 0.000 3.630 0.000 1400.000 0.900 2.220
+ 0.333 0.000 0.000 4.300 0.000 1400.000 0.900 2.050
S.0 I.0 R.0
0.999 0.001 0.000
>
> fit2 <- traj.match(bbs,est=c("beta","sigma"),transform=TRUE)
> signif(coef(fit2),3)
gamma mu iota beta beta.sd pop rho sigma
- 0.333 0.000 0.000 2.020 0.000 1400.000 0.900 0.426
+ 0.333 0.000 0.000 2.090 0.000 1400.000 0.900 0.474
S.0 I.0 R.0
0.999 0.001 0.000
>
> proc.time()
user system elapsed
- 2.784 0.068 2.880
+ 2.908 0.064 3.003
Modified: pkg/pomp/tests/dacca.Rout.save
===================================================================
--- pkg/pomp/tests/dacca.Rout.save 2014-03-20 21:51:17 UTC (rev 899)
+++ pkg/pomp/tests/dacca.Rout.save 2014-03-20 21:51:41 UTC (rev 900)
@@ -64,6 +64,8 @@
+ r
+ }
+
++ op <- options(warn=-1)
++
+ set.seed(7777+7)
+ params.tricky <- dacca.rprior(dacca.hyperparams)
+ m7 <- mif(
@@ -99,6 +101,8 @@
+ transform=TRUE
+ )
+
++ options(op)
++
+ }
Loading required package: mvtnorm
Loading required package: subplex
@@ -106,16 +110,7 @@
Loading required package: deSolve
newly created pomp object(s):
dacca
-Warning messages:
-1: 11 filtering failures occurred in 'pfilter'
-2: 4 filtering failures occurred in 'pfilter'
-3: 1 filtering failure occurred in 'pfilter'
-4: 16 filtering failures occurred in 'pfilter'
-5: 2 filtering failures occurred in 'pfilter'
-6: 2 filtering failures occurred in 'pfilter'
-7: 1 filtering failure occurred in 'pfilter'
-8: 7 filtering failures occurred in 'pfilter'
>
> proc.time()
user system elapsed
- 11.092 0.056 11.206
+ 11.068 0.104 11.227
Modified: pkg/pomp/tests/ou2-mif.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-mif.Rout.save 2014-03-20 21:51:17 UTC (rev 899)
+++ pkg/pomp/tests/ou2-mif.Rout.save 2014-03-20 21:51:41 UTC (rev 900)
@@ -229,6 +229,7 @@
+ }
Loading required package: mvtnorm
Loading required package: subplex
+Loading required package: nloptr
Loading required package: deSolve
newly created pomp object(s):
ou2
@@ -272,10 +273,9 @@
pfilter timestep 90 of 100 finished
pfilter timestep 95 of 100 finished
pfilter timestep 100 of 100 finished
-Warning messages:
-1: mif warning: the variables 'alpha.3', 'x2.0' have positive random-walk SDs specified, but are included in neither 'pars' nor 'ivps'. These random walk SDs will be ignored.
-2: mif warning: the variable 'x2.0' has positive random-walk SD specified, but is included in neither 'pars' nor 'ivps'. This random walk SD will be ignored.
+Warning message:
+mif warning: the variables 'alpha.3', 'x2.0' have positive random-walk SDs specified, but are included in neither 'pars' nor 'ivps'. These random walk SDs will be ignored.
>
> proc.time()
user system elapsed
- 11.740 0.092 11.978
+ 11.960 0.064 12.167
Modified: pkg/pomp/tests/pfilter.Rout.save
===================================================================
--- pkg/pomp/tests/pfilter.Rout.save 2014-03-20 21:51:17 UTC (rev 899)
+++ pkg/pomp/tests/pfilter.Rout.save 2014-03-20 21:51:41 UTC (rev 900)
@@ -18,6 +18,7 @@
> library(pomp)
Loading required package: mvtnorm
Loading required package: subplex
+Loading required package: nloptr
Loading required package: deSolve
>
> pompExample(ou2)
@@ -52,8 +53,8 @@
> print(coef(pf))
gamma mu iota beta1 beta2 beta3 beta.sd pop
2.60e+01 2.00e-02 1.00e-02 4.00e+02 4.80e+02 3.20e+02 1.00e-03 2.10e+06
- rho S.0 I.0 R.0
-6.00e-01 6.50e-02 1.00e-03 9.35e-01
+ rho overdisp S.0 I.0 R.0
+6.00e-01 1.00e+00 6.50e-02 1.00e-03 9.35e-01
> print(pf$loglik,digits=4)
[1] -947.4
>
@@ -64,8 +65,8 @@
> print(coef(pf))
gamma mu iota beta1 beta2 beta3 beta.sd pop
2.60e+01 2.00e-02 1.00e+00 4.00e+02 4.80e+02 3.20e+02 1.00e-03 2.10e+06
- rho S.0 I.0 R.0
-6.00e-01 6.50e-02 1.00e-03 9.35e-01
+ rho overdisp S.0 I.0 R.0
+6.00e-01 1.00e+00 6.50e-02 1.00e-03 9.35e-01
> print(logLik(pf),digits=4)
[1] -945.4
> plot(cond.loglik~time,data=as(pf,"data.frame"),type='l')
@@ -78,4 +79,4 @@
>
> proc.time()
user system elapsed
- 6.084 0.052 6.316
+ 6.124 0.060 6.351
More information about the pomp-commits
mailing list