[Pomp-commits] r104 - in pkg: R data man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Apr 24 22:58:10 CEST 2009
Author: kingaa
Date: 2009-04-24 22:58:10 +0200 (Fri, 24 Apr 2009)
New Revision: 104
Added:
pkg/data/euler.sir.R
pkg/man/euler.sir.Rd
Modified:
pkg/R/trajectory-pomp.R
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
put 'euler.sir' in as an example pomp object obtainable via 'data(euler.sir)'
make 'trajectory' a bit more user-friendly
Modified: pkg/R/trajectory-pomp.R
===================================================================
--- pkg/R/trajectory-pomp.R 2009-04-24 17:37:41 UTC (rev 103)
+++ pkg/R/trajectory-pomp.R 2009-04-24 20:58:10 UTC (rev 104)
@@ -2,12 +2,37 @@
"trajectory",
"pomp",
function (object, params, times, ...) {
+ if (missing(params)) {
+ params <- coef(object)
+ if (length(params)==0) {
+ stop("trajectory error: ",sQuote("params")," must be supplied",call.=FALSE)
+ }
+ }
+ nrep <- NCOL(params)
+ if (is.null(dim(params))) {
+ params <- matrix(
+ params,
+ nrow=length(params),
+ ncol=nrep,
+ dimnames=list(
+ names(params),
+ NULL
+ )
+ )
+ }
+ paramnames <- rownames(params)
+ if (is.null(paramnames))
+ stop("pfilter error: ",sQuote("params")," must have rownames",call.=FALSE)
+
+ if (missing(times))
+ times <- time(object,t0=TRUE)
+
params <- as.matrix(params)
if (missing(times))
times <- time(object,t0=TRUE)
x0 <- init.state(object,params=params,t0=times[1])
x <- array(
- dim=c(nrow(x0),ncol(x0),length(times)),
+ dim=c(nrow(x0),nrep,length(times)),
dimnames=list(rownames(x0),NULL,NULL)
)
switch(
@@ -24,7 +49,7 @@
}
},
vectorfield={ # integrate the vectorfield
- for (j in 1:ncol(params)) {
+ for (j in 1:nrep) {
X <- try(
lsoda(
y=x0[,j],
Added: pkg/data/euler.sir.R
===================================================================
--- pkg/data/euler.sir.R (rev 0)
+++ pkg/data/euler.sir.R 2009-04-24 20:58:10 UTC (rev 104)
@@ -0,0 +1,54 @@
+require(pomp)
+
+euler.sir <- simulate(
+ 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","B","dW"),
+ 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,
+ dens.fun="sir_euler_density",
+ dprocess=euler.density,
+ 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","B",
+ "SI","SD","IR","ID","RD","dW"
+ )
+ fracs <- p[paste(comp.names,"0",sep=".")]
+ x0 <- numeric(length(snames))
+ names(x0) <- snames
+ x0[comp.names] <- round(p['pop']*fracs/sum(fracs))
+ x0
+ }
+ ),
+ params=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
+ )
+ ),
+ nsim=1,
+ seed=329348545L
+ )
+euler.sir <- euler.sir[[1]]
Added: pkg/man/euler.sir.Rd
===================================================================
--- pkg/man/euler.sir.Rd (rev 0)
+++ pkg/man/euler.sir.Rd 2009-04-24 20:58:10 UTC (rev 104)
@@ -0,0 +1,18 @@
+\name{euler.sir}
+\alias{euler.sir}
+\docType{data}
+\title{Seasonal SIR model implemented as an Euler-multinomial model}
+\description{
+ \code{euler.sir} is a \code{pomp} object encoding a simple seasonal SIR model.
+}
+\usage{data(euler.sir)}
+\details{
+}
+\examples{
+data(euler.sir)
+plot(euler.sir)
+x <- simulate(euler.sir,nsim=10,seed=20348585)
+plot(x[[1]])
+}
+\seealso{\code{\link{pomp-class}} and the vignettes}
+\keyword{datasets}
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2009-04-24 17:37:41 UTC (rev 103)
+++ pkg/tests/sir.R 2009-04-24 20:58:10 UTC (rev 104)
@@ -184,59 +184,25 @@
)
print(h1[c("S","I","R"),,],digits=4)
-po <- pomp(
- 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","B","dW"),
- paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho"),
- covarnames=c("seas1"),
- zeronames=c("cases"),
- step.fun="sir_euler_simulator",
- rprocess=euler.simulate,
- dens.fun="sir_euler_density",
- dprocess=euler.density,
- skeleton.vectorfield="sir_ODE",
- rmeasure="binom_rmeasure",
- dmeasure="binom_dmeasure",
- PACKAGE="pomp",
- initializer=function(params,t0,...){
- p <- exp(params)
- with(
- as.list(p),
- {
- fracs <- c(S.0,I.0,R.0)
- x0 <- c(
- round(pop*fracs/sum(fracs)), # make sure the three compartments sum to 'pop' initially
- rep(0,9) # zeros for 'cases', 'W', and the transition numbers
- )
- names(x0) <- c("S","I","R","cases","W","B","SI","SD","IR","ID","RD","dW")
- x0
- }
- )
- }
- )
+data(euler.sir)
set.seed(3049953)
## simulate from the model
tic <- Sys.time()
-x <- simulate(po,params=log(params),nsim=3)
+x <- simulate(euler.sir,nsim=3)
toc <- Sys.time()
print(toc-tic)
plot(x[[1]],variables=c("S","I","R","cases","W"))
t3 <- seq(0,20,by=1/52)
tic <- Sys.time()
-X4 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
+X4 <- trajectory(euler.sir,times=t3,hmax=1/52)
toc <- Sys.time()
print(toc-tic)
plot(t3,X4['I',1,],type='l')
f2 <- dprocess(
- po,
+ euler.sir,
x=X1$states[,,31:40],
times=t1[31:40],
params=matrix(
@@ -250,7 +216,7 @@
print(apply(f2,1,sum),digits=4)
g2 <- dmeasure(
- po,
+ euler.sir,
y=rbind(measles=X1$obs[,7,]),
x=X1$states,
times=t1,
@@ -265,7 +231,7 @@
print(apply(g2,1,sum),digits=4)
h2 <- skeleton(
- po,
+ euler.sir,
x=X2$states[,1,55:70,drop=FALSE],
t=t2[55:70],
params=as.matrix(log(params))
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2009-04-24 17:37:41 UTC (rev 103)
+++ pkg/tests/sir.Rout.save 2009-04-24 20:58:10 UTC (rev 104)
@@ -148,7 +148,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 4.017606 secs
+Time difference of 4.005935 secs
>
> pdf(file='sir.pdf')
>
@@ -165,7 +165,7 @@
> X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 9.691413 secs
+Time difference of 9.746132 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
@@ -215,61 +215,27 @@
I 6165 8975 11503 15834.2 21418 27852 34849
R -30746 -26601 -22900 -16156.5 -6714 5662 22453
>
-> po <- pomp(
-+ 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","B","dW"),
-+ paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho"),
-+ covarnames=c("seas1"),
-+ zeronames=c("cases"),
-+ step.fun="sir_euler_simulator",
-+ rprocess=euler.simulate,
-+ dens.fun="sir_euler_density",
-+ dprocess=euler.density,
-+ skeleton.vectorfield="sir_ODE",
-+ rmeasure="binom_rmeasure",
-+ dmeasure="binom_dmeasure",
-+ PACKAGE="pomp",
-+ initializer=function(params,t0,...){
-+ p <- exp(params)
-+ with(
-+ as.list(p),
-+ {
-+ fracs <- c(S.0,I.0,R.0)
-+ x0 <- c(
-+ round(pop*fracs/sum(fracs)), # make sure the three compartments sum to 'pop' initially
-+ rep(0,9) # zeros for 'cases', 'W', and the transition numbers
-+ )
-+ names(x0) <- c("S","I","R","cases","W","B","SI","SD","IR","ID","RD","dW")
-+ x0
-+ }
-+ )
-+ }
-+ )
+> data(euler.sir)
>
> set.seed(3049953)
> ## simulate from the model
> tic <- Sys.time()
-> x <- simulate(po,params=log(params),nsim=3)
+> x <- simulate(euler.sir,nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.2553811 secs
+Time difference of 0.2860382 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
> tic <- Sys.time()
-> X4 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
+> X4 <- trajectory(euler.sir,times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 7.702386 secs
+Time difference of 7.802866 secs
> plot(t3,X4['I',1,],type='l')
>
> f2 <- dprocess(
-+ po,
++ euler.sir,
+ x=X1$states[,,31:40],
+ times=t1[31:40],
+ params=matrix(
@@ -284,7 +250,7 @@
[1] -51.57 -43.19 -50.14 -39.65 -41.63 -49.34 -39.60 -48.19 -42.60 -48.26
>
> g2 <- dmeasure(
-+ po,
++ euler.sir,
+ y=rbind(measles=X1$obs[,7,]),
+ x=X1$states,
+ times=t1,
@@ -300,7 +266,7 @@
[1] -Inf -Inf -Inf -Inf -Inf -Inf -263.6 -408.1 -Inf -468.1
>
> h2 <- skeleton(
-+ po,
++ euler.sir,
+ x=X2$states[,1,55:70,drop=FALSE],
+ t=t2[55:70],
+ params=as.matrix(log(params))
More information about the pomp-commits
mailing list