[Pomp-commits] r227 - pkg/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 10 15:53:46 CEST 2010
Author: kingaa
Date: 2010-05-10 15:53:46 +0200 (Mon, 10 May 2010)
New Revision: 227
Added:
pkg/tests/gillespie.R
pkg/tests/gillespie.Rout.save
Log:
- tests for new gillespie simulator
Added: pkg/tests/gillespie.R
===================================================================
--- pkg/tests/gillespie.R (rev 0)
+++ pkg/tests/gillespie.R 2010-05-10 13:53:46 UTC (rev 227)
@@ -0,0 +1,88 @@
+library(pomp)
+
+params <- c(
+ nu=log(1/70),
+ mu=log(1/70),
+ beta1=log(330),
+ beta2=log(410),
+ beta3=log(490),
+ gamma=log(24),
+ iota=log(0.1),
+ S.0=log(0.05),
+ I.0=log(1e-4),
+ R.0=log(0.95),
+ N.0=1000000
+ )
+
+seasonality <- data.frame(
+ time=seq(0,10,by=1/52/10),
+ seas=periodic.bspline.basis(seq(0,10,by=1/52/10),nbasis=3,degree=3,period=1)
+ )
+simulate(
+ pomp(
+ data=data.frame(
+ time=seq(from=0,to=2,by=1/52),
+ reports=NA
+ ),
+ times="time",
+ t0=0,
+ covar=seasonality,
+ tcovar="time",
+ rprocess=gillespie.sim(
+ rate.fun=function(j, x, t, params, covars, ...) {
+ switch(
+ j,
+ exp(params["nu"])*x["N"], # birth
+ exp(params["mu"])*x["S"], # susceptible death
+ { # infection
+ beta <- exp(params[c("beta1","beta2","beta3")])%*%covars[c("seas.1","seas.2","seas.3")]
+ (beta*x["I"]+exp(params["iota"]))*x["S"]/x["N"]
+ },
+ exp(params["mu"])*x["I"], # infected death
+ exp(params["gamma"])*x["I"], # recovery
+ exp(params["mu"])*x["R"], # recovered death
+ stop("unrecognized event ",j)
+ )
+ },
+ v=cbind(
+ birth=c(1,0,0,1,0),
+ sdeath=c(-1,0,0,-1,0),
+ infection=c(-1,1,0,0,0),
+ ideath=c(0,-1,0,-1,0),
+ recovery=c(0,-1,1,0,1),
+ rdeath=c(0,0,-1,-1,0)
+ ),
+ d=cbind(
+ birth=c(0,0,0,1,0),
+ sdeath=c(1,0,0,0,0),
+ infection=c(1,1,0,1,0),
+ ideath=c(0,1,0,0,0),
+ recovery=c(0,1,0,0,0),
+ rdeath=c(0,0,1,0,0)
+ )
+ ),
+ zeronames=c("cases"),
+ measurement.model=reports~binom(size=cases,prob=0.1),
+ initializer=function(params, t0, ...){
+ comp.names <- c("S","I","R")
+ icnames <- paste(comp.names,"0",sep=".")
+ snames <- c("S","I","R","N","cases")
+ fracs <- exp(params[icnames])
+ x0 <- numeric(length(snames))
+ names(x0) <- snames
+ x0["N"] <- params["N.0"]
+ x0[comp.names] <- round(params['N.0']*fracs/sum(fracs))
+ x0
+ }
+ ),
+ params=params,
+ nsim=1,
+ seed=1165270654L
+ ) -> gsir
+
+
+tail(as(gsir,"data.frame"))
+
+data(gillespie.sir)
+
+tail(as(simulate(gillespie.sir,times=time(gsir,t0=T),seed=1165270654L),"data.frame"))
Added: pkg/tests/gillespie.Rout.save
===================================================================
--- pkg/tests/gillespie.Rout.save (rev 0)
+++ pkg/tests/gillespie.Rout.save 2010-05-10 13:53:46 UTC (rev 227)
@@ -0,0 +1,123 @@
+
+R version 2.10.1 (2009-12-14)
+Copyright (C) 2009 The R Foundation for Statistical Computing
+ISBN 3-900051-07-0
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+> library(pomp)
+Loading required package: deSolve
+Loading required package: subplex
+Loading required package: mvtnorm
+>
+> params <- c(
++ nu=log(1/70),
++ mu=log(1/70),
++ beta1=log(330),
++ beta2=log(410),
++ beta3=log(490),
++ gamma=log(24),
++ iota=log(0.1),
++ S.0=log(0.05),
++ I.0=log(1e-4),
++ R.0=log(0.95),
++ N.0=1000000
++ )
+>
+> seasonality <- data.frame(
++ time=seq(0,10,by=1/52/10),
++ seas=periodic.bspline.basis(seq(0,10,by=1/52/10),nbasis=3,degree=3,period=1)
++ )
+> simulate(
++ pomp(
++ data=data.frame(
++ time=seq(from=0,to=2,by=1/52),
++ reports=NA
++ ),
++ times="time",
++ t0=0,
++ covar=seasonality,
++ tcovar="time",
++ rprocess=gillespie.sim(
++ rate.fun=function(j, x, t, params, covars, ...) {
++ switch(
++ j,
++ exp(params["nu"])*x["N"], # birth
++ exp(params["mu"])*x["S"], # susceptible death
++ { # infection
++ beta <- exp(params[c("beta1","beta2","beta3")])%*%covars[c("seas.1","seas.2","seas.3")]
++ (beta*x["I"]+exp(params["iota"]))*x["S"]/x["N"]
++ },
++ exp(params["mu"])*x["I"], # infected death
++ exp(params["gamma"])*x["I"], # recovery
++ exp(params["mu"])*x["R"], # recovered death
++ stop("unrecognized event ",j)
++ )
++ },
++ v=cbind(
++ birth=c(1,0,0,1,0),
++ sdeath=c(-1,0,0,-1,0),
++ infection=c(-1,1,0,0,0),
++ ideath=c(0,-1,0,-1,0),
++ recovery=c(0,-1,1,0,1),
++ rdeath=c(0,0,-1,-1,0)
++ ),
++ d=cbind(
++ birth=c(0,0,0,1,0),
++ sdeath=c(1,0,0,0,0),
++ infection=c(1,1,0,1,0),
++ ideath=c(0,1,0,0,0),
++ recovery=c(0,1,0,0,0),
++ rdeath=c(0,0,1,0,0)
++ )
++ ),
++ zeronames=c("cases"),
++ measurement.model=reports~binom(size=cases,prob=0.1),
++ initializer=function(params, t0, ...){
++ comp.names <- c("S","I","R")
++ icnames <- paste(comp.names,"0",sep=".")
++ snames <- c("S","I","R","N","cases")
++ fracs <- exp(params[icnames])
++ x0 <- numeric(length(snames))
++ names(x0) <- snames
++ x0["N"] <- params["N.0"]
++ x0[comp.names] <- round(params['N.0']*fracs/sum(fracs))
++ x0
++ }
++ ),
++ params=params,
++ nsim=1,
++ seed=1165270654L
++ ) -> gsir
+>
+>
+> tail(as(gsir,"data.frame"))
+ time reports S I R N cases
+100 1.903846 41 65075 1271 933488 999834 555
+101 1.923077 68 64761 1210 933849 999820 639
+102 1.942308 62 64434 1225 934136 999795 570
+103 1.961538 49 64055 1258 934441 999754 575
+104 1.980769 53 63765 1225 934740 999730 588
+105 2.000000 60 63451 1290 935012 999753 531
+>
+> data(gillespie.sir)
+>
+> tail(as(simulate(gillespie.sir,times=time(gsir,t0=T),seed=1165270654L),"data.frame"))
+ time reports S I R N cases
+100 1.903846 52 65867 1143 932864 999874 479
+101 1.923077 59 65570 1150 933149 999869 545
+102 1.942308 46 65300 1090 933462 999852 583
+103 1.961538 52 65047 1098 933677 999822 505
+104 1.980769 45 64735 1142 933909 999786 506
+105 2.000000 52 64487 1110 934156 999753 532
+>
More information about the pomp-commits
mailing list