[Pomp-commits] r1053 - in pkg: . pompExamples pompExamples/inst/examples pompExamples/vignettes

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri Jan 23 15:50:30 CET 2015 

Author: kingaa
Date: 2015-01-23 15:50:29 +0100 (Fri, 23 Jan 2015)
New Revision: 1053

Added:
   pkg/pompExamples/inst/examples/bsflu.R
   pkg/pompExamples/inst/examples/bsflu3.R
   pkg/pompExamples/vignettes/Makefile
   pkg/pompExamples/vignettes/bsflu-mf1.rds
   pkg/pompExamples/vignettes/bsflu-mf3.rds
   pkg/pompExamples/vignettes/bsflu-tm1.rds
   pkg/pompExamples/vignettes/bsflu-tm3.rds
   pkg/pompExamples/vignettes/bsflu.Rmd
   pkg/pompExamples/vignettes/bsflu.html
Modified:
   pkg/Makefile
   pkg/pompExamples/DESCRIPTION
Log:
- add boarding-school flu vignette to the pompExamples package

Modified: pkg/Makefile
===================================================================
--- pkg/Makefile	2015-01-22 16:10:38 UTC (rev 1052)
+++ pkg/Makefile	2015-01-23 14:50:29 UTC (rev 1053)
@@ -45,6 +45,9 @@
 	cp pomp.pdf ../www/vignettes
 	$(RCMD) Rdconv -t html pomp/inst/NEWS.Rd -o ../www/content/NEWS.html
 
+pompExamples.vignettes: pompExamples.install
+	(cd pompExamples/vignettes; make)
+
 %.data: %.install
 	cd $*/inst/data-R; make
 

Modified: pkg/pompExamples/DESCRIPTION
===================================================================
--- pkg/pompExamples/DESCRIPTION	2015-01-22 16:10:38 UTC (rev 1052)
+++ pkg/pompExamples/DESCRIPTION	2015-01-23 14:50:29 UTC (rev 1053)
@@ -1,8 +1,8 @@
 Package: pompExamples
 Type: Package
 Title: Additional pomp examples
-Version: 0.25-4
-Date: 2015-01-05
+Version: 0.26-1
+Date: 2015-01-23
 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"),
@@ -16,9 +16,10 @@
 	  person(given=c("Helen"),family="Wearing",role=c("ctb")))
 URL: http://pomp.r-forge.r-project.org
 Description: More 'pomp' examples.
-Depends: R(>= 3.0.0), stats, graphics, pomp(>= 0.56-1)
-Suggests: plyr, reshape2
+Depends: R(>= 3.0.0), stats, graphics, pomp(>= 0.58-6)
+Suggests: plyr, reshape2, knitr, ggplot2
 License: GPL (>= 2)
 LazyData: false
 BuildVignettes: true
+VignetteBuilder: knitr
 Collate: aaa.R pertussis.R

