[Pomp-commits] r1179 - pkg/pomp pkg/pomp/R pkg/pomp/tests www/vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jun 4 21:39:18 CEST 2015
Author: kingaa
Date: 2015-06-04 21:39:17 +0200 (Thu, 04 Jun 2015)
New Revision: 1179
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/R/mif2.R
pkg/pomp/tests/ou2-mif2.R
pkg/pomp/tests/ou2-mif2.Rout.save
www/vignettes/mif2.R
www/vignettes/mif2.Rmd
www/vignettes/mif2.html
www/vignettes/pomp.pdf
Log:
- fix bug in 'rw.sd'
- update 'mif2' vignette
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2015-06-04 18:07:50 UTC (rev 1178)
+++ pkg/pomp/DESCRIPTION 2015-06-04 19:39:17 UTC (rev 1179)
@@ -1,7 +1,7 @@
Package: pomp
Type: Package
Title: Statistical Inference for Partially Observed Markov Processes
-Version: 0.66-2
+Version: 0.66-3
Date: 2015-06-04
Authors at R: c(person(given=c("Aaron","A."),family="King",
role=c("aut","cre"),email="kingaa at umich.edu"),
Modified: pkg/pomp/R/mif2.R
===================================================================
--- pkg/pomp/R/mif2.R 2015-06-04 18:07:50 UTC (rev 1178)
+++ pkg/pomp/R/mif2.R 2015-06-04 19:39:17 UTC (rev 1179)
@@ -18,10 +18,13 @@
sds <- lapply(rw.sd,eval,envir=list(time=time,ivp=ivp))
for (n in names(sds)) {
len <- length(sds[[n]])
- if (len!=1 && len!=length(time))
+ if (len==1) {
+ sds[[n]] <- rep(sds[[n]],length(time))
+ } else if (len!=length(time)) {
stop(sQuote("mif2")," error: ",sQuote("rw.sd")," spec for parameter ",
sQuote(n)," does not evaluate to a vector of the correct length (",
length(time),")",call.=FALSE)
+ }
}
do.call(rbind,sds)
}
Modified: pkg/pomp/tests/ou2-mif2.R
===================================================================
--- pkg/pomp/tests/ou2-mif2.R 2015-06-04 18:07:50 UTC (rev 1178)
+++ pkg/pomp/tests/ou2-mif2.R 2015-06-04 19:39:17 UTC (rev 1179)
@@ -107,17 +107,17 @@
guess1[c('x1.0','x2.0','alpha.2','alpha.3')] <- 0.5*guess1[c('x1.0','x2.0','alpha.2','alpha.3')]
guess2[c('x1.0','x2.0','alpha.2','alpha.3')] <- 1.5*guess1[c('x1.0','x2.0','alpha.2','alpha.3')]
-m1 <- pomp:::mif2(ou2,Nmif=100,start=guess1,Np=1000,
- cooling.type="hyperbolic",cooling.fraction.50=0.05,
- rw.sd=pomp:::rw.sd(
- x1.0=ivp(0.5),x2.0=ivp(0.5),
- alpha.2=0.1,alpha.3=0.1))
+m1 <- mif2(ou2,Nmif=100,start=guess1,Np=1000,
+ cooling.type="hyperbolic",cooling.fraction.50=0.05,
+ rw.sd=rw.sd(
+ x1.0=ivp(0.5),x2.0=ivp(0.5),
+ alpha.2=0.1,alpha.3=0.1))
-m2 <- pomp:::mif2(ou2,Nmif=100,start=guess2,Np=1000,
- cooling.type="hyperbolic",cooling.fraction.50=0.05,
- rw.sd=pomp:::rw.sd(
- x1.0=ivp(0.5),x2.0=ivp(0.5),
- alpha.2=0.1,alpha.3=0.1))
+m2 <- mif2(ou2,Nmif=100,start=guess2,Np=1000,
+ cooling.type="hyperbolic",cooling.fraction.50=0.05,
+ rw.sd=rw.sd(
+ x1.0=ivp(0.5),x2.0=ivp(0.5),
+ alpha.2=0.1,alpha.3=0.1))
plot(c(m1,m2))
@@ -125,4 +125,12 @@
mle2=c(coef(m2),loglik=logLik(pfilter(m1,Np=1000))),
truth=c(coef(ou2),loglik=logLik(pfilter(m1,Np=1000))))
+m3 <- mif2(ou2,Nmif=3,start=guess1,Np=200,
+ cooling.fraction=0.2,
+ rw.sd=rw.sd(
+ x1.0=c(0.5,rep(0.2,99)),
+ x2.0=ivp(0.5),
+ alpha.2=if (time==1) 0.2 else 0.1,
+ alpha.3=0.2*(time<10)))
+
dev.off()
Modified: pkg/pomp/tests/ou2-mif2.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-mif2.Rout.save 2015-06-04 18:07:50 UTC (rev 1178)
+++ pkg/pomp/tests/ou2-mif2.Rout.save 2015-06-04 19:39:17 UTC (rev 1179)
@@ -138,17 +138,17 @@
> guess1[c('x1.0','x2.0','alpha.2','alpha.3')] <- 0.5*guess1[c('x1.0','x2.0','alpha.2','alpha.3')]
> guess2[c('x1.0','x2.0','alpha.2','alpha.3')] <- 1.5*guess1[c('x1.0','x2.0','alpha.2','alpha.3')]
>
-> m1 <- pomp:::mif2(ou2,Nmif=100,start=guess1,Np=1000,
-+ cooling.type="hyperbolic",cooling.fraction.50=0.05,
-+ rw.sd=pomp:::rw.sd(
-+ x1.0=ivp(0.5),x2.0=ivp(0.5),
-+ alpha.2=0.1,alpha.3=0.1))
+> m1 <- mif2(ou2,Nmif=100,start=guess1,Np=1000,
++ cooling.type="hyperbolic",cooling.fraction.50=0.05,
++ rw.sd=rw.sd(
++ x1.0=ivp(0.5),x2.0=ivp(0.5),
++ alpha.2=0.1,alpha.3=0.1))
>
-> m2 <- pomp:::mif2(ou2,Nmif=100,start=guess2,Np=1000,
-+ cooling.type="hyperbolic",cooling.fraction.50=0.05,
-+ rw.sd=pomp:::rw.sd(
-+ x1.0=ivp(0.5),x2.0=ivp(0.5),
-+ alpha.2=0.1,alpha.3=0.1))
+> m2 <- mif2(ou2,Nmif=100,start=guess2,Np=1000,
++ cooling.type="hyperbolic",cooling.fraction.50=0.05,
++ rw.sd=rw.sd(
++ x1.0=ivp(0.5),x2.0=ivp(0.5),
++ alpha.2=0.1,alpha.3=0.1))
>
> plot(c(m1,m2))
>
@@ -164,10 +164,21 @@
mle2 -2.696841 3.157511 -482.3154
truth -3.000000 4.000000 -481.9939
>
+> m3 <- mif2(ou2,Nmif=3,start=guess1,Np=200,
++ cooling.fraction=0.2,
++ rw.sd=rw.sd(
++ x1.0=c(0.5,rep(0.2,99)),
++ x2.0=ivp(0.5),
++ alpha.2=if (time==1) 0.2 else 0.1,
++ alpha.3=0.2*(time<10)))
+Warning message:
+In if (time == 1) 0.2 else 0.1 :
+ the condition has length > 1 and only the first element will be used
+>
> dev.off()
null device
1
>
> proc.time()
user system elapsed
- 70.828 0.064 71.333
+ 69.328 0.052 69.748
Modified: www/vignettes/mif2.R
===================================================================
--- www/vignettes/mif2.R 2015-06-04 18:07:50 UTC (rev 1178)
+++ www/vignettes/mif2.R 2015-06-04 19:39:17 UTC (rev 1179)
@@ -5,7 +5,7 @@
require(magrittr)
theme_set(theme_bw())
require(pomp)
-stopifnot(packageVersion("pomp")>="0.63-1")
+stopifnot(packageVersion("pomp")>="0.66-2")
options(
keep.source=TRUE,
stringsAsFactors=FALSE,
@@ -40,17 +40,14 @@
meanlog=log(theta.guess[estpars]),
sdlog=1
)
- m1 <- mif(
+ m1 <- mif2(
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
+ rw.sd=rw.sd(r=0.02,sigma=0.02,tau=0.05),
+ cooling.fraction.50=0.95,
+ Np=2000
)
m1 <- continue(m1,Nmif=50,cooling.fraction=0.8)
m1 <- continue(m1,Nmif=50,cooling.fraction=0.6)
Modified: www/vignettes/mif2.Rmd
===================================================================
--- www/vignettes/mif2.Rmd 2015-06-04 18:07:50 UTC (rev 1178)
+++ www/vignettes/mif2.Rmd 2015-06-04 19:39:17 UTC (rev 1179)
@@ -48,7 +48,7 @@
require(magrittr)
theme_set(theme_bw())
require(pomp)
-stopifnot(packageVersion("pomp")>="0.65-1")
+stopifnot(packageVersion("pomp")>="0.66-2")
options(
keep.source=TRUE,
stringsAsFactors=FALSE,
@@ -60,12 +60,12 @@
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.
+The iterated filtering of @Ionides2006 is implemented in the `mif` function.
@Ionides2015 describe an improvement on the original [@Ionides2006] algorithm.
-This "IF2" algorithm is currently implemented by setting the `method="mif2"` option in `mif`.
+This "IF2" algorithm is implemented in the `mif2` function.
-Let's use IF2 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$:
+The following constructs the Gompertz example that is provided with `pomp` (see `?gompertz`) and extracts the parameters at which the data were generated.
+
```{r gompertz-init,cache=FALSE}
require(pomp)
pompExample(gompertz)
@@ -73,8 +73,8 @@
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).
+Let's use IF2 to obtain an approximate MLE for these data.
+We'll initialize the algorithm at several starting points around `theta.true` and just estimate the parameters $r$, $\tau$, and $\sigma$:
```{r gompertz-multi-mif2-eval,results='hide'}
require(foreach)
@@ -96,17 +96,14 @@
meanlog=log(theta.guess[estpars]),
sdlog=1
)
- m1 <- mif(
+ m1 <- mif2(
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
+ rw.sd=rw.sd(r=0.02,sigma=0.02,tau=0.05),
+ cooling.fraction.50=0.95,
+ Np=2000
)
m1 <- continue(m1,Nmif=50,cooling.fraction=0.8)
m1 <- continue(m1,Nmif=50,cooling.fraction=0.6)
@@ -115,6 +112,15 @@
list(mif=m1,ll=ll)
}
```
+
+Note that we've set `transform=TRUE` in the above to search for the MLE with the parameters transformed to enforce their positivity.
+See the `pomp` documentation (`?pomp`) and the section on Parameter Transformations in the [Getting Started vignette](http://pomp.r-forge.r-project.org/vignettes/getting_started.html)).
+
+Each of the `r length(mf)` `mif2` runs ends up at a different point estimate.
+We focus on that with the highest estimated likelihood, having evaluated the likelihood several times to reduce the Monte Carlo error in the likelihood evaluation.
+The particle filter produces an unbiased estimate of the likelihood;
+therefore, we will average the likelihoods, not the log likelihoods.
+
```{r gompertz-post-mif2}
loglik.true <- replicate(n=10,logLik(pfilter(gompertz,Np=10000)))
loglik.true <- logmeanexp(loglik.true,se=TRUE)
@@ -129,15 +135,8 @@
) -> 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)`.
+Something like the following can be obtained by executing `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')
@@ -156,5 +155,6 @@
```{r first-mif-results-table,echo=FALSE,cache=FALSE}
print(results.table)
```
+The fact that the likelihood at the putative MLE is higher than it is at the truth is evidence that this parameter is actually close to the true MLE.
## References
Modified: www/vignettes/mif2.html
===================================================================
--- www/vignettes/mif2.html 2015-06-04 18:07:50 UTC (rev 1178)
+++ www/vignettes/mif2.html 2015-06-04 19:39:17 UTC (rev 1179)
@@ -67,16 +67,16 @@
<p>Licensed under the <a href="http://creativecommons.org/licenses/by-nc/3.0">Creative Commons attribution-noncommercial license</a>. Please share and remix noncommercially, mentioning its origin.<br /><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFgAAAAfCAYAAABjyArgAAAACXBIWXMAAAsTAAALEwEAmpwYAAAKT2lDQ1BQaG90b3Nob3AgSUNDIHByb2ZpbGUAAHjanVNnVFPpFj333vRCS4iAlEtvUhUIIFJCi4AUkSYqIQkQSoghodkVUcERRUUEG8igiAOOjoCMFVEsDIoK2AfkIaKOg6OIisr74Xuja9a89+bN/rXXPues852zzwfACAyWSDNRNYAMqUIeEeCDx8TG4eQuQIEKJHAAEAizZCFz/SMBAPh+PDwrIsAHvgABeNMLCADATZvAMByH/w/qQplcAYCEAcB0kThLCIAUAEB6jkKmAEBGAYCdmCZTAKAEAGDLY2LjAFAtAGAnf+bTAICd+Jl7AQBblCEVAaCRACATZYhEAGg7AKzPVopFAFgwABRmS8Q5ANgtADBJV2ZIALC3AMDOEAuyAAgMADBRiIUpAAR7AGDIIyN4AISZABRG8lc88SuuEOcqAAB4mbI8uSQ5RYFbCC1xB1dXLh4ozkkXKxQ2YQJhmkAuwnmZGTKBNA/g88wAAKCRFRHgg/P9eM4Ors7ONo62Dl8t6r8G/yJiYuP+5c+rcEAAAOF0ftH+LC+zGoA7BoBt/qIl7gRoXgugdfeLZrIPQLUAoOnaV/Nw+H48PEWhkLnZ2eXk5NhKxEJbYcpXff5nwl/AV/1s+X48/Pf14L7iJIEyXYFHBPjgwsz0TKUcz5IJhGLc5o9H/LcL//wd0yLESWK5WCoU41EScY5EmozzMqUiiUKSKcUl0v9k4t8s+wM+3zUAsGo+AXuRLahdYwP2SycQWHTA4vcAAPK7b8HUKAgDgGiD4c93/+8//UegJQCAZkmScQAAXkQkLlTKsz/HCAAARKCBKrBBG/TBGCzABhzBBdzBC/xgNoRCJMTCQhBCCmSAHHJgKayCQiiGzbAdKmAv1EAdNMBRaIaTcA4uwlW4Dj1wD/phCJ7BKLyBCQRByAgTYSHaiAFiilgjjggXmYX4IcFIBBKLJCDJiBRRIkuRNUgxUopUIFVIHfI9cgI5h1xGupE7yAAygvyGvEcxlIGyUT3UDLVDuag3GoRGogvQZHQxmo8WoJvQcrQaPYw2oefQq2gP2o8+Q8cwwOgYBzPEbDAuxsNCsTgsCZNjy7EirAyrxhqwVqwDu4n1Y8+xdwQSgUXACTYEd0IgYR5BSFhMWE7YSKggHCQ0EdoJNwkDhFHCJyKTqEu0JroR+cQYYjIxh1hILCPWEo8TLxB7iEPENyQSiUMyJ7mQAkmxpFTSEtJG0m5SI+ksqZs0SBojk8naZGuyBzmULCAryIXkneTD5DPkG+Qh8lsKnWJAcaT4U+IoUspqShnlEOU05QZlmDJBVaOaUt2ooVQRNY9aQq2htlKvUYeoEzR1mjnNgxZJS6WtopXTGmgXaPdpr+h0uhHdlR5Ol9BX0svpR+iX6AP0dwwNhhWDx4hnKBmbGAcYZxl3GK+YTKYZ04sZx1QwNzHrmOeZD5lvVVgqtip8FZHKCpVKlSaVGyovVKmqpqreqgtV81XLVI+pXlN9rkZVM1PjqQnUlqtVqp1Q61MbU2epO6iHqmeob1Q/pH5Z/YkGWcNMw09DpFGgsV/jvMYgC2MZs3gsIWsNq4Z1gTXEJrHN2Xx2KruY/R27iz2qqaE5QzNKM1ezUvOUZj8H45hx+Jx0TgnnKKeX836K3hTvKeIpG6Y0TLkxZVxrqpaXllirSKtRq0frvTau7aedpr1Fu1n7gQ5Bx0onXCdHZ4/OBZ3nU9lT3acKpxZNPTr1ri6qa6UbobtEd79up+6Ynr5egJ5Mb6feeb3n+hx9L/1U/W36p/VHDFgGswwkBtsMzhg8xTVxbzwdL8fb8VFDXcNAQ6VhlWGX4YSRudE8o9VGjUYPjGnGXOMk423GbcajJgYmISZLTepN7ppSTbmmKaY7TDtMx83MzaLN1pk1mz0x1zLnm+eb15vft2BaeFostqi2uGVJsuRaplnutrxuhVo5WaVYVVpds0atna0l1rutu6cRp7lOk06rntZnw7Dxtsm2qbcZsOXYBtuutm22fWFnYhdnt8Wuw+6TvZN9un2N/T0HDYfZDqsdWh1+c7RyFDpWOt6azpzuP33F9JbpL2dYzxDP2DPjthPLKcRpnVOb00dnF2e5c4PziIuJS4LLLpc+Lpsbxt3IveRKdPVxXeF60vWdm7Obwu2o26/uNu5p7ofcn8w0nymeWTNz0MPIQ+BR5dE/C5+VMGvfrH5PQ0+BZ7XnIy9jL5FXrdewt6V3qvdh7xc+9j5yn+M+4zw33jLeWV/MN8C3yLfLT8Nvnl+F30N/I/9k/3r/0QCngCUBZwOJgUGBWwL7+Hp8Ib+OPzrbZfay2e1BjKC5QRVBj4KtguXBrSFoyOyQrSH355jOkc5pDoVQfujW0Adh5mGLw34MJ4WHhVeGP45wiFga0TGXNXfR3ENz30T6RJZE3ptnMU85ry1KNSo+qi5qPNo3ujS6P8YuZlnM1VidWElsSxw5LiquNm5svt/87fOH4p3iC+N7F5gvyF1weaHOwvSFpxapLhIsOpZATIhOOJTwQRAqqBaMJfITdyWOCnnCHcJnIi/RNtGI2ENcKh5O8kgqTXqS7JG8NXkkxTOlLOW5hCepkLxMDUzdmzqeFpp2IG0yPTq9MYOSkZBxQqohTZO2Z+pn5mZ2y6xlhbL+xW6Lty8elQfJa7OQrAVZLQq2QqboVFoo1yoHsmdlV2a/zYnKOZarnivN7cyzytuQN5zvn//tEsIS4ZK2pYZLVy0dWOa9rGo5sjxxedsK4xUFK4ZWBqw8uIq2Km3VT6vtV5eufr0mek1rgV7ByoLBtQFr6wtVCuWFfevc1+1dT1gvWd+1YfqGnRs+FYmKrhTbF5cVf9go3HjlG4dvyr+Z3JS0qavEuWTPZtJm6ebeLZ5bDpaql+aXDm4N2dq0Dd9WtO319kXbL5fNKNu7g7ZDuaO/PLi8ZafJzs07P1SkVPRU+lQ27tLdtWHX+G7R7ht7vPY07NXbW7z3/T7JvttVAVVN1WbVZftJ+7P3P66Jqun4lvttXa1ObXHtxwPSA/0HIw6217nU1R3SPVRSj9Yr60cOxx++/p3vdy0NNg1VjZzG4iNwRHnk6fcJ3/ceDTradox7rOEH0x92HWcdL2pCmvKaRptTmvtbYlu6T8w+0dbq3nr8R9sfD5w0PFl5SvNUyWna6YLTk2fyz4ydlZ19fi753GDborZ752PO32oPb++6EHTh0kX/i+c7vDvOXPK4dPKy2+UTV7hXmq86X23qdOo8/pPTT8e7nLuarrlca7nuer21e2b36RueN87d9L158Rb/1tWeOT3dvfN6b/fF9/XfFt1+cif9zsu72Xcn7q28T7xf9EDtQdlD3YfVP1v+3Njv3H9qwHeg89HcR/cGhYPP/pH1jw9DBY+Zj8uGDYbrnjg+OTniP3L96fynQ89kzyaeF/6i/suuFxYvfvjV69fO0ZjRoZfyl5O/bXyl/erA6xmv28bCxh6+yXgzMV70VvvtwXfcdx3vo98PT+R8IH8o/2j5sfVT0Kf7kxmTk/8EA5jz/GMzLdsAAAAEZ0FNQQAAsY58+1GTAAAAIGNIUk0AAHolAACAgwAA+f8AAIDpAAB1MAAA6mAAADqYAAAXb5JfxUYAAAk0SURBVHja7FpdbNvWFf5IysFS1BrztA1yMBt7sQqskZMmy4Ytlta9LJ4TCnaCFkkWuQ812mCTlB+3S+3Iyk8TK/Zkb0iBYVstrCjahwZm/oDNGSLaKzBbTiIZaSM9rJCK2FiHDbArpwVmkbx7EHlF2pIty3axpjnGFX/uvR/J75577jnnmiGEWBmG+RSPZc2FEMIwAAgA3Bi+DpZlwbIsOI4Dy3LgWBYspx1ZFgzDgmUYMAwDMIyOAgICohKoRIWq5ouiKPmjqkBRVKiqQutUotL2hBD9Zej5oyD79u4HADAAiE4ux3H5wnKFc47L17GcRjIDhmGN/GrkaMSqeTIVRSvGc8VMsqqqlFgj0Y8SyRYAZnI5CyymY75cu3Id0WgUsVgMc9k5E1C1tRo7duyA68cuNO/5GRSVA8vK+QFRWDBgtLE0TB+V5GcCtDnELE3u3Yk4xMsiksmk6b7dbofQImDr9oZVkbFe+AwA8pdbf4bFYqGkWiyWfOEsGJFGEboQwvT0dFmANpsNHb/qQGPjLsiKAkWRIctaUWTIskI1vJgmL9TiT/75L1wauIRUMgUAcDqdcDgcAIBEIgFJkgAA9fZ6HPEewTe/9Y0VEbCe+Pv27s8T/NeRm7BwlgKxlipUWSwIdHVDHBJpB57nIQgCamtr0djYCAAYGRlBJpOBKIqYnZ2lbQW3gOMdx7DxiY2QZRk5WYYs5yDLstlk6ASrmi03EP0w+xDeIz5ks1kIgoBwOIza2lrTR2QyGfj9foiiCKvVinOhc2WTsN74lOBbf7uFKp3YKguqLFUmcnmeh8/ng9frBc/zJQEjkQj8fj8lut5ejz+8+Xt8beNGyHIOuVyOanJRkhfY465XTyGVTMHj8WBwcLAw7TTTYtT0SCSCtrY21NvrcebC6bIIKIX/m/5+jI+N4+1331kV/r69+8ECAKfZYAvHwcKZNdfhcCAejyMQCCxJLgB4PB6k02k6xVLJFHpDfSZbzlKPhCkU7c9I4N2JOFLJFARBQMeJE8t+jMfjgSAISCVTuDsRL8vmppIpbG1owA92ft9E7oVQCNdu3MArx09gamqqInxdWABgdbeM4zAijZrIjUaji6bNUsLzPKLRKCVZHBIxMjIKjrMUJ5nV3T6YSBYv598hHA7D/tRTC/3LogtiOBw29V1K9DafP/wMPefPw/nDH+GlF9vh9fvR3t6OkydPItTXi/GxsYrwTQTnNTivxaELIUrU4ODgslq7FMl639D5kOZLa+Qymk/Nah4GgwLJ2vRPJpNwOp2LBretrY1qfltbm6mutrYWTqdzkSdQTHT85uZm7Nu/H1NTU7g5PIzvfLsWn889xMFDB3H/ww/R0tpaEb5Zg7WPv3blOvUWfD4f1cJKhOd5OuLT09O4dvU6DVjyJJc2EboUe34kEil6vlSfUuJwOBDq68X5UA/efvcdtLS24qOPMwj19WLz5s2IvDmI5P37FeNTgnVtikajlByv10sbSZIEt9sNl8sFl8uFYDBYsq6/v99kF3Utjt6KGrS3YBoYpriJ+KLlezt3oqf3Ih48eICOY8fR8N2ncfm999C8uwkHnnseN4eHK8LNBxoMA5ZhEIvF8i6WIFBiJEmCy+UydZIkCZIkwev1wu12L6qbnJykq7IgCIhEIojFYvkI0EAsUCC34JUXsBKJRNHFTNdcj8ezqL5Yn1KysG02m8XN4WH09F6E534bmnc3AQDGx8YwPjaGmpoaMFWWSjQ4/6F6hLZlyxbawO/3U/uTTqfponf48GGqyQ6HA+l0GkNDQyYfGQA9n8vOFcwBPeq8LjYRdrsdkiQhk8mY7hvdKeO57rNKkgS73b7shxfDf+nFdpw7fQZbn96CA889j48+zqCltRU9vRdx4ODBFeGbCDYuLgvtjD7KHo+HGvl0Og2Px0Pr9OBDEARaZ1wYCu4X/Vn2xYQWwTTA5YjeVu+7Uvye3otoe+EFfPKff+Mf6TQGwmG8dqoLLa2tCJ49g4btz5SNbyb4/1C2bm9Avb0eoigu8hZKkSuKIurt9WXlDYrh19TU4LVTXTjmP4rmpib80ueD1WqtCN9MMDFHRUbbpGtzJBLB7OwsEokE6urqEAwGC76uFiYb64zTtuC/0p+yXu6Vkx2wWq2IRCJwu90Uy+gHZzIZuN1u9Pf3w2q14oj3SNkfXwr/2InjNIpbDT5d5PQXrrZWYy47h8nJSdogEAjQh2/atMlEnF6XSCQW1emiY1Vbq0GIIRwGgT6m2tWil3vS+iQGLvWj5/UQRFGEKIpwOBxwOBzgeR6SJFFlqCQZs974dN0evzOODRs2IHgqCHFIBM/zmJmZMXkGAwMDNMfgdDoRCASo9g4MDNC2xjoAqKurQyaTwbM/eRZ94V78d34eudw8cnIOuZwMRc9LqApURaVJ+IWR2pcxXUmTPWO3/46qqg14f/R9eH/hpWGhz+db1UvrCRIAOPv6Wexu+inmc/OYz+Ug53LIyQsILpLw+bIn3FkAdEvH6WqEzWYDAASDwUUu0kpkdnaWrtA2m01LxOupSRWqZot1c/GoCgsAKlGhKPl9tDPnTlOC3G63Kce7EnJdLhft2/Fqh5ZgVwokG1KUdJElj9aWEV3cAZDOQCeZuBsjiXtxcujnh/SlnjgcDhKPx0m5kk6nicPhoP0Ft0AS9+Jk4m6MnD53mt7/CpXChe+ol9yOT5DEBwkiuAV6n+d50t3dTWZmZkoSOzMzQ7q7uwnP87Tf9h3byeQHCXI7MUE6A51fRXIJs9Ap7Qp0Yq9bgMViQV+oD2/96a2Kt4yCZ7ohyzLGx2N4uf3lyqeYwWwYdzOM0efC65Xil8LSn10pNoqx3hXozGvyvTh5M/JHYrPZyh4xm81GBn47QBL34uR2YoK88bs3Vq0FhJAlz433KsVfDrfSZzClwirfUS8OHDxAdyGk6AiuXrm65Lb9HmEPnM5Gutd25cpVnA2eXf0iUUSD10JzF2KUOq5GmKXi1q3bGtAVOIWazTbDboQx3QiT36r/48n01BR+3RvG6Mjo2qzCC6bsWpmG5UzCUs9dE4J12dW4C03NTdj2zDZ83Wo153A11+rTbBZ3bt/BjWs31ozYL1qD18MGl0XwY1mFiSCEPMEwzGePqViHAIMQ5n8DAFb/49reYmyHAAAAAElFTkSuQmCC" alt="CC-BY_NC" /></p>
-<p>This document was produced using <code>pomp</code> version 0.65.1.</p>
-<p>Iterated filtering is a technique for maximizing the likelihood obtained by filtering. In <code>pomp</code>, it is the particle filter that is iterated. Iterated filtering is implemented in the <code>mif</code> function. <span class="citation">Ionides et al. (2015)</span> describe an improvement on the original <span class="citation">(Ionides et al. 2006)</span> algorithm. This “IF2” algorithm is currently implemented by setting the <code>method="mif2"</code> option in <code>mif</code>.</p>
-<p>Let’s use IF2 to obtain an approximate MLE for the data in the Gompertz model example provided with <code>pomp</code>. We’ll initialize the algorithm at several starting points around <code>theta.true</code> and just estimate the parameters <span class="math">\(r\)</span>, <span class="math">\(\tau\)</span>, and <span class="math">\(\sigma\)</span>:</p>
+<p>This document was produced using <code>pomp</code> version 0.66.3.</p>
+<p>Iterated filtering is a technique for maximizing the likelihood obtained by filtering. In <code>pomp</code>, it is the particle filter that is iterated. The iterated filtering of <span class="citation">Ionides et al. (2006)</span> is implemented in the <code>mif</code> function. <span class="citation">Ionides et al. (2015)</span> describe an improvement on the original <span class="citation">(Ionides et al. 2006)</span> algorithm. This “IF2” algorithm is implemented in the <code>mif2</code> function.</p>
+<p>The following constructs the Gompertz example that is provided with <code>pomp</code> (see <code>?gompertz</code>) and extracts the parameters at which the data were generated.</p>
<pre class="r"><code>require(pomp)
pompExample(gompertz)</code></pre>
<pre><code>## newly created object(s):
## gompertz</code></pre>
<pre class="r"><code>theta <- coef(gompertz)
theta.true <- theta</code></pre>
-<p>Note that we’ve set <code>transform=TRUE</code> in the above. This means that the likelihood maximization is done on the estimation scale (see the section above on Parameter Transformations).</p>
+<p>Let’s use IF2 to obtain an approximate MLE for these data. We’ll initialize the algorithm at several starting points around <code>theta.true</code> and just estimate the parameters <span class="math">\(r\)</span>, <span class="math">\(\tau\)</span>, and <span class="math">\(\sigma\)</span>:</p>
<pre class="r"><code>require(foreach)
require(doMC)
registerDoMC()
@@ -96,17 +96,14 @@
meanlog=log(theta.guess[estpars]),
sdlog=1
)
- m1 <- mif(
+ m1 <- mif2(
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
+ rw.sd=rw.sd(r=0.02,sigma=0.02,tau=0.05),
+ cooling.fraction.50=0.95,
+ Np=2000
)
m1 <- continue(m1,Nmif=50,cooling.fraction=0.8)
m1 <- continue(m1,Nmif=50,cooling.fraction=0.6)
@@ -114,6 +111,8 @@
ll <- replicate(n=10,logLik(pfilter(m1,Np=10000)))
list(mif=m1,ll=ll)
}</code></pre>
+<p>Note that we’ve set <code>transform=TRUE</code> in the above to search for the MLE with the parameters transformed to enforce their positivity. See the <code>pomp</code> documentation (<code>?pomp</code>) and the section on Parameter Transformations in the <a href="http://pomp.r-forge.r-project.org/vignettes/getting_started.html">Getting Started vignette</a>).</p>
+<p>Each of the 10 <code>mif2</code> runs ends up at a different point estimate. We focus on that with the highest estimated likelihood, having evaluated the likelihood several times to reduce the Monte Carlo error in the likelihood evaluation. The particle filter produces an unbiased estimate of the likelihood; therefore, we will average the likelihoods, not the log likelihoods.</p>
<pre class="r"><code>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)))
@@ -125,13 +124,13 @@
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</code></pre>
-<p>Each of the 10 <code>mif</code> 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.</p>
-<p>Convergence plots can be used to help diagnose convergence of the iterated filtering algorithm. The following shows part of the output of <code>plot(mf)</code>.</p>
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 1179
More information about the pomp-commits
mailing list