[Pomp-commits] r660 - in pkg: . data inst inst/data-R inst/examples man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Apr 15 17:43:19 CEST 2012
Author: kingaa
Date: 2012-04-15 17:43:19 +0200 (Sun, 15 Apr 2012)
New Revision: 660
Modified:
pkg/DESCRIPTION
pkg/data/bbs.rda
pkg/data/blowflies.rda
pkg/data/dacca.rda
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/gompertz.rda
pkg/data/ou2.rda
pkg/data/ricker.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/NEWS
pkg/inst/data-R/sir.R
pkg/inst/examples/sir.c
pkg/man/sir.Rd
pkg/src/sir.c
pkg/tests/gillespie.Rout.save
pkg/tests/pfilter.Rout.save
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
- change the way beta is modeled in the SIR examples
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-04-14 21:19:02 UTC (rev 659)
+++ pkg/DESCRIPTION 2012-04-15 15:43:19 UTC (rev 660)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.41-5
-Date: 2012-04-14
+Version: 0.41-6
+Date: 2012-04-15
Revision: $Rev$
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing, Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
Modified: pkg/data/bbs.rda
===================================================================
(Binary files differ)
Modified: pkg/data/blowflies.rda
===================================================================
(Binary files differ)
Modified: pkg/data/dacca.rda
===================================================================
(Binary files differ)
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gompertz.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ricker.rda
===================================================================
(Binary files differ)
Modified: pkg/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/verhulst.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2012-04-14 21:19:02 UTC (rev 659)
+++ pkg/inst/NEWS 2012-04-15 15:43:19 UTC (rev 660)
@@ -1,13 +1,21 @@
NEWS
+0.41-6
+ o The 'gompertz' example parameter transformations have been changed.
+ No longer are the names of the parameter vector changed in the transformation.
+ This change is not backward-compatible, but only codes using this example are affected.
+
+ o The 'euler.sir' and 'gillespie.sir' examples have been changed.
+ The transmission rate beta(t) is now the arithmetic sum of the seasonality basis functions.
+ Before, it was the geometric sum.
+ R0 is now given by the arithmetic average of the beta parameters divided by (gamma+mu).
+ This change is not backward-compatible, but only codes using these examples are affected.
+
0.41-5
o An experimental facility for constructing pomp objects with native C routines is now included.
- o The 'gompertz' example parameter transformations have been changed.
- No longer are the names of the parameter vector changed in the transformation.
-
0.41-4
o The 'blowflies', 'euler.sir', 'gillespie.sir', 'bbs', and 'ricker' data()-loadable examples have been changed to make the parameterization simpler and more natural.
- Although these changes are not backward compatible, they make for much simpler explanations.
+ This change is not backward-compatible, but only codes using this example are affected.
0.41-3
o In 'trajectory', all vectorfield and map evaluation is now done in C for speed.
@@ -24,8 +32,10 @@
o New arguments in 'mif', 'nlf', 'bsmc', 'pmcmc', 'probe-match', and 'traj-match' allow the estimation to be done on a transformed parameter space.
When 'transform=TRUE' in these commands ('transform.params=TRUE' for 'nlf'), estimation is performed on the transformed parameter space.
This is described and demonstrated in the 'intro_to_pomp' vignette.
- The data()-loadable examples have been re-implemented to make use of this facility.
+
+ o The data()-loadable examples have been re-implemented to make use of the above-mentioned facility.
Note that this new functionality makes it unnecessary to "un-transform" model parameters within the user-specified 'rprocess', 'dprocess', 'rmeasure', 'dmeasure', 'skeleton', and 'initializer' codes.
+ This change is not backward-compatible, but only codes using these data()-loadable example are affected.
o The Bayesian sequential Monte Carlo command 'bsmc' now returns not a list but an object of class 'bsmcd.pomp'.
An experimental 'plot' method for objects of this class now exists.
Modified: pkg/inst/data-R/sir.R
===================================================================
--- pkg/inst/data-R/sir.R 2012-04-14 21:19:02 UTC (rev 659)
+++ pkg/inst/data-R/sir.R 2012-04-15 15:43:19 UTC (rev 660)
@@ -43,14 +43,14 @@
coef(po) <- c(
gamma=26,mu=0.02,iota=0.01,
nbasis=3,degree=3,period=1,
- beta1=1200,beta2=1800,beta3=600,
+ beta1=400,beta2=480,beta3=320,
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
+ S.0=26/400,I.0=0.001,R.0=1-26/400
)
-simulate(po,nsim=1,seed=329348545L) -> euler.sir
+simulate(po,nsim=1,seed=329343545L) -> euler.sir
save(euler.sir,file="euler.sir.rda",compress="xz")
Modified: pkg/inst/examples/sir.c
===================================================================
--- pkg/inst/examples/sir.c 2012-04-14 21:19:02 UTC (rev 659)
+++ pkg/inst/examples/sir.c 2012-04-15 15:43:19 UTC (rev 660)
@@ -85,10 +85,9 @@
// compute transmission rate from seasonality
if (nbasis <= 0) return;
- (*eval_basis)(t,PERIOD,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
+ eval_basis(t,PERIOD,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
for (k = 0, beta = 0; k < nbasis; k++)
- beta += log(BETA[k])*seasonality[k];
- beta = exp(beta);
+ beta += BETA[k]*seasonality[k];
// compute the environmental stochasticity
dW = rgammawn(BETA_SD,dt);
@@ -136,10 +135,9 @@
// compute transmission rate from seasonality
if (nbasis <= 0) return;
- (*eval_basis)(t,PERIOD,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
+ eval_basis(t,PERIOD,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
for (k = 0, beta = 0; k < nbasis; k++)
- beta += log(BETA[k])*seasonality[k];
- beta = exp(beta);
+ beta += BETA[k]*seasonality[k];
// compute the transition rates
rate[0] = MU*POPSIZE; // birth into susceptible class
Modified: pkg/man/sir.Rd
===================================================================
--- pkg/man/sir.Rd 2012-04-14 21:19:02 UTC (rev 659)
+++ pkg/man/sir.Rd 2012-04-15 15:43:19 UTC (rev 660)
@@ -16,8 +16,12 @@
data(bbs)
}
\details{
+ This example is discussed extensively in the \dQuote{Introduction to \pkg{pomp}} and \dQuote{Advanced topics in \pkg{pomp}} vignettes.
+
The codes that construct these \code{pomp} objects can be found in the \dQuote{data-R} directory in the installed package.
Do \code{file.show(system.file("data-R/sir.R",package="pomp"))} to view these codes.
+ For the basic \code{rprocess}, \code{dmeasure}, \code{rmeasure}, and \code{skeleton} functions, these codes use compiled native routines built into the package's library.
+ View \dQuote{src/sir.c} in the package source or \code{file.show("examples/sir.c")} from an \R session to view these codes.
The boarding school influenza outbreak is described in Anonymous (1978).
}
Modified: pkg/src/sir.c
===================================================================
--- pkg/src/sir.c 2012-04-14 21:19:02 UTC (rev 659)
+++ pkg/src/sir.c 2012-04-15 15:43:19 UTC (rev 660)
@@ -132,8 +132,7 @@
if (nbasis <= 0) return;
periodic_bspline_basis_eval(t,PERIOD,deg,nbasis,&seasonality[0]);
for (k = 0, beta = 0; k < nbasis; k++)
- beta += seasonality[k]*log(BETA[k]);
- beta = exp(beta);
+ beta += seasonality[k]*BETA[k];
// test to make sure the parameters and state variable values are sane
if (!(R_FINITE(beta)) ||
@@ -190,8 +189,7 @@
if (nbasis <= 0) return;
periodic_bspline_basis_eval(t,PERIOD,deg,nbasis,&seasonality[0]);
for (k = 0, beta = 0; k < nbasis; k++)
- beta += seasonality[k]*log(BETA[k]);
- beta = exp(beta);
+ beta += seasonality[k]*BETA[k];
// compute the transition rates
rate[0] = MU*POPSIZE; // birth into susceptible class
@@ -250,8 +248,7 @@
case 3: // infection
periodic_bspline_basis_eval(t,PERIOD,deg,nbasis,&seasonality[0]);
for (k = 0, beta = 0; k < nbasis; k++)
- beta += seasonality[k]*log(BETA[k]);
- beta = exp(beta);
+ beta += seasonality[k]*BETA[k];
rate = (beta*INFD+IOTA)*SUSC/POPSIZE;
break;
case 4: // infected death
Modified: pkg/tests/gillespie.Rout.save
===================================================================
--- pkg/tests/gillespie.Rout.save 2012-04-14 21:19:02 UTC (rev 659)
+++ pkg/tests/gillespie.Rout.save 2012-04-15 15:43:19 UTC (rev 660)
@@ -122,13 +122,13 @@
>
> tail(as.data.frame(simulate(gillespie.sir,times=time(gsir),t0=timezero(gsir),seed=1165270654L)))
time reports S I R N cases
-100 1.903846 49 65867 1145 932868 999880 478
-101 1.923077 49 65570 1151 933153 999874 545
-102 1.942308 57 65302 1092 933464 999858 581
-103 1.961538 41 65048 1097 933683 999828 509
-104 1.980769 44 64737 1140 933910 999787 504
-105 2.000000 47 64491 1109 934157 999757 530
+100 1.903846 50 65088 1273 933471 999832 551
+101 1.923077 50 64777 1200 933837 999814 647
+102 1.942308 59 64442 1230 934114 999786 559
+103 1.961538 62 64060 1258 934433 999751 586
+104 1.980769 56 63775 1219 934732 999726 587
+105 2.000000 49 63461 1286 935000 999747 531
>
> proc.time()
user system elapsed
- 2.280 0.036 2.333
+ 2.248 0.044 2.308
Modified: pkg/tests/pfilter.Rout.save
===================================================================
--- pkg/tests/pfilter.Rout.save 2012-04-14 21:19:02 UTC (rev 659)
+++ pkg/tests/pfilter.Rout.save 2012-04-15 15:43:19 UTC (rev 660)
@@ -45,29 +45,25 @@
> data(euler.sir)
> pf <- pfilter(euler.sir,Np=100,seed=394343L)
> print(coef(pf))
- gamma mu iota nbasis degree period
-2.600000e+01 2.000000e-02 1.000000e-02 3.000000e+00 3.000000e+00 1.000000e+00
- beta1 beta2 beta3 beta.sd pop rho
-1.200000e+03 1.800000e+03 6.000000e+02 1.000000e-03 2.100000e+06 6.000000e-01
- S.0 I.0 R.0
-2.166667e-02 1.000000e-03 9.773333e-01
+ gamma mu iota nbasis degree period beta1 beta2
+2.60e+01 2.00e-02 1.00e-02 3.00e+00 3.00e+00 1.00e+00 4.00e+02 4.80e+02
+ beta3 beta.sd pop rho S.0 I.0 R.0
+3.20e+02 1.00e-03 2.10e+06 6.00e-01 6.50e-02 1.00e-03 9.35e-01
> print(pf$loglik,digits=4)
-[1] -892.7
+[1] -964.3
>
> p <- coef(euler.sir)
> euler.sir at params <- numeric(0)
> p["iota"] <- 1
> pf <- pfilter(euler.sir,params=p,Np=100,seed=394343L)
> print(coef(pf))
- gamma mu iota nbasis degree period
-2.600000e+01 2.000000e-02 1.000000e+00 3.000000e+00 3.000000e+00 1.000000e+00
- beta1 beta2 beta3 beta.sd pop rho
-1.200000e+03 1.800000e+03 6.000000e+02 1.000000e-03 2.100000e+06 6.000000e-01
- S.0 I.0 R.0
-2.166667e-02 1.000000e-03 9.773333e-01
+ gamma mu iota nbasis degree period beta1 beta2
+2.60e+01 2.00e-02 1.00e+00 3.00e+00 3.00e+00 1.00e+00 4.00e+02 4.80e+02
+ beta3 beta.sd pop rho S.0 I.0 R.0
+3.20e+02 1.00e-03 2.10e+06 6.00e-01 6.50e-02 1.00e-03 9.35e-01
> print(logLik(pf),digits=4)
-[1] -907.6
+[1] -964
>
> proc.time()
user system elapsed
- 5.772 0.048 5.858
+ 5.732 0.080 5.850
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2012-04-14 21:19:02 UTC (rev 659)
+++ pkg/tests/sir.R 2012-04-15 15:43:19 UTC (rev 660)
@@ -7,11 +7,11 @@
## some parameters
params <- c(
gamma=26,mu=0.02,iota=0.01,
- beta1=1200,beta2=1800,beta3=600,
+ beta1=400,beta2=480,beta3=320,
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
+ S.0=26/400,I.0=0.001,R.0=1-26/400
)
## set up the pomp object
@@ -26,13 +26,11 @@
rprocess=euler.sim(
delta.t=1/52/20,
step.fun=function(t,x,params,covars,delta.t,...) {
- params <- exp(params)
with(
as.list(c(x,params)),
{
- beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
- beta.var <- beta.sd^2
- dW <- rgamma(n=1,shape=delta.t/beta.var,scale=beta.var)
+ beta <- sum(c(beta1,beta2,beta3)*covars)
+ dW <- rgammawn(n=1,sigma=beta.sd,dt=delta.t)
foi <- (iota+beta*I*dW/delta.t)/pop
trans <- c(
rpois(n=1,lambda=mu*pop*delta.t),
@@ -60,12 +58,11 @@
),
dprocess=onestep.dens(
dens.fun=function(t1,t2,params,x1,x2,covars,...) {
- params <- exp(params)
with(
as.list(params),
{
dt <- t2-t1
- beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
+ beta <- sum(c(beta1,beta2,beta3)*covars)
beta.var <- beta.sd^2
dW <- x2['dW']
foi <- (iota+beta*x1["I"]*dW/dt)/pop
@@ -84,11 +81,10 @@
skeleton.type="vectorfield",
skeleton=function(x,t,params,covars,...) {
xdot <- rep(0,length(x))
- params <- exp(params)
with(
as.list(c(x,params)),
{
- beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
+ beta <- sum(c(beta1,beta2,beta3)*covars)
foi <- (iota+beta*I)/pop
terms <- c(
mu*pop,
@@ -108,10 +104,10 @@
}
)
},
-# measurement.model=reports~binom(size=cases,prob=exp(rho)),
+# measurement.model=reports~binom(size=cases,prob=rho),
rmeasure=function(x,t,params,covars,...){
with(
- as.list(c(x,exp(params))),
+ as.list(c(x,params)),
{
rep <- round(rnorm(n=1,mean=rho*cases,sd=sqrt(rho*(1-rho)*cases)))
if (rep<0) rep <- 0
@@ -121,7 +117,7 @@
},
dmeasure=function(y,x,t,params,log,covars,...){
with(
- as.list(c(x,exp(params))),
+ as.list(c(x,params)),
{
if (y > 0)
f <- diff(pnorm(q=y+c(-0.5,0.5),mean=rho*cases,sd=sqrt(rho*(1-rho)*cases),lower.tail=TRUE,log.p=FALSE))
@@ -132,9 +128,8 @@
)
},
initializer=function(params,t0,...){
- p <- exp(params)
with(
- as.list(p),
+ as.list(params),
{
fracs <- c(S.0,I.0,R.0)
x0 <- c(
@@ -154,7 +149,7 @@
set.seed(3049953)
## simulate from the model
tic <- Sys.time()
-x <- simulate(po,params=log(params),nsim=3)
+x <- simulate(po,params=params,nsim=3)
toc <- Sys.time()
print(toc-tic)
@@ -163,14 +158,14 @@
plot(x[[1]],variables=c("S","I","R","cases","W"))
t1 <- seq(0,4/52,by=1/52/25)
-X1 <- simulate(po,params=log(params),nsim=10,states=TRUE,obs=TRUE,times=t1)
+X1 <- simulate(po,params=params,nsim=10,states=TRUE,obs=TRUE,times=t1)
t2 <- seq(0,2,by=1/52)
-X2 <- simulate(po,params=log(params),nsim=1,states=TRUE,obs=TRUE,times=t2)
+X2 <- simulate(po,params=params,nsim=1,states=TRUE,obs=TRUE,times=t2)
t3 <- seq(0,20,by=1/52)
tic <- Sys.time()
-X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
+X3 <- trajectory(po,params=params,times=t3,hmax=1/52)
toc <- Sys.time()
print(toc-tic)
plot(t3,X3['I',1,],type='l')
@@ -179,12 +174,7 @@
po,
x=X1$states[,,31:40],
times=t1[31:40],
- params=matrix(
- log(params),
- nrow=length(params),
- ncol=10,
- dimnames=list(names(params),NULL)
- ),
+ params=parmat(params,nrep=10),
log=TRUE
)
print(apply(f1,1,sum),digits=4)
@@ -194,12 +184,7 @@
y=rbind(reports=X1$obs[,7,]),
x=X1$states,
times=t1,
- params=matrix(
- log(params),
- nrow=length(params),
- ncol=10,
- dimnames=list(names(params),NULL)
- ),
+ params=parmat(params,nrep=10),
log=TRUE
)
print(apply(g1,1,sum),digits=4)
@@ -208,7 +193,7 @@
po,
x=X2$states[,1,55:70,drop=FALSE],
t=t2[55:70],
- params=as.matrix(log(params))
+ params=params
)
print(h1[c("S","I","R"),,],digits=4)
@@ -236,12 +221,7 @@
y=rbind(reports=X1$obs[,7,]),
x=X1$states,
times=t1,
- params=matrix(
- coef(po),
- nrow=length(params)+3,
- ncol=10,
- dimnames=list(names(coef(po)),NULL)
- ),
+ params=parmat(coef(po),nrep=10),
log=TRUE
)
print(apply(g2,1,sum),digits=4)
@@ -250,7 +230,7 @@
po,
x=X2$states[,1,55:70,drop=FALSE],
t=t2[55:70],
- params=as.matrix(coef(po))
+ params=coef(po)
)
print(h2[c("S","I","R"),,],digits=4)
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2012-04-14 21:19:02 UTC (rev 659)
+++ pkg/tests/sir.Rout.save 2012-04-15 15:43:19 UTC (rev 660)
@@ -28,11 +28,11 @@
> ## some parameters
> params <- c(
+ gamma=26,mu=0.02,iota=0.01,
-+ beta1=1200,beta2=1800,beta3=600,
++ beta1=400,beta2=480,beta3=320,
+ 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
++ S.0=26/400,I.0=0.001,R.0=1-26/400
+ )
>
> ## set up the pomp object
@@ -47,13 +47,11 @@
+ rprocess=euler.sim(
+ delta.t=1/52/20,
+ step.fun=function(t,x,params,covars,delta.t,...) {
-+ params <- exp(params)
+ with(
+ as.list(c(x,params)),
+ {
-+ beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
-+ beta.var <- beta.sd^2
-+ dW <- rgamma(n=1,shape=delta.t/beta.var,scale=beta.var)
++ beta <- sum(c(beta1,beta2,beta3)*covars)
++ dW <- rgammawn(n=1,sigma=beta.sd,dt=delta.t)
+ foi <- (iota+beta*I*dW/delta.t)/pop
+ trans <- c(
+ rpois(n=1,lambda=mu*pop*delta.t),
@@ -81,12 +79,11 @@
+ ),
+ dprocess=onestep.dens(
+ dens.fun=function(t1,t2,params,x1,x2,covars,...) {
-+ params <- exp(params)
+ with(
+ as.list(params),
+ {
+ dt <- t2-t1
-+ beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
++ beta <- sum(c(beta1,beta2,beta3)*covars)
+ beta.var <- beta.sd^2
+ dW <- x2['dW']
+ foi <- (iota+beta*x1["I"]*dW/dt)/pop
@@ -105,11 +102,10 @@
+ skeleton.type="vectorfield",
+ skeleton=function(x,t,params,covars,...) {
+ xdot <- rep(0,length(x))
-+ params <- exp(params)
+ with(
+ as.list(c(x,params)),
+ {
-+ beta <- exp(sum(log(c(beta1,beta2,beta3))*covars))
++ beta <- sum(c(beta1,beta2,beta3)*covars)
+ foi <- (iota+beta*I)/pop
+ terms <- c(
+ mu*pop,
@@ -129,10 +125,10 @@
+ }
+ )
+ },
-+ # measurement.model=reports~binom(size=cases,prob=exp(rho)),
++ # measurement.model=reports~binom(size=cases,prob=rho),
+ rmeasure=function(x,t,params,covars,...){
+ with(
-+ as.list(c(x,exp(params))),
++ as.list(c(x,params)),
+ {
+ rep <- round(rnorm(n=1,mean=rho*cases,sd=sqrt(rho*(1-rho)*cases)))
+ if (rep<0) rep <- 0
@@ -142,7 +138,7 @@
+ },
+ dmeasure=function(y,x,t,params,log,covars,...){
+ with(
-+ as.list(c(x,exp(params))),
++ as.list(c(x,params)),
+ {
+ if (y > 0)
+ f <- diff(pnorm(q=y+c(-0.5,0.5),mean=rho*cases,sd=sqrt(rho*(1-rho)*cases),lower.tail=TRUE,log.p=FALSE))
@@ -153,9 +149,8 @@
+ )
+ },
+ initializer=function(params,t0,...){
-+ p <- exp(params)
+ with(
-+ as.list(p),
++ as.list(params),
+ {
+ fracs <- c(S.0,I.0,R.0)
+ x0 <- c(
@@ -394,7 +389,7 @@
zeronames = zeronames, tcovar = tcovar, covar = covar,
args = pairlist(...))
}
-<environment: 0x37f5b80>
+<environment: 0x25543e0>
process model density, dprocess =
function (x, times, params, ..., statenames = character(0), paramnames = character(0),
covarnames = character(0), tcovar, covar, log = FALSE)
@@ -404,11 +399,11 @@
covarnames = covarnames, tcovar = tcovar, covar = covar,
log = log, args = pairlist(...))
}
-<environment: 0x3921ab8>
+<environment: 0x24e7238>
measurement model simulator, rmeasure =
function (x, t, params, covars, ...)
{
- with(as.list(c(x, exp(params))), {
+ with(as.list(c(x, params)), {
rep <- round(rnorm(n = 1, mean = rho * cases, sd = sqrt(rho *
(1 - rho) * cases)))
if (rep < 0)
@@ -419,7 +414,7 @@
measurement model density, dmeasure =
function (y, x, t, params, log, covars, ...)
{
- with(as.list(c(x, exp(params))), {
+ with(as.list(c(x, params)), {
if (y > 0)
f <- diff(pnorm(q = y + c(-0.5, 0.5), mean = rho *
cases, sd = sqrt(rho * (1 - rho) * cases), lower.tail = TRUE,
@@ -435,9 +430,8 @@
function (x, t, params, covars, ...)
{
xdot <- rep(0, length(x))
- params <- exp(params)
with(as.list(c(x, params)), {
- beta <- exp(sum(log(c(beta1, beta2, beta3)) * covars))
+ beta <- sum(c(beta1, beta2, beta3) * covars)
foi <- (iota + beta * I)/pop
terms <- c(mu * pop, foi * S, mu * S, gamma * I, mu *
I, mu * R)
@@ -449,8 +443,7 @@
initializer =
function (params, t0, ...)
{
- p <- exp(params)
- with(as.list(p), {
+ with(as.list(params), {
fracs <- c(S.0, I.0, R.0)
x0 <- c(round(pop * fracs/sum(fracs)), rep(0, 9))
names(x0) <- c("S", "I", "R", "cases", "W", "B", "SI",
@@ -467,75 +460,65 @@
> set.seed(3049953)
> ## simulate from the model
> tic <- Sys.time()
-> x <- simulate(po,params=log(params),nsim=3)
+> x <- simulate(po,params=params,nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 1.090771 secs
+Time difference of 1.024599 secs
>
> pdf(file='sir.pdf')
>
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t1 <- seq(0,4/52,by=1/52/25)
-> X1 <- simulate(po,params=log(params),nsim=10,states=TRUE,obs=TRUE,times=t1)
+> X1 <- simulate(po,params=params,nsim=10,states=TRUE,obs=TRUE,times=t1)
>
> t2 <- seq(0,2,by=1/52)
-> X2 <- simulate(po,params=log(params),nsim=1,states=TRUE,obs=TRUE,times=t2)
+> X2 <- simulate(po,params=params,nsim=1,states=TRUE,obs=TRUE,times=t2)
>
> t3 <- seq(0,20,by=1/52)
> tic <- Sys.time()
-> X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
+> X3 <- trajectory(po,params=params,times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 1.077192 secs
+Time difference of 0.7358234 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
+ po,
+ x=X1$states[,,31:40],
+ times=t1[31:40],
-+ params=matrix(
-+ log(params),
-+ nrow=length(params),
-+ ncol=10,
-+ dimnames=list(names(params),NULL)
-+ ),
++ params=parmat(params,nrep=10),
+ log=TRUE
+ )
> print(apply(f1,1,sum),digits=4)
- [1] -41.29 -40.62 -58.44 -40.74 -48.97 -46.98 -40.59 -44.77 -49.61 -42.44
+ [1] -47.60 -63.21 -50.66 -48.84 -48.28 -52.25 -45.21 -45.75 -51.05 -49.98
>
> g1 <- dmeasure(
+ po,
+ y=rbind(reports=X1$obs[,7,]),
+ x=X1$states,
+ times=t1,
-+ params=matrix(
-+ log(params),
-+ nrow=length(params),
-+ ncol=10,
-+ dimnames=list(names(params),NULL)
-+ ),
++ params=parmat(params,nrep=10),
+ log=TRUE
+ )
> print(apply(g1,1,sum),digits=4)
- [1] -363.9 -400.0 -418.2 -387.1 -404.6 -417.2 -251.7 -440.1 -404.5 -416.1
+ [1] -406.8 -396.6 -Inf -383.8 -447.9 -380.6 -255.1 -382.1 -460.4 -Inf
>
> h1 <- skeleton(
+ po,
+ x=X2$states[,1,55:70,drop=FALSE],
+ t=t2[55:70],
-+ params=as.matrix(log(params))
++ params=params
+ )
> print(h1[c("S","I","R"),,],digits=4)
- [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
-S 39608 39770.1 39222.4 39000 38470.0 37635 36130 34844 33808 30538
-I 243 246.8 415.5 519 722.8 1047 1628 2172 2649 4006
-R -39852 -40019.5 -39639.0 -39520 -39191.5 -38682 -37758 -37016 -36456 -34544
- [,11] [,12] [,13] [,14] [,15] [,16]
-S 27082 20266 13528 2809 -8808 -26739.3
-I 5473 8300 11076 15339 19697 25951.9
-R -32554 -28565 -24604 -18148 -10889 787.2
+ [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
+S 32161.6 31582.4 32142.21 31758.36 31639.0 32006 31953.5 31666
+I -318.4 -218.9 -88.93 21.52 128.7 218 309.2 403
+R -31855.2 -31375.0 -32064.98 -31792.40 -31780.0 -32236 -32276.8 -32083
+ [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
+S 31017.5 29782.5 29443.9 28916.7 27997.1 27494 26434 25448
+I 518.5 675.3 768.6 868.8 996.8 1074 1187 1278
+R -31549.7 -30472.3 -30226.0 -29798.8 -29007.2 -28580 -27634 -26740
>
> ## now repeat using the compiled native codes built into the package
> data(euler.sir)
@@ -547,7 +530,7 @@
> x <- simulate(po,nsim=100)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 1.203072 secs
+Time difference of 1.200051 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
@@ -555,7 +538,7 @@
> X4 <- trajectory(po,times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.01984811 secs
+Time difference of 0.01725101 secs
> plot(t3,X4['I',1,],type='l')
>
> g2 <- dmeasure(
@@ -563,60 +546,55 @@
+ y=rbind(reports=X1$obs[,7,]),
+ x=X1$states,
+ times=t1,
-+ params=matrix(
-+ coef(po),
-+ nrow=length(params)+3,
-+ ncol=10,
-+ dimnames=list(names(coef(po)),NULL)
-+ ),
++ params=parmat(coef(po),nrep=10),
+ log=TRUE
+ )
> print(apply(g2,1,sum),digits=4)
- [1] -363.9 -400.0 -418.2 -387.1 -404.6 -417.2 -251.7 -440.1 -404.5 -416.1
+ [1] -406.8 -396.6 -Inf -383.8 -447.9 -380.6 -255.1 -382.1 -460.4 -Inf
>
> h2 <- skeleton(
+ po,
+ x=X2$states[,1,55:70,drop=FALSE],
+ t=t2[55:70],
-+ params=as.matrix(coef(po))
++ params=coef(po)
+ )
> print(h2[c("S","I","R"),,],digits=4)
- [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
-S 39608 39770.1 39222.4 39000 38470.0 37635 36130 34844 33808 30538
-I 243 246.8 415.5 519 722.8 1047 1628 2172 2649 4006
-R -39852 -40019.5 -39639.0 -39520 -39191.5 -38682 -37758 -37016 -36456 -34544
- [,11] [,12] [,13] [,14] [,15] [,16]
-S 27082 20266 13528 2809 -8808 -26739.3
-I 5473 8300 11076 15339 19697 25951.9
-R -32554 -28565 -24604 -18148 -10889 787.2
+ [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
+S 32161.6 31582.4 32142.21 31758.36 31639.0 32006 31953.5 31666
+I -318.4 -218.9 -88.93 21.52 128.7 218 309.2 403
+R -31855.2 -31375.0 -32064.98 -31792.40 -31780.0 -32236 -32276.8 -32083
+ [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
+S 31017.5 29782.5 29443.9 28916.7 27997.1 27494 26434 25448
+I 518.5 675.3 768.6 868.8 996.8 1074 1187 1278
+R -31549.7 -30472.3 -30226.0 -29798.8 -29007.2 -28580 -27634 -26740
>
> print(max(abs(g2-g1),na.rm=T),digits=4)
[1] 7.105e-15
> print(max(abs(h2-h1),na.rm=T),digits=4)
-[1] 9.459e-11
+[1] 7.276e-12
>
> states(po)[,1:2]
- [,1] [,2]
-S 4.531800e+04 4.511300e+04
-I 2.038000e+03 2.020000e+03
-R 2.052647e+06 2.052886e+06
-cases 1.045000e+03 9.960000e+02
-W -2.334331e-01 -2.382536e-02
+ [,1] [,2]
+S 1.361010e+05 1.357700e+05
+I 2.064000e+03 2.065000e+03
+R 1.961799e+06 1.962133e+06
+cases 1.062000e+03 1.060000e+03
+W 6.596952e-02 9.166551e-02
> time(po) <- seq(0,1,by=1/52)
> states(po)[,1:3]
- [,1] [,2] [,3]
-S NA 4.531800e+04 4.511300e+04
-I NA 2.038000e+03 2.020000e+03
-R NA 2.052647e+06 2.052886e+06
-cases NA 1.045000e+03 9.960000e+02
-W NA -2.334331e-01 -2.382536e-02
+ [,1] [,2] [,3]
+S NA 1.361010e+05 1.357700e+05
+I NA 2.064000e+03 2.065000e+03
+R NA 1.961799e+06 1.962133e+06
+cases NA 1.062000e+03 1.060000e+03
+W NA 6.596952e-02 9.166551e-02
> states(simulate(po))[,1:3]
- [,1] [,2] [,3]
-S 45500 4.526000e+04 4.506900e+04
-I 2100 2.091000e+03 2.141000e+03
-R 2052400 2.052645e+06 2.052811e+06
-cases 0 1.013000e+03 9.660000e+02
-W 0 9.854644e-03 -1.055328e-01
+ [,1] [,2] [,3]
+S 136364 1.360240e+05 1.357070e+05
+I 2098 2.164000e+03 2.218000e+03
+R 1961538 1.961821e+06 1.962178e+06
+cases 0 9.990000e+02 1.076000e+03
+W 0 -1.381354e-02 3.473754e-01
>
> po <- window(euler.sir,start=1,end=2)
> plot(simulate(po))
@@ -632,4 +610,4 @@
>
> proc.time()
user system elapsed
- 4.304 0.056 4.595
+ 3.920 0.036 4.195
More information about the pomp-commits
mailing list