[Pomp-commits] r1145 - in www: content vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 23 13:59:09 CET 2015
Author: kingaa
Date: 2015-03-23 13:59:09 +0100 (Mon, 23 Mar 2015)
New Revision: 1145
Added:
www/vignettes/mif2.R
www/vignettes/mif2.Rmd
www/vignettes/mif2.html
Modified:
www/content/vignettes.htm
www/vignettes/Makefile
www/vignettes/pomp.pdf
Log:
- add mif2 vignette
Modified: www/content/vignettes.htm
===================================================================
--- www/content/vignettes.htm 2015-03-16 12:21:59 UTC (rev 1144)
+++ www/content/vignettes.htm 2015-03-23 12:59:09 UTC (rev 1145)
@@ -11,6 +11,11 @@
<td><a target="_blank" href="vignettes/getting_started.html">(HTML)</a></td>
<td><a target="_blank" href="vignettes/getting_started.R">(R code)</a></td>
</tr>
+<tr>
+<td>An example using IF2</td>
+<td><a target="_blank" href="vignettes/mif2.html">(HTML)</a></td>
+<td><a target="_blank" href="vignettes/mif2.R">(R code)</a></td>
+</tr>
<td><strong>pomp</strong> package manual</td>
<td><a target="_blank" href="vignettes/pomp.pdf">(PDF)</a></td>
</tr>
Modified: www/vignettes/Makefile
===================================================================
--- www/vignettes/Makefile 2015-03-16 12:21:59 UTC (rev 1144)
+++ www/vignettes/Makefile 2015-03-23 12:59:09 UTC (rev 1145)
@@ -1,5 +1,5 @@
RSCRIPT = $(R_HOME)/bin/Rscript --vanilla
-REXE = $(R_HOME)/bin/R --vanilla
+REXE = $(R_HOME)/bin/R --vanilla --slave
LATEX = latex
BIBTEX = bibtex
PDFLATEX = pdflatex
Added: www/vignettes/mif2.R
===================================================================
--- www/vignettes/mif2.R (rev 0)
+++ www/vignettes/mif2.R 2015-03-23 12:59:09 UTC (rev 1145)
@@ -0,0 +1,115 @@
+## ----options,include=FALSE,cache=FALSE-----------------------------------
+library(knitr)
+prefix <- "mif2"
+opts_chunk$set(
+ progress=TRUE,
+ prompt=FALSE,tidy=FALSE,highlight=TRUE,
+ strip.white=TRUE,
+ warning=FALSE,
+ message=FALSE,
+ error=FALSE,
+ echo=TRUE,
+ cache=TRUE,
+ results='markup',
+ fig.show='asis',
+ size='small',
+ fig.lp="fig:",
+ fig.path=paste0("figure/",prefix,"-"),
+ cache.path=paste0("cache/",prefix,"-"),
+ fig.pos="h!",
+ fig.align='center',
+ fig.height=4,fig.width=6.83,
+ dpi=300,
+ dev='png',
+ dev.args=list(bg='transparent')
+ )
+require(ggplot2)
+require(plyr)
+require(reshape2)
+require(magrittr)
+require(pomp)
+theme_set(theme_bw())
+stopifnot(packageVersion("pomp")>="0.63-1")
+options(
+ keep.source=TRUE,
+ stringsAsFactors=FALSE,
+ encoding="UTF-8",
+ scipen=5,
+ cores=5
+ )
+
+## ----gompertz-init,cache=FALSE-------------------------------------------
+require(pomp)
+pompExample(gompertz)
+theta <- coef(gompertz)
+theta.true <- theta
+
+## ----gompertz-multi-mif2-eval,results='hide'-----------------------------
+require(foreach)
+require(doMC)
+registerDoMC()
+
+save.seed <- .Random.seed
+set.seed(334388458L,kind="L'Ecuyer")
+
+estpars <- c("r","sigma","tau")
+mf <- foreach(i=1:10,
+ .inorder=FALSE,
+ .options.multicore=list(set.seed=TRUE)
+ ) %dopar%
+ {
+ theta.guess <- theta.true
+ theta.guess[estpars] <- rlnorm(
+ n=length(estpars),
+ meanlog=log(theta.guess[estpars]),
+ sdlog=1
+ )
+ m1 <- mif(
+ gompertz,
+ Nmif=50,
+ method="mif2",
+ start=theta.guess,
+ transform=TRUE,
+ rw.sd=c(r=0.02,sigma=0.02,tau=0.05),
+ Np=2000,
+ var.factor=2,
+ cooling.type="hyperbolic",
+ cooling.fraction=0.95
+ )
+ m1 <- continue(m1,Nmif=50,cooling.fraction=0.8)
+ m1 <- continue(m1,Nmif=50,cooling.fraction=0.6)
+ m1 <- continue(m1,Nmif=50,cooling.fraction=0.2)
+ ll <- replicate(n=10,logLik(pfilter(m1,Np=10000)))
+ list(mif=m1,ll=ll)
+ }
+
+rbind(
+ mle=c(signif(theta.mif[estpars],3),loglik=round(loglik.mif,2)),
+ truth=c(signif(theta.true[estpars],3),loglik=round(loglik.true,2))
+ ) -> results.table
+
+## ----gompertz-post-mif2--------------------------------------------------
+theta.true <- coef(gompertz)
+loglik.true <- replicate(n=10,logLik(pfilter(gompertz,Np=10000)))
+loglik.true <- logmeanexp(loglik.true,se=TRUE)
+theta.mif <- t(sapply(mf,function(x)coef(x$mif)))
+loglik.mif <- t(sapply(mf,function(x)logmeanexp(x$ll,se=TRUE)))
+best <- which.max(loglik.mif[,1])
+theta.mif <- theta.mif[best,]
+loglik.mif <- loglik.mif[best,]
+
+## ----mif2-plot,echo=FALSE,cache=FALSE------------------------------------
+op <- par(mfrow=c(4,1),mar=c(3,3,0,0),mgp=c(2,1,0),bty='l')
+loglik <- sapply(mf,function(x)conv.rec(x$mif,"loglik"))
+log.r <- sapply(mf,function(x)conv.rec(x$mif,"r"))
+log.sigma <- sapply(mf,function(x)conv.rec(x$mif,"sigma"))
+log.tau <- sapply(mf,function(x)conv.rec(x$mif,"tau"))
+matplot(loglik,type='l',lty=1,xlab="",ylab=expression(log~L),xaxt='n',ylim=max(loglik,na.rm=T)+c(-12,3))
+matplot(log.r,type='l',lty=1,xlab="",ylab=expression(log~r),xaxt='n')
+matplot(log.sigma,type='l',lty=1,xlab="",ylab=expression(log~sigma),xaxt='n')
+matplot(log.tau,type='l',lty=1,xlab="MIF iteration",ylab=expression(log~tau))
+par(op)
+
+## ----first-mif-results-table,echo=FALSE,cache=FALSE----------------------
+print(results.table)
+
Added: www/vignettes/mif2.Rmd
===================================================================
--- www/vignettes/mif2.Rmd (rev 0)
+++ www/vignettes/mif2.Rmd 2015-03-23 12:59:09 UTC (rev 1145)
@@ -0,0 +1,152 @@
+---
+title: "An IF2 example"
+author: "Aaron A. King"
+output:
+ html_document:
+ theme: flatly
+ toc: yes
+bibliography: pomp.bib
+csl: ecology.csl
+---
+
+```{r options,include=FALSE,cache=FALSE}
+library(knitr)
+prefix <- "mif2"
+opts_chunk$set(
+ progress=TRUE,
+ prompt=FALSE,tidy=FALSE,highlight=TRUE,
+ strip.white=TRUE,
+ warning=FALSE,
+ message=FALSE,
+ error=FALSE,
+ echo=TRUE,
+ cache=TRUE,
+ results='markup',
+ fig.show='asis',
+ size='small',
+ fig.lp="fig:",
+ fig.path=paste0("figure/",prefix,"-"),
+ cache.path=paste0("cache/",prefix,"-"),
+ fig.pos="h!",
+ fig.align='center',
+ fig.height=4,fig.width=6.83,
+ dpi=300,
+ dev='png',
+ dev.args=list(bg='transparent')
+ )
+require(ggplot2)
+require(plyr)
+require(reshape2)
+require(magrittr)
+require(pomp)
+theme_set(theme_bw())
+stopifnot(packageVersion("pomp")>="0.63-1")
+options(
+ keep.source=TRUE,
+ stringsAsFactors=FALSE,
+ encoding="UTF-8",
+ scipen=5,
+ cores=5
+ )
+```
+
+Iterated filtering is a technique for maximizing the likelihood obtained by filtering.
+In `pomp`, it is the particle filter that is iterated.
+Iterated filtering is implemented in the `mif` function.
+ at Ionides2015 describe an improvement on the original [@Ionides2006] algorithm.
+This "IF2" algorithm is currently implemented by setting the `method="mif2"` option in `mif`.
+
+Let's use iterated filtering to obtain an approximate MLE for the data in the Gompertz model example provided with `pomp`.
+We'll initialize the algorithm at several starting points around `theta.true` and just estimate the parameters $r$, $\tau$, and $\sigma$:
+```{r gompertz-init,cache=FALSE}
+require(pomp)
+pompExample(gompertz)
+theta <- coef(gompertz)
+theta.true <- theta
+```
+
+Note that we've set `transform=TRUE` in the above.
+This means that the likelihood maximization is done on the estimation scale (see the section above on Parameter Transformations).
+
+```{r gompertz-multi-mif2-eval,results='hide'}
+require(foreach)
+require(doMC)
+registerDoMC()
+
+save.seed <- .Random.seed
+set.seed(334388458L,kind="L'Ecuyer")
+
+estpars <- c("r","sigma","tau")
+mf <- foreach(i=1:10,
+ .inorder=FALSE,
+ .options.multicore=list(set.seed=TRUE)
+ ) %dopar%
+ {
+ theta.guess <- theta.true
+ theta.guess[estpars] <- rlnorm(
+ n=length(estpars),
+ meanlog=log(theta.guess[estpars]),
+ sdlog=1
+ )
+ m1 <- mif(
+ gompertz,
+ Nmif=50,
+ method="mif2",
+ start=theta.guess,
+ transform=TRUE,
+ rw.sd=c(r=0.02,sigma=0.02,tau=0.05),
+ Np=2000,
+ var.factor=2,
+ cooling.type="hyperbolic",
+ cooling.fraction=0.95
+ )
+ m1 <- continue(m1,Nmif=50,cooling.fraction=0.8)
+ m1 <- continue(m1,Nmif=50,cooling.fraction=0.6)
+ m1 <- continue(m1,Nmif=50,cooling.fraction=0.2)
+ ll <- replicate(n=10,logLik(pfilter(m1,Np=10000)))
+ list(mif=m1,ll=ll)
+ }
+```
+```{r gompertz-post-mif2}
+loglik.true <- replicate(n=10,logLik(pfilter(gompertz,Np=10000)))
+loglik.true <- logmeanexp(loglik.true,se=TRUE)
+theta.mif <- t(sapply(mf,function(x)coef(x$mif)))
+loglik.mif <- t(sapply(mf,function(x)logmeanexp(x$ll,se=TRUE)))
+best <- which.max(loglik.mif[,1])
+theta.mif <- theta.mif[best,]
+loglik.mif <- loglik.mif[best,]
+rbind(
+ mle=c(signif(theta.mif[estpars],3),loglik=round(loglik.mif,2)),
+ truth=c(signif(theta.true[estpars],3),loglik=round(loglik.true,2))
+ ) -> results.table
+```
+
+Each of the `r length(mf)` `mif` runs ends up at a different place.
+In this case (but by no means in every case), we can average across these parameter estimates to get an approximate maximum likelihood estimate.
+We'll evaluate the likelihood several times at this estimate to get an idea of the Monte Carlo error in our likelihood estimate.
+The particle filter produces an unbiased estimate of the likelihood;
+therefore, we'll average the likelihoods, not the log likelihoods.
+
+
+Convergence plots can be used to help diagnose convergence of the iterated filtering algorithm.
+The following shows part of the output of `plot(mf)`.
+
+```{r mif2-plot,echo=FALSE,cache=FALSE,fig.height=6}
+op <- par(mfrow=c(4,1),mar=c(3,3,0,0),mgp=c(2,1,0),bty='l')
+loglik <- sapply(mf,function(x)conv.rec(x$mif,"loglik"))
+log.r <- sapply(mf,function(x)conv.rec(x$mif,"r"))
+log.sigma <- sapply(mf,function(x)conv.rec(x$mif,"sigma"))
+log.tau <- sapply(mf,function(x)conv.rec(x$mif,"tau"))
+matplot(loglik,type='l',lty=1,xlab="",ylab=expression(log~L),xaxt='n',ylim=max(loglik,na.rm=T)+c(-12,3))
+matplot(log.r,type='l',lty=1,xlab="",ylab=expression(log~r),xaxt='n')
+matplot(log.sigma,type='l',lty=1,xlab="",ylab=expression(log~sigma),xaxt='n')
+matplot(log.tau,type='l',lty=1,xlab="MIF iteration",ylab=expression(log~tau))
+par(op)
+```
+
+Here is a summary of the results.
+```{r first-mif-results-table,echo=FALSE,cache=FALSE}
+print(results.table)
+```
+
+## References
Added: www/vignettes/mif2.html
===================================================================
--- www/vignettes/mif2.html (rev 0)
+++ www/vignettes/mif2.html 2015-03-23 12:59:09 UTC (rev 1145)
@@ -0,0 +1,162 @@
+<!DOCTYPE html>
+
+<html xmlns="http://www.w3.org/1999/xhtml">
+
+<head>
+
+<meta charset="utf-8">
+<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
+<meta name="generator" content="pandoc" />
+
+<meta name="author" content="Aaron A. King" />
+
+
+<title>An IF2 example</title>
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 1145
More information about the pomp-commits
mailing list