Added: pkg/pompExamples/inst/examples/bsflu.R
===================================================================
--- pkg/pompExamples/inst/examples/bsflu.R	                        (rev 0)
+++ pkg/pompExamples/inst/examples/bsflu.R	2015-01-23 14:50:29 UTC (rev 1053)
@@ -0,0 +1,109 @@
+## data read from graph in Anonymous (1978) 'Influenza in a boarding school' Brit. Med. J. 1:578.
+## cases are recorded with error of +/- 1 case
+## 763 boys were at risk, 512 boys spent time away from class
+flu <- read.csv(text="
+date,confined,convalescent
+1978-01-22,1,0
+1978-01-23,6,0
+1978-01-24,26,0
+1978-01-25,73,1
+1978-01-26,222,8
+1978-01-27,293,16
+1978-01-28,258,99
+1978-01-29,236,160
+1978-01-30,191,173
+1978-01-31,124,162
+1978-02-01,69,150
+1978-02-02,26,89
+1978-02-03,11,44
+1978-02-04,4,22
+",colClasses=c(date='Date'))
+flu$day <- flu$date-min(flu$date)+1
+units(flu$day) <- "days"
+flu$day <- as.numeric(flu$day)
+
+partrans <- "
+ TBeta = exp(Beta);
+ Tinf_pd = exp(inf_pd);
+ Trho = expit(rho);
+ Tsfrac = expit(sfrac);
+"
+
+paruntrans <- "
+ TBeta = log(Beta);
+ Tinf_pd = log(inf_pd);
+ Trho = logit(rho);
+ Tsfrac = logit(sfrac);
+"
+
+dmeas <- "
+  lik = dpois(confined,rho*R+1e-6,give_log);
+"
+
+rmeas <- "
+  confined = rpois(rho*R+1e-6);
+  convalescent = rpois(rho*C);
+"
+
+stochsim <- "
+  double t1 = rbinom(S,1-exp(-Beta*I*dt));
+  double t2 = rbinom(I,1-exp(-dt/inf_pd));
+  double t3 = rbinom(R,1-exp(-dt/conf_pd));
+  double t4 = rbinom(C,1-exp(-dt/conv_pd));
+  S -= t1;
+  I += t1 - t2;
+  R += t2 - t3;
+  C += t3 - t4;
+"
+
+skel <- "
+  double dt = 1.0/24.0;
+  double t1 = S*(1-exp(-Beta*I*dt));
+  double t2 = I*(1-exp(-dt/inf_pd));
+  double t3 = R*(1-exp(-dt/conf_pd));
+  double t4 = C*(1-exp(-dt/conv_pd));
+  DS = S - t1;
+  DI = I + t1 - t2;
+  DR = R + t2 - t3;
+  DC = C + t3 - t4;
+"
+
+pomp(
+     data=flu[c("day","confined","convalescent")],
+     times="day",
+     t0=0,
+     params=c(
+       Beta=0.004,
+       inf.pd=0.7,
+       conf.pd=sum(flu$confined)/512,
+       conv.pd=sum(flu$convalescent)/512,
+       rho=0.9,
+       sfrac=762/763
+       ),
+     rprocess=euler.sim(
+       step.fun=Csnippet(stochsim),
+       delta.t=1/24
+       ),
+     skeleton=Csnippet(skel),
+     skelmap.delta.t=1/24,
+     skeleton.type="map",
+     rmeasure=Csnippet(rmeas),
+     dmeasure=Csnippet(dmeas),
+     parameter.transform=Csnippet(partrans),
+     parameter.inv.transform=Csnippet(paruntrans),
+     obsnames = c("confined","convalescent"),
+     statenames=c("S","I","R","C"),
+     paramnames=c(
+       "Beta",
+       "inf.pd","conf.pd","conv.pd",
+       "rho","sfrac"
+       ),
+     initializer=function(params, t0, ...) {
+       x0 <- setNames(numeric(4),c("S","I","R","C"))
+       S.0 <- round(763*params["sfrac"])
+       x0[c("S","I")] <- c(S.0,763-S.0)
+       x0
+     }
+     ) -> bsflu
+
+c("bsflu")

Added: pkg/pompExamples/inst/examples/bsflu3.R
===================================================================
--- pkg/pompExamples/inst/examples/bsflu3.R	                        (rev 0)
+++ pkg/pompExamples/inst/examples/bsflu3.R	2015-01-23 14:50:29 UTC (rev 1053)
@@ -0,0 +1,117 @@
+## data read from graph in Anonymous (1978) 'Influenza in a boarding school' Brit. Med. J. 1:578.
+## cases are recorded with error of +/- 1 case
+## 763 boys were at risk, 512 boys spent time away from class
+flu <- read.csv(text="
+date,confined,convalescent
+1978-01-22,1,0
+1978-01-23,6,0
+1978-01-24,26,0
+1978-01-25,73,1
+1978-01-26,222,8
+1978-01-27,293,16
+1978-01-28,258,99
+1978-01-29,236,160
+1978-01-30,191,173
+1978-01-31,124,162
+1978-02-01,69,150
+1978-02-02,26,89
+1978-02-03,11,44
+1978-02-04,4,22
+",colClasses=c(date='Date'))
+flu$day <- flu$date-min(flu$date)+1
+units(flu$day) <- "days"
+flu$day <- as.numeric(flu$day)
+
+partrans <- "
+ TBeta = exp(Beta);
+ Tinf_pd = exp(inf_pd);
+ Trho = expit(rho);
+ Tsfrac = expit(sfrac);
+"
+
+paruntrans <- "
+ TBeta = log(Beta);
+ Tinf_pd = log(inf_pd);
+ Trho = logit(rho);
+ Tsfrac = logit(sfrac);
+"
+
+dmeas <- "
+  lik = dpois(confined,rho*(R1+R2+R3)+1e-6,give_log);
+"
+
+rmeas <- "
+  confined = rpois(rho*(R1+R2+R3)+1e-6);
+  convalescent = rpois(rho*C);
+"
+
+stochsim <- "
+  double t1 = rbinom(S,1-exp(-Beta*I*dt));
+  double t2 = rbinom(I,1-exp(-dt/inf_pd));
+  double t3a = rbinom(R1,1-exp(-3*dt/conf_pd));
+  double t3b = rbinom(R2,1-exp(-3*dt/conf_pd));
+  double t3c = rbinom(R3,1-exp(-3*dt/conf_pd));
+  double t4 = rbinom(C,1-exp(-dt/conv_pd));
+  S -= t1;
+  I += t1 - t2;
+  R1 += t2 - t3a;
+  R2 += t3a - t3b;
+  R3 += t3b - t3c;
+  C += t3c - t4;
+"
+
+skel <- "
+  double dt = 1.0/24.0;
+  double t1 = S*(1-exp(-Beta*I*dt));
+  double t2 = I*(1-exp(-dt/inf_pd));
+  double t3a = R1*(1-exp(-3*dt/conf_pd));
+  double t3b = R2*(1-exp(-3*dt/conf_pd));
+  double t3c = R3*(1-exp(-3*dt/conf_pd));
+  double t4 = C*(1-exp(-dt/conv_pd));
+  DS = S - t1;
+  DI = I + t1 - t2;
+  DR1 = R1 + t2 - t3a;
+  DR2 = R2 + t3a - t3b;
+  DR3 = R3 + t3b - t3c;
+  DC = C + t3c - t4;
+"
+
+pomp(
+     data=flu[c("day","confined","convalescent")],
+     times="day",
+     t0=0,
+     params=c(
+       Beta=0.004,
+       inf.pd=0.7,
+       conf.pd=sum(flu$confined)/512,
+       conv.pd=sum(flu$convalescent)/512,
+       rho=0.9,
+       sfrac=762/763
+       ),
+     rprocess=euler.sim(
+       step.fun=Csnippet(stochsim),
+       delta.t=1/24
+       ),
+     skeleton=Csnippet(skel),
+     skelmap.delta.t=1/24,
+     skeleton.type="map",
+     rmeasure=Csnippet(rmeas),
+     dmeasure=Csnippet(dmeas),
+     parameter.transform=Csnippet(partrans),
+     parameter.inv.transform=Csnippet(paruntrans),
+     obsnames = c("confined","convalescent"),
+     statenames=c("S","I","R1","R2","R3","C"),
+     paramnames=c(
+       "Beta",
+       "inf.pd","conf.pd","conv.pd",
+       "rho","sfrac"
+       ),
+     initializer=function(params, t0, ...) {
+       x0 <- setNames(numeric(6),c("S","I","R1","R2","R3","C"))
+       S.0 <- round(763*params["sfrac"])
+       x0[c("S","I")] <- c(S.0,763-S.0)
+       x0
+     }
+     ) -> bsflu3
+
+c("bsflu3")

Added: pkg/pompExamples/vignettes/Makefile
===================================================================
--- pkg/pompExamples/vignettes/Makefile	                        (rev 0)
+++ pkg/pompExamples/vignettes/Makefile	2015-01-23 14:50:29 UTC (rev 1053)
@@ -0,0 +1,10 @@
+RSCRIPT = Rscript --vanilla
+RM = rm -f
+
+vignettes: bsflu.html
+	$(RM) -r figure
+
+%.html: %.Rmd
+	$(RSCRIPT) -e 'knitr::knit("$*.Rmd")'
+	$(RSCRIPT) -e 'markdown::markdownToHTML("$*.md","$*.html")'
+	$(RM) $*.md

Added: pkg/pompExamples/vignettes/bsflu-mf1.rds
===================================================================
--- pkg/pompExamples/vignettes/bsflu-mf1.rds	                        (rev 0)
+++ pkg/pompExamples/vignettes/bsflu-mf1.rds	2015-01-23 14:50:29 UTC (rev 1053)
@@ -0,0 +1 @@
+‹      ‹àb```b`fdb`b2™…ó1だm[Y{`|ûo®nÙŽM†óóŽnJó™êý¥5Æwà`€ˆüŸ¨ü·!¬ÛßFyÂùï¿:¤þe``añX84k^bnj1! vD-'?='3Êã„ðôŠS¡,N©%‰0¥™yiz)P{r>:·Áe.ÊȇY[œV”˜Œæ΢ür=˜{xA× $þº£“s‹aŽ†	r¥$–$êÍMòþ ý¹„&~  
\ No newline at end of file

Added: pkg/pompExamples/vignettes/bsflu-mf3.rds
===================================================================
--- pkg/pompExamples/vignettes/bsflu-mf3.rds	                        (rev 0)
+++ pkg/pompExamples/vignettes/bsflu-mf3.rds	2015-01-23 14:50:29 UTC (rev 1053)
@@ -0,0 +1,4 @@
+‹      ‹àb```b`fdb`b2™…ó1ãïŠ|ôcÓ`|û]_„¸VÍ…óóŠ¿|>翵ú1ëßþ>߁C€
+ ò òŸOòÛòÎPÎÿÕ¹ õ(ˆÏÀÂÀ	¤YósS‹°ã ‚l9ùé9™ÙP'„§Wœ
+`qJ-I„)ÍÌKÓ+HòØ“óѹe.sQF>ÌÚâ´¢Äd4·på—ëÁÜÃ
+¼ ñНœ“Xs4L+%±$Qhln*÷ ¦‚ý®~  
\ No newline at end of file

Added: pkg/pompExamples/vignettes/bsflu-tm1.rds
===================================================================
--- pkg/pompExamples/vignettes/bsflu-tm1.rds	                        (rev 0)
+++ pkg/pompExamples/vignettes/bsflu-tm1.rds	2015-01-23 14:50:29 UTC (rev 1053)
@@ -0,0 +1 @@
+‹      ‹àb```b`fdb`b2™…ó1ãØ…UÑFÚ7`| óíóds'MuóŸËKì`mm„ñ8PÕÿy€Ê]è,üþÁM8ÿýWç‚Ôÿ¡,L`ûX84k^bnj1! vD-'?='3Êã„ðôŠS¡,N©%‰0¥™yiz)P{r>:·Áe.ÊȇY[œV”˜Œæ΢ür=˜{xA× $þº£“s‹aŽ†	r¥$–$êÍMòþ î³£¡~  
\ No newline at end of file

Added: pkg/pompExamples/vignettes/bsflu-tm3.rds
===================================================================
--- pkg/pompExamples/vignettes/bsflu-tm3.rds	                        (rev 0)
+++ pkg/pompExamples/vignettes/bsflu-tm3.rds	2015-01-23 14:50:29 UTC (rev 1053)
@@ -0,0 +1,4 @@
+‹      ‹àb```b`fdb`b2™…ó1ã?Ï—*S,Ua| óí³ß+­v7σó_éþx»TýŒïÀ!€ªþÏTþKÛFÿ©2pþû¯Î©ÿCX˜Àö±0piÖ¼ÄÜÔb C ì8ˆ [N~zNf6”Ç	áé§BXœRKaJ3óÒô
+R <öä|tn‚Ë\”‘³¶8­(1Í-œEùåz0÷ð‚¯HütG'ç$Ã
+äJI,IÔ››
+äý Q&~  
\ No newline at end of file

Added: pkg/pompExamples/vignettes/bsflu.Rmd
===================================================================
--- pkg/pompExamples/vignettes/bsflu.Rmd	                        (rev 0)
+++ pkg/pompExamples/vignettes/bsflu.Rmd	2015-01-23 14:50:29 UTC (rev 1053)
@@ -0,0 +1,183 @@
+%\VignetteIndexEntry{Boarding-school influenza example}
+%\VignetteEngine{knitr::knitr}
+
+# Boarding-School Flu Outbreak Analysis  
+***Aaron A. King***
+
+```{r setup,include=FALSE}
+require(pomp)
+require(plyr)
+require(reshape2)
+options(stringsAsFactors=FALSE,keep.source=TRUE,encoding="UTF-8")
+
+require(ggplot2)
+theme_set(theme_bw())
+
+require(knitr)
+opts_knit$set(out.format="html")
+opts_chunk$set(
+  progress=TRUE,
+  prompt=FALSE,tidy=FALSE,highlight=TRUE,
+  strip.white=TRUE,
+  warning=FALSE,message=FALSE,error=FALSE,
+  echo=TRUE,cache=FALSE,
+  results='markup',
+  fig.show='asis',
+  fig.height=5,fig.width=10,
+  dpi=100
+  )
+
+require(pompExamples)
+set.seed(862663052L)
+```
+
+First, a little function to cache the results of expensive computations.
+
+```{r juliaChild}
+juliaChild <- function (file, expr) {
+  if (file.exists(file)) {
+    readRDS(file)
+    } else {
+      val <- eval(expr)
+      saveRDS(val,file=file)
+      val
+      }
+  }
+```
+
+## Flu model with exponentially-distributed waiting times
+
+```{r bsflu1,cache=FALSE}
+pompExample(bsflu)
+coef(bsflu)
+```
+```{r plot1}
+simdat <- simulate(bsflu,nsim=10,obs=TRUE,as.data.frame=TRUE,include.data=TRUE)
+ggplot(data=melt(simdat,id=c("sim","time")),
+       mapping=aes(x=time,y=value,group=interaction(variable,sim),
+                   color=variable,size=sim=="data",alpha=sim=="data"))+
+  geom_line()+
+  scale_alpha_manual(values=c(`TRUE`=1,`FALSE`=0.5))+
+  scale_size_manual(values=c(`TRUE`=2,`FALSE`=1))
+```
+
+### Trajectory matching
+
+```{r tm1}
+juliaChild("bsflu-tm1.rds",{
+  traj.match(bsflu,est=c("Beta","inf.pd","rho"),transform=TRUE) -> tm
+  traj.match(tm,method='subplex') -> tm
+  data.frame(loglik=logLik(tm),loglik.se=0,as.list(coef(tm)))
+}) -> tm1
+```
+```{r tm1-simplot}
+simdat <- simulate(bsflu,params=unlist(tm1),nsim=10,
+                   obs=TRUE,as.data.frame=TRUE,include.data=TRUE)
+ggplot(data=melt(simdat,id=c("sim","time")),
+       mapping=aes(x=time,y=value,group=interaction(variable,sim),
+                   color=variable,size=sim=="data",alpha=sim=="data"))+
+  geom_line()+
+  scale_alpha_manual(values=c(`TRUE`=1,`FALSE`=0.5))+
+  scale_size_manual(values=c(`TRUE`=2,`FALSE`=1))
+```
+
+### Iterated filtering
+
+```{r mf1}
+juliaChild("bsflu-mf1.rds",{
+  mif(bsflu,rw.sd=c(Beta=0.05,inf.pd=0.05,rho=0.05),
+      cooling.fraction=0.9,var.factor=2,
+      Nmif=50,Np=1000,method='mif2',transform=TRUE) -> mf
+  mif(mf,Nmif=50,cooling.fraction=0.5) -> mf
+  mif(mf,Nmif=50,cooling.fraction=0.1) -> mf
+  ll <- unname(logmeanexp(raply(5,logLik(pfilter(mf))),se=TRUE))
+  data.frame(loglik=ll[1],loglik.se=ll[2],as.list(coef(mf)))
+}) -> mf1
+```
+```{r mf1-simplot}
+simdat <- simulate(bsflu,params=unlist(mf1),nsim=10,
+                   obs=TRUE,as.data.frame=TRUE,include.data=TRUE)
+ggplot(data=melt(simdat,id=c("sim","time")),
+       mapping=aes(x=time,y=value,group=interaction(variable,sim),
+                   color=variable,size=sim=="data",alpha=sim=="data"))+
+  geom_line()+
+  scale_alpha_manual(values=c(`TRUE`=1,`FALSE`=0.5))+
+  scale_size_manual(values=c(`TRUE`=2,`FALSE`=1))
+```
+
+## Model with Erlang(3) confinement period
+
+```{r bsflu3,cache=FALSE}
+pompExample(bsflu3)
+coef(bsflu3)
+```
+```{r plot3}
+simdat <- simulate(bsflu3,nsim=10,obs=TRUE,as.data.frame=TRUE,include.data=TRUE)
+ggplot(data=melt(simdat,id=c("sim","time")),
+       mapping=aes(x=time,y=value,group=interaction(variable,sim),
+                   color=variable,size=sim=="data",alpha=sim=="data"))+
+  geom_line()+
+  scale_alpha_manual(values=c(`TRUE`=1,`FALSE`=0.5))+
+  scale_size_manual(values=c(`TRUE`=2,`FALSE`=1))
+```
+
+### Trajectory matching
+
+```{r tm3}
+juliaChild("bsflu-tm3.rds",{
+  traj.match(bsflu3,est=c("Beta","inf.pd","rho"),transform=TRUE) -> tm
+  traj.match(tm,method='subplex') -> tm
+  data.frame(loglik=logLik(tm),loglik.se=0,as.list(coef(tm)))
+}) -> tm3
+```
+```{r tm3-simplot}
+simdat <- simulate(bsflu3,params=unlist(tm3),nsim=10,
+                   obs=TRUE,as.data.frame=TRUE,include.data=TRUE)
+ggplot(data=melt(simdat,id=c("sim","time")),
+       mapping=aes(x=time,y=value,group=interaction(variable,sim),
+                   color=variable,size=sim=="data",alpha=sim=="data"))+
+  geom_line()+
+  scale_alpha_manual(values=c(`TRUE`=1,`FALSE`=0.5))+
+  scale_size_manual(values=c(`TRUE`=2,`FALSE`=1))
+```
+
+### Iterated filtering
+
+```{r mf3}
+juliaChild("bsflu-mf3.rds",{
+  mif(bsflu3,rw.sd=c(Beta=0.05,inf.pd=0.05,rho=0.05),
+      cooling.fraction=0.9,var.factor=2,
+      Nmif=50,Np=1000,method='mif2',transform=TRUE) -> mf
+  mif(mf,Nmif=50,cooling.fraction=0.5) -> mf
+  mif(mf,Nmif=50,cooling.fraction=0.1) -> mf
+  ll <- unname(logmeanexp(raply(5,logLik(pfilter(mf))),se=TRUE))
+  data.frame(loglik=ll[1],loglik.se=ll[2],as.list(coef(mf)))
+}) -> mf3
+```
+```{r mf3-simplot}
+simdat <- simulate(bsflu3,params=unlist(mf3),nsim=10,
+                   obs=TRUE,as.data.frame=TRUE,include.data=TRUE)
+ggplot(data=melt(simdat,id=c("sim","time")),
+       mapping=aes(x=time,y=value,group=interaction(variable,sim),
+                   color=variable,size=sim=="data",alpha=sim=="data"))+
+  geom_line()+
+  scale_alpha_manual(values=c(`TRUE`=1,`FALSE`=0.5))+
+  scale_size_manual(values=c(`TRUE`=2,`FALSE`=1))
+```
+
+## Model comparison
+
+```{r comp}
+tb <- ldply(list(det1=tm1,stoch1=mf1,det3=tm3,stoch3=mf3),.id="model")
+kable(tb)
+```
+
+```{r probes,fig.height=6,fig.width=6}
+plist <- list(probe.acf("confined",lags=c(2),type="cov",transform=sqrt),
+              probe.ccf(c("confined","convalescent"),lags=c(0,3),transform=sqrt),
+              tot=function (x) apply(x,1,sum))
+plot(probe(bsflu,params=unlist(tm1),nsim=500,probes=plist))
+plot(probe(bsflu,params=unlist(mf1),nsim=500,probes=plist))
+plot(probe(bsflu3,params=unlist(tm3),nsim=500,probes=plist))
+plot(probe(bsflu3,params=unlist(mf3),nsim=500,probes=plist))
+```

Added: pkg/pompExamples/vignettes/bsflu.html
===================================================================
--- pkg/pompExamples/vignettes/bsflu.html	                        (rev 0)
+++ pkg/pompExamples/vignettes/bsflu.html	2015-01-23 14:50:29 UTC (rev 1053)
@@ -0,0 +1,460 @@
+<!DOCTYPE html>
+<html>
+<head>
+<meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
+
+<title>Boarding-School Flu Outbreak Analysis</title>
+
+<script type="text/javascript">
+window.onload = function() {
+  var imgs = document.getElementsByTagName('img'), i, img;
+  for (i = 0; i < imgs.length; i++) {
+    img = imgs[i];
+    // center an image if it is the only element of its parent
+    if (img.parentElement.childElementCount === 1)
+      img.parentElement.style.textAlign = 'center';
+  }
+};
+</script>
+
+<!-- Styles for R syntax highlighter -->
+<style type="text/css">
+   pre .operator,
+   pre .paren {
+     color: rgb(104, 118, 135)
+   }
+
+   pre .literal {
+     color: #990073
+   }
+
+   pre .number {
+     color: #099;
+   }
+
+   pre .comment {
+     color: #998;
+     font-style: italic
+   }
+
+   pre .keyword {
+     color: #900;
+     font-weight: bold
+   }
+
+   pre .identifier {
+     color: rgb(0, 0, 0);
+   }
+
+   pre .string {
+     color: #d14;
+   }
+</style>
+
+<!-- R syntax highlighter -->
+<script type="text/javascript">
+var hljs=new function(){function m(p){return p.replace(/&/gm,"&").replace(/</gm,"<")}function f(r,q,p){return RegExp(q,"m"+(r.cI?"i":"")+(p?"g":""))}function b(r){for(var p=0;p<r.childNodes.length;p++){var q=r.childNodes[p];if(q.nodeName=="CODE"){return q}if(!(q.nodeType==3&&q.nodeValue.match(/\s+/))){break}}}function h(t,s){var p="";for(var r=0;r<t.childNodes.length;r++){if(t.childNodes[r].nodeType==3){var q=t.childNodes[r].nodeValue;if(s){q=q.replace(/\n/g,"")}p+=q}else{if(t.childNodes[r].nodeName=="BR"){p+="\n"}else{p+=h(t.childNodes[r])}}}if(/MSIE [678]/.test(navigator.userAgent)){p=p.replace(/\r/g,"\n")}return p}function a(s){var r=s.className.split(/\s+/);r=r.concat(s.parentNode.className.split(/\s+/));for(var q=0;q<r.length;q++){var p=r[q].replace(/^language-/,"");if(e[p]){return p}}}function c(q){var p=[];(function(s,t){for(var r=0;r<s.childNodes.length;r++){if(s.childNodes[r].nodeType==3){t+=s.childNodes[r].nodeValue.length}else{if(s.childNodes[r].nodeName=="BR"){t+=1}else{if(s.childNodes[r].nodeType==1){p.push({event:"start",offset:t,node:s.childNodes[r]});t=arguments.callee(s.childNodes[r],t);p.push({event:"stop",offset:t,node:s.childNodes[r]})}}}}return t})(q,0);return p}function k(y,w,x){var q=0;var z="";var s=[];function u(){if(y.length&&w.length){if(y[0].offset!=w[0].offset){return(y[0].offset<w[0].offset)?y:w}else{return w[0].event=="start"?y:w}}else{return y.length?y:w}}function t(D){var A="<"+D.nodeName.toLowerCase();for(var B=0;B<D.attributes.length;B++){var C=D.attributes[B];A+=" "+C.nodeName.toLowerCase();if(C.value!==undefined&&C.value!==false&&C.value!==null){A+='="'+m(C.value)+'"'}}return A+">"}while(y.length||w.length){var v=u().splice(0,1)[0];z+=m(x.substr(q,v.offset-q));q=v.offset;if(v.event=="start"){z+=t(v.node);s.push(v.node)}else{if(v.event=="stop"){var p,r=s.length;do{r--;p=s[r];z+=("</"+p.nodeName.toLowerCase()+">")}while(p!=v.node);s.splice(r,1);while(r<s.length){z+=t(s[r]);r++}}}}return z+m(x.substr(q))}function j(){function q(x,y,v){if(x.compiled){return}var u;var s=[];if(x.k){x.lR=f(y,x.l||hljs.IR,true);for(var w in x.k){if(!x.k.hasOwnProperty(w)){continue}if(x.k[w] instanceof Object){u=x.k[w]}else{u=x.k;w="keyword"}for(var r in u){if(!u.hasOwnProperty(r)){continue}x.k[r]=[w,u[r]];s.push(r)}}}if(!v){if(x.bWK){x.b="\\b("+s.join("|")+")\\s"}x.bR=f(y,x.b?x.b:"\\B|\\b");if(!x.e&&!x.eW){x.e="\\B|\\b"}if(x.e){x.eR=f(y,x.e)}}if(x.i){x.iR=f(y,x.i)}if(x.r===undefined){x.r=1}if(!x.c){x.c=[]}x.compiled=true;for(var t=0;t<x.c.length;t++){if(x.c[t]=="self"){x.c[t]=x}q(x.c[t],y,false)}if(x.starts){q(x.starts,y,false)}}for(var p in e){if(!e.hasOwnProperty(p)){continue}q(e[p].dM,e[p],true)}}function d(B,C){if(!j.called){j();j.called=true}function q(r,M){for(var L=0;L<M.c.length;L++){if((M.c[L].bR.exec(r)||[null])[0]==r){return M.c[L]}}}function v(L,r){if(D[L].e&&D[L].eR.test(r)){return 1}if(D[L].eW){var M=v(L-1,r);return M?M+1:0}return 0}function w(r,L){return L.i&&L.iR.test(r)}function K(N,O){var M=[];for(var L=0;L<N.c.length;L++){M.push(N.c[L].b)}var r=D.length-1;do{if(D[r].e){M.push(D[r].e)}r--}while(D[r+1].eW);if(N.i){M.push(N.i)}return f(O,M.join("|"),true)}function p(M,L){var N=D[D.length-1];if(!N.t){N.t=K(N,E)}N.t.lastIndex=L;var r=N.t.exec(M);return r?[M.substr(L,r.index-L),r[0],false]:[M.substr(L),"",true]}function z(N,r){var L=E.cI?r[0].toLowerCase():r[0];var M=N.k[L];if(M&&M instanceof Array){return M}return false}function F(L,P){L=m(L);if(!P.k){return L}var r="";var O=0;P.lR.lastIndex=0;var M=P.lR.exec(L);while(M){r+=L.substr(O,M.index-O);var N=z(P,M);if(N){x+=N[1];r+='<span class="'+N[0]+'">'+M[0]+"</span>"}else{r+=M[0]}O=P.lR.lastIndex;M=P.lR.exec(L)}return r+L.substr(O,L.length-O)}function J(L,M){if(M.sL&&e[M.sL]){var r=d(M.sL,L);x+=r.keyword_count;return r.value}else{return F(L,M)}}function I(M,r){var L=M.cN?'<span class="'+M.cN+'">':"";if(M.rB){y+=L;M.buffer=""}else{if(M.eB){y+=m(r)+L;M.buffer=""}else{y+=L;M.buffer=r}}D.push(M);A+=M.r}function G(N,M,Q){var R=D[D.length-1];if(Q){y+=J(R.buffer+N,R);return false}var P=q(M,R);if(P){y+=J(R.buffer+N,R);I(P,M);return P.rB}var L=v(D.length-1,M);if(L){var O=R.cN?"</span>":"";if(R.rE){y+=J(R.buffer+N,R)+O}else{if(R.eE){y+=J(R.buffer+N,R)+O+m(M)}else{y+=J(R.buffer+N+M,R)+O}}while(L>1){O=D[D.length-2].cN?"</span>":"";y+=O;L--;D.length--}var r=D[D.length-1];D.length--;D[D.length-1].buffer="";if(r.starts){I(r.starts,"")}return R.rE}if(w(M,R)){throw"Illegal"}}var E=e[B];var D=[E.dM];var A=0;var x=0;var y="";try{var s,u=0;E.dM.buffer="";do{s=p(C,u);var t=G(s[0],s[1],s[2]);u+=s[0].length;if(!t){u+=s[1].length}}while(!s[2]);if(D.length>1){throw"Illegal"}return{r:A,keyword_count:x,value:y}}catch(H){if(H=="Illegal"){return{r:0,keyword_count:0,value:m(C)}}else{throw H}}}function g(t){var p={keyword_count:0,r:0,value:m(t)};var r=p;for(var q in e){if(!e.hasOwnProperty(q)){continue}var s=d(q,t);s.language=q;if(s.keyword_count+s.r>r.keyword_count+r.r){r=s}if(s.keyword_count+s.r>p.keyword_count+p.r){r=p;p=s}}if(r.language){p.second_best=r}return p}function i(r,q,p){if(q){r=r.replace(/^((<[^>]+>|\t)+)/gm,function(t,w,v,u){return w.replace(/\t/g,q)})}if(p){r=r.replace(/\n/g,"<br>")}return r}function n(t,w,r){var x=h(t,r);var v=a(t);var y,s;if(v){y=d(v,x)}else{return}var q=c(t);if(q.length){s=document.createElement("pre");s.innerHTML=y.value;y.value=k(q,c(s),x)}y.value=i(y.value,w,r);var u=t.className;if(!u.match("(\\s|^)(language-)?"+v+"(\\s|$)")){u=u?(u+" "+v):v}if(/MSIE [678]/.test(navigator.userAgent)&&t.tagName=="CODE"&&t.parentNode.tagName=="PRE"){s=t.parentNode;var p=document.createElement("div");p.innerHTML="<pre><code>"+y.value+"</code></pre>";t=p.firstChild.firstChild;p.firstChild.cN=s.cN;s.parentNode.replaceChild(p.firstChild,s)}else{t.innerHTML=y.value}t.className=u;t.result={language:v,kw:y.keyword_count,re:y.r};if(y.second_best){t.second_best={language:y.second_best.language,kw:y.second_best.keyword_count,re:y.second_best.r}}}function o(){if(o.called){return}o.called=true;var r=document.getElementsByTagName("pre");for(var p=0;p<r.length;p++){var q=b(r[p]);if(q){n(q,hljs.tabReplace)}}}function l(){if(window.addEventListener){window.addEventListener("DOMContentLoaded",o,false);window.addEventListener("load",o,false)}else{if(window.attachEvent){window.attachEvent("onload",o)}else{window.onload=o}}}var e={};this.LANGUAGES=e;this.highlight=d;this.highlightAuto=g;this.fixMarkup=i;this.highlightBlock=n;this.initHighlighting=o;this.initHighlightingOnLoad=l;this.IR="[a-zA-Z][a-zA-Z0-9_]*";this.UIR="[a-zA-Z_][a-zA-Z0-9_]*";this.NR="\\b\\d+(\\.\\d+)?";this.CNR="\\b(0[xX][a-fA-F0-9]+|(\\d+(\\.\\d*)?|\\.\\d+)([eE][-+]?\\d+)?)";this.BNR="\\b(0b[01]+)";this.RSR="!|!=|!==|%|%=|&|&&|&=|\\*|\\*=|\\+|\\+=|,|\\.|-|-=|/|/=|:|;|<|<<|<<=|<=|=|==|===|>|>=|>>|>>=|>>>|>>>=|\\?|\\[|\\{|\\(|\\^|\\^=|\\||\\|=|\\|\\||~";this.ER="(?![\\s\\S])";this.BE={b:"\\\\.",r:0};this.ASM={cN:"string",b:"'",e:"'",i:"\\n",c:[this.BE],r:0};this.QSM={cN:"string",b:'"',e:'"',i:"\\n",c:[this.BE],r:0};this.CLCM={cN:"comment",b:"//",e:"$"};this.CBLCLM={cN:"comment",b:"/\\*",e:"\\*/"};this.HCM={cN:"comment",b:"#",e:"$"};this.NM={cN:"number",b:this.NR,r:0};this.CNM={cN:"number",b:this.CNR,r:0};this.BNM={cN:"number",b:this.BNR,r:0};this.inherit=function(r,s){var p={};for(var q in r){p[q]=r[q]}if(s){for(var q in s){p[q]=s[q]}}return p}}();hljs.LANGUAGES.cpp=function(){var a={keyword:{"false":1,"int":1,"float":1,"while":1,"private":1,"char":1,"catch":1,"export":1,virtual:1,operator:2,sizeof:2,dynamic_cast:2,typedef:2,const_cast:2,"const":1,struct:1,"for":1,static_cast:2,union:1,namespace:1,unsigned:1,"long":1,"throw":1,"volatile":2,"static":1,"protected":1,bool:1,template:1,mutable:1,"if":1,"public":1,friend:2,"do":1,"return":1,"goto":1,auto:1,"void":2,"enum":1,"else":1,"break":1,"new":1,extern:1,using:1,"true":1,"class":1,asm:1,"case":1,typeid:1,"short":1,reinterpret_cast:2,"default":1,"double":1,register:1,explicit:1,signed:1,typename:1,"try":1,"this":1,"switch":1,"continue":1,wchar_t:1,inline:1,"delete":1,alignof:1,char16_t:1,char32_t:1,constexpr:1,decltype:1,noexcept:1,nullptr:1,static_assert:1,thread_local:1,restrict:1,_Bool:1,complex:1},built_in:{std:1,string:1,cin:1,cout:1,cerr:1,clog:1,stringstream:1,istringstream:1,ostringstream:1,auto_ptr:1,deque:1,list:1,queue:1,stack:1,vector:1,map:1,set:1,bitset:1,multiset:1,multimap:1,unordered_set:1,unordered_map:1,unordered_multiset:1,unordered_multimap:1,array:1,shared_ptr:1}};return{dM:{k:a,i:"</",c:[hljs.CLCM,hljs.CBLCLM,hljs.QSM,{cN:"string",b:"'\\\\?.",e:"'",i:"."},{cN:"number",b:"\\b(\\d+(\\.\\d*)?|\\.\\d+)(u|U|l|L|ul|UL|f|F)"},hljs.CNM,{cN:"preprocessor",b:"#",e:"$"},{cN:"stl_container",b:"\\b(deque|list|queue|stack|vector|map|set|bitset|multiset|multimap|unordered_map|unordered_set|unordered_multiset|unordered_multimap|array)\\s*<",e:">",k:a,r:10,c:["self"]}]}}}();hljs.LANGUAGES.r={dM:{c:[hljs.HCM,{cN:"number",b:"\\b0[xX][0-9a-fA-F]+[Li]?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+(?:[eE][+\\-]?\\d*)?L\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+\\.(?!\\d)(?:i\\b)?",e:hljs.IMMEDIATE_RE,r:1},{cN:"number",b:"\\b\\d+(?:\\.\\d*)?(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\.\\d+(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"keyword",b:"(?:tryCatch|library|setGeneric|setGroupGeneric)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\.",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\d+(?![\\w.])",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\b(?:function)",e:hljs.IMMEDIATE_RE,r:2},{cN:"keyword",b:"(?:if|in|break|next|repeat|else|for|return|switch|while|try|stop|warning|require|attach|detach|source|setMethod|setClass)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"literal",b:"(?:NA|NA_integer_|NA_real_|NA_character_|NA_complex_)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"literal",b:"(?:NULL|TRUE|FALSE|T|F|Inf|NaN)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"identifier",b:"[a-zA-Z.][a-zA-Z0-9._]*\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"<\\-(?!\\s*\\d)",e:hljs.IMMEDIATE_RE,r:2},{cN:"operator",b:"\\->|<\\-",e:hljs.IMMEDIATE_RE,r:1},{cN:"operator",b:"%%|~",e:hljs.IMMEDIATE_RE},{cN:"operator",b:">=|<=|==|!=|\\|\\||&&|=|\\+|\\-|\\*|/|\\^|>|<|!|&|\\||\\$|:",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"%",e:"%",i:"\\n",r:1},{cN:"identifier",b:"`",e:"`",r:0},{cN:"string",b:'"',e:'"',c:[hljs.BE],r:0},{cN:"string",b:"'",e:"'",c:[hljs.BE],r:0},{cN:"paren",b:"[[({\\])}]",e:hljs.IMMEDIATE_RE,r:0}]}};
+hljs.initHighlightingOnLoad();
+</script>
+
+
+
+<style type="text/css">
+body, td {
+   font-family: sans-serif;
+   background-color: white;
+   font-size: 13px;
+}
+
+body {
+  max-width: 800px;
+  margin: auto;
+  padding: 1em;
+  line-height: 20px;
+}
+
+tt, code, pre {
+   font-family: 'DejaVu Sans Mono', 'Droid Sans Mono', 'Lucida Console', Consolas, Monaco, monospace;
+}
+
+h1 {
+   font-size:2.2em;
+}
+
+h2 {
+   font-size:1.8em;
+}
+
+h3 {
+   font-size:1.4em;
+}
+
+h4 {
+   font-size:1.0em;
+}
+
+h5 {
+   font-size:0.9em;
+}
+
+h6 {
+   font-size:0.8em;
+}
+
+a:visited {
+   color: rgb(50%, 0%, 50%);
+}
+
+pre, img {
+  max-width: 100%;
+}
+pre {
+  overflow-x: auto;
+}
+pre code {
+   display: block; padding: 0.5em;
+}
+
+code {
+  font-size: 92%;
+  border: 1px solid #ccc;
+}
+
+code[class] {
+  background-color: #F8F8F8;
+}
+
+table, td, th {
+  border: none;
+}
+
+blockquote {
+   color:#666666;
+   margin:0;
+   padding-left: 1em;
+   border-left: 0.5em #EEE solid;
+}
+
+hr {
+   height: 0px;
+   border-bottom: none;
+   border-top-width: thin;
+   border-top-style: dotted;
+   border-top-color: #999999;
+}
+
+ at media print {
+   * {
+      background: transparent !important;
+      color: black !important;
+      filter:none !important;
+      -ms-filter: none !important;
+   }
+
+   body {
+      font-size:12pt;
+      max-width:100%;
+   }
+
+   a, a:visited {
+      text-decoration: underline;
+   }
+
+   hr {
+      visibility: hidden;
+      page-break-before: always;
+   }
+
+   pre, blockquote {
+      padding-right: 1em;
+      page-break-inside: avoid;
+   }
+
+   tr, img {
+      page-break-inside: avoid;
+   }
+
+   img {
+      max-width: 100% !important;
+   }
+
+   @page :left {
+      margin: 15mm 20mm 15mm 10mm;
+   }
+
+   @page :right {
+      margin: 15mm 10mm 15mm 20mm;
+   }
+
+   p, h2, h3 {
+      orphans: 3; widows: 3;
+   }
+
+   h2, h3 {
+      page-break-after: avoid;
+   }
+}
+</style>
+
+
+
+</head>
+
+<body>
+<h1>Boarding-School Flu Outbreak Analysis</h1>
+
+<p><strong><em>Aaron A. King</em></strong></p>
+
+<p>First, a little function to cache the results of expensive computations.</p>
+
+<pre><code class="r">juliaChild <- function (file, expr) {
+  if (file.exists(file)) {
+    readRDS(file)
+    } else {
+      val <- eval(expr)
+      saveRDS(val,file=file)
+      val
+      }
+  }
+</code></pre>
+
+<h2>Flu model with exponentially-distributed waiting times</h2>
+
+<pre><code class="r">pompExample(bsflu)
+</code></pre>
+
+<pre><code>## newly created object(s):
+##  bsflu
+</code></pre>
+
+<pre><code class="r">coef(bsflu)
+</code></pre>
+
+<pre><code>##      Beta    inf.pd   conf.pd   conv.pd       rho     sfrac 
+## 0.0040000 0.7000000 3.0078125 1.8046875 0.9000000 0.9986894
+</code></pre>
+
+<pre><code class="r">simdat <- simulate(bsflu,nsim=10,obs=TRUE,as.data.frame=TRUE,include.data=TRUE)
+ggplot(data=melt(simdat,id=c("sim","time")),
+       mapping=aes(x=time,y=value,group=interaction(variable,sim),
+                   color=variable,size=sim=="data",alpha=sim=="data"))+
+  geom_line()+
+  scale_alpha_manual(values=c(`TRUE`=1,`FALSE`=0.5))+
+  scale_size_manual(values=c(`TRUE`=2,`FALSE`=1))
+</code></pre>
+
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/pomp -r 1053

More information about the pomp-commits mailing list