[Pomp-commits] r197 - in pkg: . R data inst inst/doc inst/examples inst/include src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jan 25 23:23:05 CET 2010
Author: kingaa
Date: 2010-01-25 23:23:04 +0100 (Mon, 25 Jan 2010)
New Revision: 197
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/bsplines.R
pkg/data/euler.sir.R
pkg/inst/ChangeLog
pkg/inst/doc/advanced_topics_in_pomp.Rnw
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/examples/euler_sir.R
pkg/inst/examples/euler_sir.c
pkg/inst/include/pomp.h
pkg/src/bspline.c
pkg/src/euler_sir.c
pkg/src/pomp.h
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
- add a new interface at the C level for evaluation of periodic B-splines bases.
- rewrite the periodically-forced SIR model to use this facility.
- rewrite the R level 'periodic.bspline.basis' function using the new C codes.
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/DESCRIPTION 2010-01-25 22:23:04 UTC (rev 197)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.27-1
-Date: 2010-01-19
+Version: 0.27-2
+Date: 2010-01-25
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall
Maintainer: Aaron A. King <kingaa at umich.edu>
Description: Inference methods for partially-observed Markov processes
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/NAMESPACE 2010-01-25 22:23:04 UTC (rev 197)
@@ -1,6 +1,7 @@
useDynLib(
pomp,
bspline_basis,
+ periodic_bspline_basis,
bspline_basis_function,
systematic_resampling,
euler_model_simulator,
Modified: pkg/R/bsplines.R
===================================================================
--- pkg/R/bsplines.R 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/R/bsplines.R 2010-01-25 22:23:04 UTC (rev 197)
@@ -13,13 +13,5 @@
}
periodic.bspline.basis <- function (x, nbasis, degree = 3, period = 1) {
- if (nbasis<degree)
- stop("periodic.bspline.basis error: must have ",sQuote("nbasis")," >= ",sQuote("degree"),call.=FALSE)
- dx <- period/nbasis
- knots <- seq(-degree*dx,period+degree*dx,by=dx)
- y <- .Call(bspline_basis,x%%period,degree,knots)
- if (degree>0)
- y[,1:degree] <- y[,1:degree]+y[,-(1:nbasis)]
- shift <- floor((degree-1)/2)
- y[,c(seq(from=shift+1,to=nbasis,by=1),seq_len(shift))]
+ .Call(periodic_bspline_basis,x,nbasis,degree,period)
}
Modified: pkg/data/euler.sir.R
===================================================================
--- pkg/data/euler.sir.R 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/data/euler.sir.R 2010-01-25 22:23:04 UTC (rev 197)
@@ -6,16 +6,13 @@
times=seq(1/52,4,by=1/52),
data=rbind(measles=numeric(52*4)),
t0=0,
- tcovar=seq(0,50,by=1/52),
- covar=matrix(
- periodic.bspline.basis(seq(0,50,by=1/52),nbasis=3,period=1,degree=3),
- ncol=3,
- dimnames=list(NULL,paste("seas",1:3,sep=''))
- ),
delta.t=1/52/20,
statenames=c("S","I","R","cases","W"),
- paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho"),
- covarnames=c("seas1"),
+ paramnames=c(
+ "gamma","mu","iota",
+ "beta1","beta.sd","pop","rho",
+ "nbasis","degree","period"
+ ),
zeronames=c("cases"),
comp.names=c("S","I","R"),
step.fun="sir_euler_simulator",
@@ -35,16 +32,15 @@
}
)
- coef(po) <- log(
- c(
- gamma=26,mu=0.02,iota=0.01,
- beta1=1200,beta2=1800,beta3=600,
- beta.sd=1e-3,
- pop=2.1e6,
- rho=0.6,
- S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200
- )
- )
+ coef(po) <- c(
+ gamma=log(26),mu=log(0.02),iota=log(0.01),
+ nbasis=3,degree=3,period=1, # NB: all parameters are log-transformed but these
+ beta1=log(1200),beta2=log(1800),beta3=log(600),
+ beta.sd=log(1e-3),
+ pop=log(2.1e6),
+ rho=log(0.6),
+ S.0=log(26/1200),I.0=log(0.001),R.0=log(1-0.001-26/1200)
+ )
simulate(po,nsim=1,seed=329348545L)
}
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/inst/ChangeLog 2010-01-25 22:23:04 UTC (rev 197)
@@ -1,5 +1,8 @@
2010-01-19 kingaa
+ * [r196] inst/ChangeLog, inst/doc/advanced_topics_in_pomp.pdf,
+ inst/doc/intro_to_pomp.pdf: - update the vignettes for version
+ 0.27-1
* [r195] DESCRIPTION, inst/ChangeLog: - version 0.27-1
* [r194] DESCRIPTION, R/mif.R, R/pfilter-mif.R, R/pfilter.R,
inst/ChangeLog, inst/doc/advanced_topics_in_pomp.pdf,
Modified: pkg/inst/doc/advanced_topics_in_pomp.Rnw
===================================================================
--- pkg/inst/doc/advanced_topics_in_pomp.Rnw 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/inst/doc/advanced_topics_in_pomp.Rnw 2010-01-25 22:23:04 UTC (rev 197)
@@ -238,51 +238,44 @@
As an additonal wrinkle, we'll assume that the rate of the infection process $\beta\,I/N$ is perturbed by white noise.
<<>>=
-euler.sir <- pomp(
- times=seq(1/52,4,by=1/52),
- data=rbind(measles=numeric(52*4)),
- t0=0,
- tcovar=seq(0,25,by=1/52),
- covar=matrix(
- periodic.bspline.basis(seq(0,25,by=1/52),nbasis=3,period=1,degree=3),
- ncol=3,
- dimnames=list(NULL,paste("seas",1:3,sep=''))
- ),
- delta.t=1/52/20,
- statenames=c("S","I","R","cases","W"),
- paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho"),
- covarnames=c("seas1"),
- zeronames=c("cases"),
- comp.names=c("S","I","R"),
- step.fun="sir_euler_simulator",
- rprocess=euler.simulate,
- skeleton.vectorfield="sir_ODE",
- rmeasure="binom_rmeasure",
- dmeasure="binom_dmeasure",
- PACKAGE="pomp",
- initializer=function(params, t0, comp.names, ...){
- p <- exp(params)
- snames <- c("S","I","R","cases","W")
- fracs <- p[paste(comp.names,"0",sep=".")]
- x0 <- numeric(length(snames))
- names(x0) <- snames
- x0[comp.names] <- round(p['pop']*fracs/sum(fracs))
- x0
- }
- )
+pomp(
+ times=seq(1/52,4,by=1/52),
+ data=rbind(measles=numeric(52*4)),
+ t0=0,
+ delta.t=1/52/20,
+ statenames=c("S","I","R","cases","W"),
+ paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho","nbasis","degree","period"),
+ zeronames=c("cases"),
+ comp.names=c("S","I","R"),
+ step.fun="sir_euler_simulator",
+ rprocess=euler.simulate,
+ skeleton.vectorfield="sir_ODE",
+ rmeasure="binom_rmeasure",
+ dmeasure="binom_dmeasure",
+ PACKAGE="pomp",
+ initializer=function(params, t0, comp.names, ...){
+ p <- exp(params)
+ snames <- c("S","I","R","cases","W")
+ fracs <- p[paste(comp.names,"0",sep=".")]
+ x0 <- numeric(length(snames))
+ names(x0) <- snames
+ x0[comp.names] <- round(p['pop']*fracs/sum(fracs))
+ x0
+ }
+ ) -> euler.sir
@
<<>>=
-coef(euler.sir) <- log(
- c(
- gamma=26,mu=0.02,iota=0.01,
- beta1=1200,beta2=1800,beta3=600,
- beta.sd=1e-3,
- pop=2.1e6,
- rho=0.6,
- S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200
- )
- )
+coef(euler.sir) <- c(
+ gamma=log(26),mu=log(0.02),iota=log(0.01),
+ beta1=log(1200),beta2=log(1800),beta3=log(600),
+ beta.sd=log(1e-3),
+ pop=log(2.1e6),
+ rho=log(0.6),
+ S.0=log(26/1200),I.0=log(0.001),R.0=log(1-0.001-26/1200),
+ nbasis=3,degree=3,period=1
+ )
+@
<<>>=
euler.sir <- simulate(euler.sir,nsim=1,seed=329348545L)
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/examples/euler_sir.R
===================================================================
--- pkg/inst/examples/euler_sir.R 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/inst/examples/euler_sir.R 2010-01-25 22:23:04 UTC (rev 197)
@@ -118,12 +118,9 @@
times=seq(1/52,4,by=1/52),
data=rbind(measles=numeric(52*4)),
t0=0,
- tcovar=tbasis,
- covar=basis,
delta.t=1/52/20,
statenames=c("S","I","R","cases","W"),
- paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho"),
- covarnames=c("seas1"),
+ paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho","nbasis","degree","period"),
zeronames=c("cases"),
step.fun="sir_euler_simulator",
rprocess=euler.simulate,
@@ -152,13 +149,13 @@
## compute a trajectory of the deterministic skeleton
tic <- Sys.time()
- X <- trajectory(po,params=log(params),hmax=1/52)
+ X <- trajectory(po,c(log(params),nbasis=3,degree=3,period=1),hmax=1/52)
toc <- Sys.time()
print(toc-tic)
## simulate from the model
tic <- Sys.time()
- x <- simulate(po,params=log(params),nsim=3)
+ x <- simulate(po,params=c(log(params),nbasis=3,degree=3,period=1),nsim=3)
toc <- Sys.time()
print(toc-tic)
Modified: pkg/inst/examples/euler_sir.c
===================================================================
--- pkg/inst/examples/euler_sir.c 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/inst/examples/euler_sir.c 2010-01-25 22:23:04 UTC (rev 197)
@@ -30,6 +30,9 @@
#define LOGBETA_SD (p[parindex[4]]) // environmental stochasticity SD in transmission rate
#define LOGPOPSIZE (p[parindex[5]]) // population size
#define LOGRHO (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
@@ -37,8 +40,6 @@
#define CASE (x[stateindex[3]]) // number of cases (accumulated per reporting period)
#define W (x[stateindex[4]]) // integrated white noise
-#define SEASBASIS (covar[covindex[0]]) // first column of seasonality basis in lookup table
-
#define MEASLES (y[0])
void binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
@@ -70,6 +71,9 @@
double gamma, mu, iota, beta_sd, beta_var, popsize;
double beta;
double dW;
+ int nseas = (int) NBASIS; // number of seasonal basis functions
+ int deg = (int) DEGREE; // degree of seasonal basis functions
+ double seasonality[nseas];
// untransform the parameters
gamma = exp(LOGGAMMA);
@@ -79,7 +83,9 @@
popsize = exp(LOGPOPSIZE);
beta_var = beta_sd*beta_sd;
- beta = exp(dot_product(covdim,&SEASBASIS,&LOGBETA));
+ if (nseas <= 0) return;
+ periodic_bspline_basis_eval(t,PERIOD,deg,nseas,&seasonality[0]);
+ beta = exp(dot_product(nseas,&seasonality[0],&LOGBETA));
// test to make sure the parameters and state variable values are sane
if (!(R_FINITE(beta)) ||
@@ -141,6 +147,9 @@
double term[nrate]; // terms in the equations
double gamma, mu, iota, popsize;
double beta;
+ int nseas = (int) NBASIS; // number of seasonal basis functions
+ int deg = (int) DEGREE; // degree of seasonal basis functions
+ double seasonality[nseas];
// untransform the parameters
gamma = exp(LOGGAMMA);
@@ -148,7 +157,9 @@
iota = exp(LOGIOTA);
popsize = exp(LOGPOPSIZE);
- beta = exp(dot_product(covdim,&SEASBASIS,&LOGBETA));
+ if (nseas <= 0) return;
+ periodic_bspline_basis_eval(t,PERIOD,deg,nseas,&seasonality[0]);
+ beta = exp(dot_product(nseas,&seasonality[0],&LOGBETA));
// compute the transition rates
rate[0] = mu*popsize; // birth into susceptible class
@@ -185,8 +196,6 @@
#undef CASE
#undef W
-#undef SEASBASIS
-
#undef LOGGAMMA
#undef LOGMU
#undef LOGIOTA
@@ -194,3 +203,6 @@
#undef LOGBETA_SD
#undef LOGPOPSIZE
#undef LOGRHO
+#undef NBASIS
+#undef DEGREE
+#undef PERIOD
Modified: pkg/inst/include/pomp.h
===================================================================
--- pkg/inst/include/pomp.h 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/inst/include/pomp.h 2010-01-25 22:23:04 UTC (rev 197)
@@ -15,6 +15,9 @@
// facility for dotting a vector of parameters ('coef') against a vector of basis-function values ('basis')
double dot_product (int dim, const double *basis, const double *coef);
+// facility for computing evaluating a basis of periodic bsplines
+void periodic_bspline_basis_eval (double x, double period, int degree, int nbasis, double *y);
+
// Prototype for one-step simulator, as used by "euler.simulate" and "onestep.simulate":
typedef void pomp_onestep_sim(double *x, const double *p,
const int *stateindex, const int *parindex, const int *covindex,
Modified: pkg/src/bspline.c
===================================================================
--- pkg/src/bspline.c 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/src/bspline.c 2010-01-25 22:23:04 UTC (rev 197)
@@ -2,28 +2,28 @@
#include "pomp_internal.h"
+void periodic_bspline_basis_eval (double x, double period, int degree, int nbasis, double *y);
static void bspline_internal (double *y, const double *x, int nx, int i, int p, const double *knots, int nknots);
// B-spline basis
SEXP bspline_basis (SEXP x, SEXP degree, SEXP knots) {
int nprotect = 0;
- SEXP y, xr, kr, di;
+ SEXP y, xr, kr;
int nx = length(x);
int nknots = length(knots);
- int deg;
- int nbasis;
+ int deg = INTEGER_VALUE(degree);
+ int nb;
double *ydata;
int i;
- PROTECT(di = AS_INTEGER(degree)); nprotect++;
- deg = INTEGER_VALUE(di);
- if (deg < 0) error("must have degree > 0 in 'bspline.basis'");
- nbasis = nknots-deg-1;
+ nb = nknots-deg-1;
+ if (deg < 0) error("bspline.basis error: must have degree > 0");
+ if (nb<=deg) error("bspline.basis error: must have nbasis > degree");
PROTECT(xr = AS_NUMERIC(x)); nprotect++;
PROTECT(kr = AS_NUMERIC(knots)); nprotect++;
- PROTECT(y = allocMatrix(REALSXP,nx,nbasis)); nprotect++;
+ PROTECT(y = allocMatrix(REALSXP,nx,nb)); nprotect++;
ydata = REAL(y);
- for (i = 0; i < nbasis; i++) {
+ for (i = 0; i < nb; i++) {
bspline_internal(ydata,REAL(xr),nx,i,deg,REAL(kr),nknots);
ydata += nx;
}
@@ -49,6 +49,61 @@
}
+SEXP periodic_bspline_basis (SEXP x, SEXP nbasis, SEXP degree, SEXP period) {
+ int nprotect = 0;
+ SEXP y, xr;
+ int nx = length(x);
+ int nb = INTEGER_VALUE(nbasis);
+ int deg = INTEGER_VALUE(degree);
+ double pd = NUMERIC_VALUE(period);
+ int j, k;
+ double *xrd, *ydata, *val;
+ if (deg < 0) error("periodic_bspline_basis error: must have degree >= 0");
+ if (nb <= 0) error("periodic_bspline_basis error: must have nbasis > 0");
+ if (deg > nb) error("periodic_bspline_basis error: must have degree <= nbasis");
+ if (pd <= 0.0) error("periodic_bspline_basis error: must have period > 0");
+ PROTECT(xr = AS_NUMERIC(x)); nprotect++;
+ xrd = REAL(xr);
+ PROTECT(y = allocMatrix(REALSXP,nx,nb)); nprotect++;
+ ydata = REAL(y);
+ val = (double *) Calloc(nb,double);
+ for (j = 0; j < nx; j++) {
+ periodic_bspline_basis_eval(xrd[j],pd,deg,nb,val);
+ for (k = 0; k < nb; k++) ydata[j+nx*k] = val[k];
+ }
+ Free(val);
+ UNPROTECT(nprotect);
+ return y;
+}
+
+void periodic_bspline_basis_eval (double x, double period, int degree, int nbasis, double *y)
+{
+ int nknots = nbasis+2*degree+1;
+ int shift = (degree-1)/2;
+ double knots[nknots];
+ double yy[nknots];
+ double dx;
+ int j, k;
+ if (period <= 0.0) error("periodic_bspline_basis_eval error: must have period > 0");
+ if (nbasis <= 0) error("periodic_bspline_basis_eval error: must have nbasis > 0");
+ if (degree < 0) error("periodic_bspline_basis_eval error: must have degree >= 0");
+ if (nbasis < degree) error("periodic_bspline_basis_eval error: must have nbasis >= degree");
+ dx = period/((double) nbasis);
+ for (k = -degree; k <= nbasis+degree; k++) {
+ knots[degree+k] = k*dx;
+ }
+ x = fmod(x,period);
+ if (x < 0.0) x += period;
+ for (k = 0; k < nknots; k++) {
+ bspline_internal(&yy[k],&x,1,k,degree,&knots[0],nknots);
+ }
+ for (k = 0; k < degree; k++) yy[k] += yy[nbasis+k];
+ for (k = 0; k < nbasis; k++) {
+ j = (shift+k)%nbasis;
+ y[k] = yy[j];
+ }
+}
+
static void bspline_internal (double *y, const double *x, int nx, int i, int p, const double *knots, int nknots)
{
int j;
Modified: pkg/src/euler_sir.c
===================================================================
--- pkg/src/euler_sir.c 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/src/euler_sir.c 2010-01-25 22:23:04 UTC (rev 197)
@@ -30,6 +30,9 @@
#define LOGBETA_SD (p[parindex[4]]) // environmental stochasticity SD in transmission rate
#define LOGPOPSIZE (p[parindex[5]]) // population size
#define LOGRHO (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
@@ -37,8 +40,6 @@
#define CASE (x[stateindex[3]]) // number of cases (accumulated per reporting period)
#define W (x[stateindex[4]]) // integrated white noise
-#define SEASBASIS (covar[covindex[0]]) // first column of seasonality basis in lookup table
-
#define MEASLES (y[0])
void binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
@@ -70,6 +71,9 @@
double gamma, mu, iota, beta_sd, beta_var, popsize;
double beta;
double dW;
+ int nseas = (int) NBASIS; // number of seasonal basis functions
+ int deg = (int) DEGREE; // degree of seasonal basis functions
+ double seasonality[nseas];
// untransform the parameters
gamma = exp(LOGGAMMA);
@@ -79,7 +83,9 @@
popsize = exp(LOGPOPSIZE);
beta_var = beta_sd*beta_sd;
- beta = exp(dot_product(covdim,&SEASBASIS,&LOGBETA));
+ if (nseas <= 0) return;
+ periodic_bspline_basis_eval(t,PERIOD,deg,nseas,&seasonality[0]);
+ beta = exp(dot_product(nseas,&seasonality[0],&LOGBETA));
// test to make sure the parameters and state variable values are sane
if (!(R_FINITE(beta)) ||
@@ -141,6 +147,9 @@
double term[nrate]; // terms in the equations
double gamma, mu, iota, popsize;
double beta;
+ int nseas = (int) NBASIS; // number of seasonal basis functions
+ int deg = (int) DEGREE; // degree of seasonal basis functions
+ double seasonality[nseas];
// untransform the parameters
gamma = exp(LOGGAMMA);
@@ -148,7 +157,9 @@
iota = exp(LOGIOTA);
popsize = exp(LOGPOPSIZE);
- beta = exp(dot_product(covdim,&SEASBASIS,&LOGBETA));
+ if (nseas <= 0) return;
+ periodic_bspline_basis_eval(t,PERIOD,deg,nseas,&seasonality[0]);
+ beta = exp(dot_product(nseas,&seasonality[0],&LOGBETA));
// compute the transition rates
rate[0] = mu*popsize; // birth into susceptible class
@@ -185,8 +196,6 @@
#undef CASE
#undef W
-#undef SEASBASIS
-
#undef LOGGAMMA
#undef LOGMU
#undef LOGIOTA
@@ -194,3 +203,6 @@
#undef LOGBETA_SD
#undef LOGPOPSIZE
#undef LOGRHO
+#undef NBASIS
+#undef DEGREE
+#undef PERIOD
Modified: pkg/src/pomp.h
===================================================================
--- pkg/src/pomp.h 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/src/pomp.h 2010-01-25 22:23:04 UTC (rev 197)
@@ -15,6 +15,9 @@
// facility for dotting a vector of parameters ('coef') against a vector of basis-function values ('basis')
double dot_product (int dim, const double *basis, const double *coef);
+// facility for computing evaluating a basis of periodic bsplines
+void periodic_bspline_basis_eval (double x, double period, int degree, int nbasis, double *y);
+
// Prototype for one-step simulator, as used by "euler.simulate" and "onestep.simulate":
typedef void pomp_onestep_sim(double *x, const double *p,
const int *stateindex, const int *parindex, const int *covindex,
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/tests/sir.R 2010-01-25 22:23:04 UTC (rev 197)
@@ -211,10 +211,10 @@
x=X1$states,
times=t1,
params=matrix(
- log(params),
- nrow=length(params),
+ c(log(params),nbasis=3,degree=3,period=1),
+ nrow=length(params)+3,
ncol=10,
- dimnames=list(names(params),NULL)
+ dimnames=list(c(names(params),"nbasis","degree","period"),NULL)
),
log=TRUE
)
@@ -224,7 +224,7 @@
po,
x=X2$states[,1,55:70,drop=FALSE],
t=t2[55:70],
- params=as.matrix(log(params))
+ params=as.matrix(c(log(params),nbasis=3,degree=3,period=1))
)
print(h2[c("S","I","R"),,],digits=4)
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2010-01-19 12:41:13 UTC (rev 196)
+++ pkg/tests/sir.Rout.save 2010-01-25 22:23:04 UTC (rev 197)
@@ -423,7 +423,7 @@
}
y
}
-<environment: 0x185b800>
+<environment: 0x28d8c68>
measurement model density, dmeasure =
function (y, x, t, params, log, covars, ...)
{
@@ -436,7 +436,7 @@
f
else exp(f)
}
-<environment: 0x185b800>
+<environment: 0x28d8c68>
initializer =
function (params, t0, ...)
{
@@ -507,7 +507,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 2.009385 secs
+Time difference of 2.041221 secs
>
> pdf(file='sir.pdf')
>
@@ -524,7 +524,7 @@
> X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 4.488558 secs
+Time difference of 4.547439 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
@@ -584,7 +584,7 @@
> x <- simulate(po,nsim=100)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 1.264836 secs
+Time difference of 2.148207 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
@@ -592,7 +592,7 @@
> X4 <- trajectory(po,times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 3.623464 secs
+Time difference of 0.70203 secs
> plot(t3,X4['I',1,],type='l')
>
> g2 <- dmeasure(
@@ -601,10 +601,10 @@
+ x=X1$states,
+ times=t1,
+ params=matrix(
-+ log(params),
-+ nrow=length(params),
++ c(log(params),nbasis=3,degree=3,period=1),
++ nrow=length(params)+3,
+ ncol=10,
-+ dimnames=list(names(params),NULL)
++ dimnames=list(c(names(params),"nbasis","degree","period"),NULL)
+ ),
+ log=TRUE
+ )
@@ -615,7 +615,7 @@
+ po,
+ x=X2$states[,1,55:70,drop=FALSE],
+ t=t2[55:70],
-+ params=as.matrix(log(params))
++ params=as.matrix(c(log(params),nbasis=3,degree=3,period=1))
+ )
> print(h2[c("S","I","R"),,],digits=4)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
@@ -648,12 +648,12 @@
cases 0 0 1.045000e+03
W 0 0 -2.334331e-01
> states(simulate(po))[,1:3]
- [,1] [,2] [,3]
-S 45500 45500 4.527200e+04
-I 2100 2100 2.030000e+03
-R 2052400 2052400 2.052644e+06
-cases 0 0 1.034000e+03
-W 0 0 1.415062e-01
+ [,1] [,2] [,3]
+S 45500 45500 4.528700e+04
+I 2100 2100 2.062000e+03
+R 2052400 2052400 2.052647e+06
+cases 0 0 1.007000e+03
+W 0 0 -1.602887e-01
>
> dev.off()
null device
More information about the pomp-commits
mailing list