[Pomp-commits] r1038 - in pkg/pompExamples: . inst/examples src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jan 1 22:20:25 CET 2015
Author: kingaa
Date: 2015-01-01 22:20:25 +0100 (Thu, 01 Jan 2015)
New Revision: 1038
Added:
pkg/pompExamples/inst/examples/ebola.R
pkg/pompExamples/src/cholmodel.c
pkg/pompExamples/src/ebola.c
Modified:
pkg/pompExamples/DESCRIPTION
pkg/pompExamples/tests/examples.R
pkg/pompExamples/tests/pertussis.R
pkg/pompExamples/tests/pertussis.Rout.save
Log:
- Ebola example
- put Matteo Fasiolo's variant cholera model source code into the package
Modified: pkg/pompExamples/DESCRIPTION
===================================================================
--- pkg/pompExamples/DESCRIPTION 2014-12-31 12:36:42 UTC (rev 1037)
+++ pkg/pompExamples/DESCRIPTION 2015-01-01 21:20:25 UTC (rev 1038)
@@ -1,8 +1,8 @@
Package: pompExamples
Type: Package
Title: Additional pomp examples
-Version: 0.24-1
-Date: 2014-12-31
+Version: 0.25-1
+Date: 2015-01-01
Maintainer: Aaron A. King <kingaa at umich.edu>
Authors at R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),
email="kingaa at umich.edu"),
Added: pkg/pompExamples/inst/examples/ebola.R
===================================================================
--- pkg/pompExamples/inst/examples/ebola.R (rev 0)
+++ pkg/pompExamples/inst/examples/ebola.R 2015-01-01 21:20:25 UTC (rev 1038)
@@ -0,0 +1,122 @@
+require(pomp)
+require(plyr)
+require(reshape2)
+
+WHO.situation.report.Oct.1 <- '
+Week,Guinea,Liberia,SierraLeone
+1,2.244,,
+2,2.244,,
+3,0.073,,
+4,5.717,,
+5,3.954,,
+6,5.444,,
+7,3.274,,
+8,5.762,,
+9,7.615,,
+10,7.615,,
+11,27.392,,
+12,17.387,,
+13,27.115,,
+14,29.29,,
+15,27.84,,
+16,16.345,,
+17,10.917,,
+18,11.959,,
+19,11.959,,
+20,8.657,,
+21,26.537,,
+22,47.764,3.517,
+23,26.582,1.043,5.494
+24,32.967,18,57.048
+25,18.707,16.34,76.022
+26,24.322,13.742,36.768
+27,4.719,10.155,81.929
+28,7.081,24.856,102.632
+29,8.527,53.294,69.823
+30,92.227,70.146,81.783
+31,26.423,139.269,99.775
+32,16.549,65.66,88.17
+33,36.819,240.645,90.489
+34,92.08,274.826,161.54
+35,101.03,215.56,168.966
+36,102.113,388.553,186.144
+37,83.016,410.299,220.442
+38,106.674,300.989,258.693
+39,55.522,240.237,299.546
+'
+
+## Population sizes in Guinea, Liberia, and Sierra Leone (census 2014)
+populations <- c(Guinea=10628972,Liberia=4092310,SierraLeone=6190280)
+populations["WestAfrica"] <- sum(populations)
+
+dat <- read.csv(text=WHO.situation.report.Oct.1,stringsAsFactors=FALSE)
+rename(dat,c(Week="week")) -> dat
+dat <- melt(dat,id="week",variable.name="country",value.name="cases")
+mutate(dat,deaths=NA) -> dat
+
+ebolaModel <- function (country=c("Guinea", "SierraLeone", "Liberia", "WestAfrica"),
+ timestep = 0.01, nstageE = 3L,
+ type = c("raw","cum"), na.rm = FALSE, least.sq = FALSE) {
+
+ type <- match.arg(type)
+ ctry <- match.arg(country)
+ pop <- unname(populations[ctry])
+
+ ## Incubation period is supposed to be Gamma distributed with shape parameter 3 and mean 11.4 days
+ ## The discrete-time formula is used to calculate the corresponding alpha (cf He et al., Interface 2010)
+ ## Case-fatality ratio is fixed at 0.7 (cf WHO Ebola response team, NEJM 2014)
+ incubation_period <- 11.4/7
+ infectious_period <- 7/7
+ index_case <- 10/pop
+ dt <- timestep
+
+ theta <- c(N=pop,R0=1.4,
+ alpha=-1/(nstageE*dt)*log(1-nstageE*dt/incubation_period),
+ gamma=-log(1-dt/infectious_period)/dt,
+ rho=0.2,cfr=0.7,
+ k=0,
+ S_0=1-index_case,E_0=index_case/2-5e-9,
+ I_0=index_case/2-5e-9,R_0=1e-8)
+
+ dat <- subset(dat,country==ctry,select=-country)
+ if (na.rm) {
+ dat <- mutate(subset(dat,!is.na(cases)),week=week-min(week)+1)
+ }
+ if (type=="cum") {
+ dat <- mutate(dat,cases=cumsum(cases),deaths=cumsum(deaths))
+ }
+
+ print(dat)
+
+ ## Create the pomp object
+ pomp(
+ data=dat,
+ times="week",
+ t0=0,
+ params=theta,
+ obsnames=c("cases","deaths"),
+ statenames=c("S","E1","I","R","N_EI","N_IR"),
+ zeronames=if (type=="raw") c("N_EI","N_IR") else character(0),
+ paramnames=c("N","R0","alpha","gamma","rho","k","cfr",
+ "S_0","E_0","I_0","R_0"),
+ nstageE=nstageE,
+ PACKAGE="pompExamples",
+ dmeasure=if (least.sq) "_ebola_dObsLS" else "_ebola_dObs",
+ rmeasure=if (least.sq) "_ebola_rObsLS" else "_ebola_rObs",
+ rprocess=discrete.time.sim(step.fun="_ebola_rSim",delta.t=timestep),
+ skeleton="_ebola_skel",
+ skeleton.type="vectorfield",
+ parameter.transform="_ebola_par_trans",
+ parameter.inv.transform="_ebola_par_untrans",
+ initializer=function (params, t0, nstageE, ...) {
+ all.state.names <- c("S",paste0("E",1:nstageE),"I","R","N_EI","N_IR")
+ comp.names <- c("S",paste0("E",1:nstageE),"I","R")
+ x0 <- setNames(numeric(length(all.state.names)),all.state.names)
+ frac <- c(params["S_0"],rep(params["E_0"]/nstageE,nstageE),params["I_0"],params["R_0"])
+ x0[comp.names] <- round(params["N"]*frac/sum(frac))
+ x0
+ }
+ ) -> po
+}
+
+c("ebolaModel")
Added: pkg/pompExamples/src/cholmodel.c
===================================================================
--- pkg/pompExamples/src/cholmodel.c (rev 0)
+++ pkg/pompExamples/src/cholmodel.c 2015-01-01 21:20:25 UTC (rev 1038)
@@ -0,0 +1,238 @@
+// -*- C++ -*-
+
+#include <R.h>
+#include <Rmath.h>
+#include "pomp.h"
+
+#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 LOGOMEGA (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 S0 (p[parindex[14]])
+#define I0 (p[parindex[15]])
+#define RS0 (p[parindex[16]])
+#define RL0 (p[parindex[17]])
+
+#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_untrans (double *pt, double *p, int *parindex)
+{
+ int k, nrstage = (int) NRSTAGE;
+ pt[parindex[0]] = log(TAU);
+ pt[parindex[1]] = log(GAMMA);
+ pt[parindex[2]] = log(EPS);
+ pt[parindex[3]] = log(DELTA);
+ pt[parindex[4]] = log(DELTA_I);
+ pt[parindex[6]] = log(SD_BETA);
+ pt[parindex[9]] = log(ALPHA);
+ pt[parindex[10]] = log(RHO);
+ pt[parindex[11]] = logit(CLIN);
+
+ to_log_barycentric(&pt[parindex[14]],&S0,3+nrstage);
+}
+
+void _cholmodel_trans (double *pt, double *p, int *parindex)
+{
+ int k, nrstage = (int) NRSTAGE;
+ pt[parindex[0]] = exp(TAU);
+ pt[parindex[1]] = exp(GAMMA);
+ pt[parindex[2]] = exp(EPS);
+ pt[parindex[3]] = exp(DELTA);
+ pt[parindex[4]] = exp(DELTA_I);
+ pt[parindex[6]] = exp(SD_BETA);
+ pt[parindex[9]] = exp(ALPHA);
+ pt[parindex[10]] = exp(RHO);
+ pt[parindex[11]] = expit(CLIN);
+
+ from_log_barycentric(&pt[parindex[14]],&S0,3+nrstage);
+}
+
+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 tol = 1.0e-18;
+ if ((COUNT > 0) || (!(R_FINITE(DEATHS)))) {
+ DATADEATHS = R_NaReal;
+ } else {
+ DATADEATHS = rnbinom_mu(1 / (TAU*TAU), DEATHS + 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 tol = 1.0e-18;
+ if ((COUNT>0.0) || (!(R_FINITE(DEATHS)))) {
+ *lik = tol;
+ } else {
+ *lik = dnbinom_mu(DATADEATHS+tol, 1 / (TAU*TAU), DEATHS + tol, 0) + tol;
+ }
+ if (give_log) *lik = log(*lik);
+}
+
+#undef DATADEATHS
+
+
+inline double out_flow(double init, double rate_by_dt)
+{
+
+ return init * (1.0 - exp( - rate_by_dt ));
+
+}
+
+// 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_two (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 *pt_r, *pt_w;
+ int j;
+
+ double k_eps = NRSTAGE * EPS;
+
+ double gamma_delta_deltaI = GAMMA + DELTA + DELTA_I;
+
+ double keps_div_keps_delta = k_eps / (k_eps + DELTA);
+
+ /*
+ * Preliminaries
+ */
+
+ double beta = exp(dot_product(nbasis,&SEASBASIS,&LOGBETA)+BETATREND*TREND);
+
+ double omega = exp(dot_product(nbasis,&SEASBASIS,&LOGOMEGA));
+
+ double dw = rgammawn(SD_BETA, dt); // gamma noise, mean=dt, variance=(beta_sd^2 dt)
+
+ double lambda = omega + (beta * dw/dt) * INFECT/POP; // Time-dependent force of infection
+
+ /*
+ * Out-flows from compartments using aboundances at previous time step
+ */
+
+ double S_o = out_flow(SUSCEP, (lambda + DELTA) * dt);
+ double I_o = out_flow(INFECT, gamma_delta_deltaI * dt);
+ double RSH_o = out_flow(RSHORT, (RHO + DELTA) * dt);
+
+ double R_o[nrstage];
+ for (pt_r = &RLONG, pt_w = R_o, j = 0;
+ j < nrstage;
+ j++, pt_r++, pt_w++) {
+
+ *pt_w = out_flow(*pt_r, (k_eps + DELTA) * dt);
+
+ }
+
+ double new_infect = S_o * (lambda / (lambda + DELTA));
+
+ /*
+ * Cholera deaths = (I_o * DELTA_I) / gamma_delta_deltaI are not reborn in SUSCEPT.
+ * To offset this downward bias on total population, we add the average number of (observed) deaths:
+ *
+ * average_chol_death = tot_cholc_deaths / (n_months * n_step) <-----> 29.522 = 354275.0 / (600 * 200)
+ */
+ double chol_deaths_dt = 29.523;
+ double births = DPOPDT * dt + chol_deaths_dt + S_o * (DELTA / (lambda + DELTA)) + I_o * DELTA / gamma_delta_deltaI + RSH_o * (DELTA / (DELTA + RHO));
+ for (pt_r = R_o, j = 0;
+ j < nrstage;
+ j++, pt_r++) {
+
+ births += *pt_r * (DELTA / (k_eps + DELTA));
+
+ }
+
+ /*
+ * Updating all compartments
+ */
+
+ SUSCEP += - S_o + births + R_o[nrstage-1] * keps_div_keps_delta + RSH_o * (RHO / (RHO + DELTA));
+
+ INFECT += - I_o + CLIN * new_infect;
+
+ RSHORT += - RSH_o + (1 - CLIN) * new_infect;
+
+ for (pt_w = &RLONG, pt_r = R_o, j = 0;
+ j < nrstage;
+ j++, pt_r++, pt_w++) {
+
+ if(j == 0)
+ {
+
+ *pt_w += - *pt_r + I_o * (GAMMA / gamma_delta_deltaI);
+
+ } else {
+
+ *pt_w += - *pt_r + *(pt_r-1) * keps_div_keps_delta;
+
+ }
+
+ }
+
+ DEATHS += I_o * (DELTA_I / gamma_delta_deltaI); // cumulative deaths due to disease
+
+ NOISE += (dw-dt) / SD_BETA;
+
+}
+
+#undef GAMMA
+#undef EPS
+#undef DELTA
+#undef DELTA_I
+#undef LOGOMEGA
+#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
Added: pkg/pompExamples/src/ebola.c
===================================================================
--- pkg/pompExamples/src/ebola.c (rev 0)
+++ pkg/pompExamples/src/ebola.c 2015-01-01 21:20:25 UTC (rev 1038)
@@ -0,0 +1,187 @@
+// SEIR Ebola model
+
+#include <R.h>
+#include <Rmath.h>
+#include <R_ext/Rdynload.h>
+#include <pomp.h>
+
+// State variables
+#define S (x[stateindex[0]]) // Susceptible
+#define E(J) (x[stateindex[1] + (J)]) // Exposed
+#define I (x[stateindex[2]]) // Infected
+#define R (x[stateindex[3]]) // Removed
+#define N_EI (x[stateindex[4]]) // Number of transitions from E to I
+#define N_IR (x[stateindex[5]]) // Number of transitions from I to R
+
+// Variations
+#define DS (f[stateindex[0]]) // Susceptible
+#define DE(J) (f[stateindex[1] + (J)]) // Exposed
+#define DI (f[stateindex[2]]) // Infected
+#define DR (f[stateindex[3]]) // Removed
+#define DN_EI (f[stateindex[4]]) // Number of transitions from E to I
+#define DN_IR (f[stateindex[5]]) // Number of transitions from I to R
+
+// Parameters on the natural scale (all rates are per day)
+#define N (p[parindex[0]]) // Population size
+#define R0 (p[parindex[1]]) // Basic reproduction number
+#define alpha (p[parindex[2]]) // Inverse of latency period
+#define gamma (p[parindex[3]]) // Inverse of duration of infection
+#define rho (p[parindex[4]]) // Reporting probability
+#define k (p[parindex[5]]) // Reporting overdispersion
+#define cfr (p[parindex[6]]) // Case-fatality ratio
+#define IC(J) (p[parindex[7] + (J)]) // Initial conditions
+
+// Parameters on the transformed scale (all rates are per day)
+#define TN (pt[parindex[0]]) // Population size
+#define TR0 (pt[parindex[1]]) // Basic reproduction number
+#define Talpha (pt[parindex[2]]) // Inverse of latency period
+#define Tgamma (pt[parindex[3]]) // Inverse of duration of infection
+#define Trho (pt[parindex[4]]) // Reporting probability
+#define Tk (pt[parindex[5]]) // Reporting overdispersion
+#define Tcfr (pt[parindex[6]]) // Case-fatality ratio
+#define TIC(J) (pt[parindex[7] + (J)]) // Initial conditions
+
+// Observations
+#define cases (y[obsindex[0]]) // Number of reported cases
+#define deaths (y[obsindex[1]]) // Number of reported deaths
+
+// Transforms the parameters to the transformed scale
+void _ebola_par_untrans (double *pt, double *p, int *parindex){
+ TN = log(N);
+ TR0 = log(R0);
+ Talpha = log(alpha);
+ Tgamma = log(gamma);
+ Trho = logit(rho);
+ Tk = log(k);
+ Tcfr = logit(cfr);
+ to_log_barycentric(&(TIC(0)),&(IC(0)),4);
+}
+
+// Transforms the parameters to the natural scale
+void _ebola_par_trans (double *pt, double *p, int *parindex){
+ TN = exp(N);
+ TR0 = exp(R0);
+ Talpha = exp(alpha);
+ Tgamma = exp(gamma);
+ Trho = expit(rho);
+ Tk = exp(k);
+ Tcfr = expit(cfr);
+ from_log_barycentric(&(TIC(0)),&(IC(0)),4);
+}
+
+// Observation model: hierarchical model for cases and deaths
+// p(R_t, D_t| C_t) = p(R_t | C_t) * p(D_t | C_t, R_t)
+// p(R_t | C_t): Negative binomial with mean rho * C_t and dispersion parameter 1 / k
+// p(D_t | C_t, R_t): Binomial B(R_t, cfr)
+void _ebola_dObs (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 f;
+ if (k > 0.0)
+ f = dnbinom_mu(nearbyint(cases),1.0/k,rho*N_EI,1);
+ // f += dnbinom_mu(nearbyint(deaths), 1.0 / k, rho * cfr * N_IR, 1);
+ else
+ f = dpois(nearbyint(cases),rho*N_EI,1);
+ *lik = (give_log) ? f : exp(f);
+}
+
+// For least-squares trajectory-matching:
+void _ebola_dObsLS (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 f;
+ f = dnorm(cases,rho*N_EI,k,1);
+ *lik = (give_log) ? f : exp(f);
+}
+
+void _ebola_rObs (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t)
+{
+ if (k > 0) {
+ cases = rnbinom_mu(1.0/k,rho*N_EI);
+ deaths = rnbinom_mu(1.0/k,rho*cfr*N_IR);
+ } else {
+ cases = rpois(rho*N_EI);
+ deaths = rpois(rho*cfr*N_IR);
+ }
+}
+
+// For least-squares trajectory-matching:
+void _ebola_rObsLS (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t)
+{
+ cases = rnorm(rho*N_EI,k);
+ deaths = NA_REAL;
+}
+
+
+// Process model
+void _ebola_rSim (double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covars,
+ double t, double dt) {
+
+ // Retrieve user data in the pomp object
+ int *(*get_pomp_userdata_int)(const char *);
+ get_pomp_userdata_int = (int *(*)(const char *)) R_GetCCallable("pomp","get_pomp_userdata_int");
+ int nstageE = *(get_pomp_userdata_int("nstageE")); // Number of stages in the E class
+
+ // Other parameters
+ double lambda, beta;
+ beta = R0 * gamma; // Transmission rate
+ lambda = beta * I / N; // Force of infection
+ int i;
+
+ // Transitions
+
+ // From class S
+ double transS = rbinom(S, 1.0 - exp(- lambda * dt)); // No of infections
+
+ // From class E
+ double transE[nstageE]; // No of transitions between classes E
+ for(i = 0; i < nstageE; i++){
+ transE[i] = rbinom(E(i), 1.0 - exp(- nstageE * alpha * dt));
+ }
+
+ // From class I
+ double transI = rbinom(I, 1.0 - exp(- gamma * dt)); // No of transitions I->R
+
+ // Balance the equations
+ S -= transS;
+ E(0) += transS - transE[0];
+ for(i=1; i < nstageE; i++) {
+ E(i) += transE[i-1] - transE[i];
+ }
+ I += transE[nstageE - 1] - transI;
+ R += transI;
+ N_EI += transE[nstageE - 1]; // No of transitions from E to I
+ N_IR += transI; // No of transitions from I to R
+}
+
+// Continuous-time deterministic skeleton
+void _ebola_skel (double *f, double *x, double *p,
+ int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+
+ // Retrieve user data in the pomp object
+ int *(*get_pomp_userdata_int)(const char *);
+ get_pomp_userdata_int = (int *(*)(const char *)) R_GetCCallable("pomp","get_pomp_userdata_int");
+ int nstageE = *(get_pomp_userdata_int("nstageE")); // Number of stages in the E class
+
+ // Other parameters
+ double lambda, beta;
+ beta = R0 * gamma; // Transmission rate
+ lambda = beta * I / N; // Force of infection
+ int i;
+
+ // Balance the equations
+ DS = - lambda * S;
+ DE(0) = lambda * S - nstageE * alpha * E(0);
+ for(i=1; i < nstageE; i++) DE(i) = nstageE * alpha * (E(i - 1) - E(i));
+ DI = nstageE * alpha * E(nstageE - 1) - gamma * I;
+ DR = gamma * I;
+ DN_EI = nstageE * alpha * E(nstageE - 1);
+ DN_IR = gamma * I;
+}
Modified: pkg/pompExamples/tests/examples.R
===================================================================
--- pkg/pompExamples/tests/examples.R 2014-12-31 12:36:42 UTC (rev 1037)
+++ pkg/pompExamples/tests/examples.R 2015-01-01 21:20:25 UTC (rev 1038)
@@ -22,3 +22,12 @@
pompExample(bbp)
pf <- pfilter(simulate(bbp),Np=100,max.fail=Inf)
tj <- trajectory(bbp)
+
+pompExample(ebola)
+ebolaModel(country="Guinea") -> po
+pf <- pfilter(simulate(po),Np=100)
+tj <- trajectory(po)
+
+ebolaModel(country="SierraLeone",na.rm=TRUE,type='cum') -> po
+pf <- pfilter(simulate(po),Np=100)
+tj <- trajectory(po)
Modified: pkg/pompExamples/tests/pertussis.R
===================================================================
--- pkg/pompExamples/tests/pertussis.R 2014-12-31 12:36:42 UTC (rev 1037)
+++ pkg/pompExamples/tests/pertussis.R 2015-01-01 21:20:25 UTC (rev 1038)
@@ -15,7 +15,7 @@
y <- trajectory(pertussis.sim(SEIRS.small),as.data.frame=TRUE)
tail(y)
-system.time(pf <- pfilter(pertussis.sim(full.small),seed=3445886L,Np=1000))
+pf <- pfilter(pertussis.sim(full.small),seed=3445886L,Np=1000)
logLik(pf)
pttest <- function (po, digits = 15) {
Modified: pkg/pompExamples/tests/pertussis.Rout.save
===================================================================
--- pkg/pompExamples/tests/pertussis.Rout.save 2014-12-31 12:36:42 UTC (rev 1037)
+++ pkg/pompExamples/tests/pertussis.Rout.save 2015-01-01 21:20:25 UTC (rev 1038)
@@ -1,5 +1,5 @@
-R Under development (unstable) (2014-12-14 r67168) -- "Unsuffered Consequences"
+R version 3.1.2 (2014-10-31) -- "Pumpkin Helmet"
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)
@@ -144,9 +144,7 @@
1040 19.98077 1
1041 20.00000 1
>
-> system.time(pf <- pfilter(pertussis.sim(full.small),seed=3445886L,Np=1000))
- user system elapsed
- 17.801 0.004 17.861
+> pf <- pfilter(pertussis.sim(full.small),seed=3445886L,Np=1000)
> logLik(pf)
[1] -3829.33
>
@@ -170,4 +168,4 @@
>
> proc.time()
user system elapsed
- 18.657 0.060 18.801
+ 23.494 0.034 23.544
More information about the pomp-commits
mailing list