[Pomp-commits] r1182 - in pkg/pomp: . R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 5 01:29:30 CEST 2015
Author: kingaa
Date: 2015-06-05 01:29:29 +0200 (Fri, 05 Jun 2015)
New Revision: 1182
Added:
pkg/pomp/R/mif2-methods.R
pkg/pomp/man/mif2.Rd
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/R/generics.R
pkg/pomp/R/mif-methods.R
pkg/pomp/R/mif2.R
pkg/pomp/man/mif.Rd
pkg/pomp/man/pfilter.Rd
pkg/pomp/man/pmcmc.Rd
pkg/pomp/tests/ou2-mif2.R
pkg/pomp/tests/ou2-mif2.Rout.save
Log:
- improve documentation
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2015-06-04 23:29:13 UTC (rev 1181)
+++ pkg/pomp/DESCRIPTION 2015-06-04 23:29:29 UTC (rev 1182)
@@ -1,7 +1,7 @@
Package: pomp
Type: Package
Title: Statistical Inference for Partially Observed Markov Processes
-Version: 0.66-3
+Version: 0.66-4
Date: 2015-06-04
Authors at R: c(person(given=c("Aaron","A."),family="King",
role=c("aut","cre"),email="kingaa at umich.edu"),
@@ -34,7 +34,7 @@
simulate-pomp.R trajectory-pomp.R plot-pomp.R
pfilter.R pfilter-methods.R minim.R traj-match.R
bsmc.R bsmc2.R
- mif.R mif-methods.R mif2.R
+ mif.R mif-methods.R mif2.R mif2-methods.R
proposals.R pmcmc.R pmcmc-methods.R
nlf-funcs.R nlf-guts.R nlf-objfun.R nlf.R
probe.R probe-match.R basic-probes.R spect.R spect-match.R
Modified: pkg/pomp/R/generics.R
===================================================================
--- pkg/pomp/R/generics.R 2015-06-04 23:29:13 UTC (rev 1181)
+++ pkg/pomp/R/generics.R 2015-06-04 23:29:29 UTC (rev 1182)
@@ -74,6 +74,7 @@
## iterated filtering
setGeneric('mif',function(object,...)standardGeneric("mif"))
+setGeneric("mif2",function(object,...)standardGeneric("mif2"))
## generate new particles
setGeneric('particles',function(object,...)standardGeneric("particles"))
Modified: pkg/pomp/R/mif-methods.R
===================================================================
--- pkg/pomp/R/mif-methods.R 2015-06-04 23:29:13 UTC (rev 1181)
+++ pkg/pomp/R/mif-methods.R 2015-06-04 23:29:29 UTC (rev 1182)
@@ -40,10 +40,6 @@
}
)
-
-## extract the estimated log likelihood
-setMethod('logLik','mif',function(object,...)object at loglik)
-
## extract the convergence record
conv.rec.internal <- function (object, pars, transform = FALSE, ...) {
if (transform) {
Added: pkg/pomp/R/mif2-methods.R
===================================================================
--- pkg/pomp/R/mif2-methods.R (rev 0)
+++ pkg/pomp/R/mif2-methods.R 2015-06-04 23:29:29 UTC (rev 1182)
@@ -0,0 +1,209 @@
+## ancillary methods for working with 'mif2d.pomp' objects
+
+setMethod('conv.rec','mif2d.pomp',
+ function (object, pars, transform = FALSE, ...) {
+ conv.rec.internal(object=object,pars=pars,transform=transform,...)
+ }
+ )
+
+## mif2List class
+setClass(
+ 'mif2List',
+ contains='list',
+ validity=function (object) {
+ if (!all(sapply(object,is,'mif2d.pomp'))) {
+ retval <- paste0(
+ "error in ",sQuote("c"),
+ ": dissimilar objects cannot be combined"
+ )
+ return(retval)
+ }
+ d <- sapply(object,function(x)dim(x at conv.rec))
+ if (!all(apply(d,1,diff)==0)) {
+ retval <- paste0(
+ "error in ",sQuote("c"),
+ ": to be combined, ",sQuote("mif2d.pomp"),
+ " objects must equal numbers of iterations"
+ )
+ return(retval)
+ }
+ TRUE
+ }
+ )
+
+setMethod(
+ 'c',
+ signature=signature(x='mif2d.pomp'),
+ definition=function (x, ...) {
+ y <- list(...)
+ if (length(y)==0) {
+ new("mif2List",list(x))
+ } else {
+ p <- sapply(y,is,'mif2d.pomp')
+ pl <- sapply(y,is,'mif2List')
+ if (any(!(p||pl)))
+ stop("cannot mix ",sQuote("mif2d.pomp"),
+ " and non-",sQuote("mif2d.pomp")," objects")
+ y[p] <- lapply(y[p],list)
+ y[pl] <- lapply(y[pl],as,"list")
+ new("mif2List",c(list(x),y,recursive=TRUE))
+ }
+ }
+ )
+
+setMethod(
+ 'c',
+ signature=signature(x='mif2List'),
+ definition=function (x, ...) {
+ y <- list(...)
+ if (length(y)==0) {
+ x
+ } else {
+ p <- sapply(y,is,'mif2d.pomp')
+ pl <- sapply(y,is,'mif2List')
+ if (any(!(p||pl)))
+ stop("cannot mix ",sQuote("mif2d.pomp"),
+ " and non-",sQuote("mif2d.pomp")," objects")
+ y[p] <- lapply(y[p],list)
+ y[pl] <- lapply(y[pl],as,"list")
+ new("mif2List",c(as(x,"list"),y,recursive=TRUE))
+ }
+ }
+ )
+
+setMethod(
+ "[",
+ signature=signature(x="mif2List"),
+ definition=function(x, i, ...) {
+ new('mif2List',as(x,"list")[i])
+ }
+ )
+
+setMethod(
+ 'conv.rec',
+ signature=signature(object='mif2List'),
+ definition=function (object, ...) {
+ lapply(object,conv.rec,...)
+ }
+ )
+
+mif2.diagnostics <- function (z) {
+ ## assumes that z is a list of mif2d.pomps with identical structure
+ mar.multi <- c(0,5.1,0,2.1)
+ oma.multi <- c(6,0,5,0)
+ xx <- z[[1]]
+ parnames <- names(coef(xx,transform=xx at transform))
+ estnames <- parnames
+
+ ## plot filter means
+ filt.diag <- rbind("eff. sample size"=xx at eff.sample.size,filter.mean(xx))
+ filtnames <- rownames(filt.diag)
+ plotnames <- filtnames
+ lognames <- filtnames[1] # eff. sample size
+ nplots <- length(plotnames)
+ n.per.page <- min(nplots,10)
+ if(n.per.page<=4) nc <- 1 else nc <- 2
+ nr <- ceiling(n.per.page/nc)
+ oldpar <- par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc),ask=dev.interactive(orNone=TRUE))
+ on.exit(par(oldpar))
+ low <- 1
+ hi <- 0
+ time <- time(xx)
+ while (hi<nplots) {
+ hi <- min(low+n.per.page-1,nplots)
+ for (i in seq(from=low,to=hi,by=1)) {
+ n <- i-low+1
+ logplot <- if (plotnames[i]%in%lognames) "y" else ""
+ dat <- sapply(
+ z,
+ function(po, label) {
+ if (label=="eff. sample size")
+ po at eff.sample.size
+ else
+ filter.mean(po,label)
+ },
+ label=plotnames[i]
+ )
+ matplot(
+ y=dat,
+ x=time,
+ axes = FALSE,
+ xlab = "",
+ log=logplot,
+ ylab = "",
+ type = "l"
+ )
+ box()
+ y.side <- 2
+ axis(y.side, xpd = NA)
+ mtext(plotnames[i], y.side, line = 3)
+ do.xax <- (n%%nr==0||n==n.per.page)
+ if (do.xax) axis(1,xpd=NA)
+ if (do.xax) mtext("time",side=1,line=3)
+ }
+ low <- hi+1
+ mtext("Filter diagnostics (last iteration)",3,line=2,outer=TRUE)
+ }
+
+ ## plot mif convergence diagnostics
+ other.diagnostics <- c("loglik", "nfail")
+ plotnames <- c(other.diagnostics,estnames)
+ nplots <- length(plotnames)
+ n.per.page <- min(nplots,10)
+ nc <- if (n.per.page<=4) 1 else 2
+ nr <- ceiling(n.per.page/nc)
+ par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc))
+ ## on.exit(par(oldpar))
+ low <- 1
+ hi <- 0
+ iteration <- seq(0,xx at Nmif)
+ while (hi<nplots) {
+ hi <- min(low+n.per.page-1,nplots)
+ for (i in seq(from=low,to=hi,by=1)) {
+ n <- i-low+1
+ dat <- sapply(z,function(po,label) conv.rec(po,label),label=plotnames[i])
+ matplot(
+ y=dat,
+ x=iteration,
+ axes = FALSE,
+ xlab = "",
+ ylab = "",
+ type = "l"
+ )
+ box()
+ y.side <- 2
+ axis(y.side,xpd=NA)
+ mtext(plotnames[i],y.side,line=3)
+ do.xax <- (n%%nr==0||n==n.per.page)
+ if (do.xax) axis(1,xpd=NA)
+ if (do.xax) mtext("MIF iteration",side=1,line=3)
+ }
+ low <- hi+1
+ mtext("MIF2 convergence diagnostics",3,line=2,outer=TRUE)
+ }
+ invisible(NULL)
+}
+
+setMethod(
+ "plot",
+ "mif2d.pomp",
+ function (x, y, ...) {
+ if (!missing(y)) {
+ y <- substitute(y)
+ warning(sQuote(y)," is ignored")
+ }
+ mif2.diagnostics(list(x))
+ }
+ )
+
+setMethod(
+ "plot",
+ signature=signature(x='mif2List'),
+ definition=function (x, y, ...) {
+ if (!missing(y)) {
+ y <- substitute(y)
+ warning(sQuote(y)," is ignored")
+ }
+ mif2.diagnostics(x)
+ }
+ )
Modified: pkg/pomp/R/mif2.R
===================================================================
--- pkg/pomp/R/mif2.R 2015-06-04 23:29:13 UTC (rev 1181)
+++ pkg/pomp/R/mif2.R 2015-06-04 23:29:29 UTC (rev 1182)
@@ -284,8 +284,6 @@
)
}
-setGeneric("mif2",function(object,...)standardGeneric("mif2"))
-
setMethod(
"mif2",
signature=signature(object="pomp"),
@@ -423,215 +421,3 @@
obj
}
)
-
-## extract the estimated log likelihood
-setMethod('logLik','mif2d.pomp',function(object,...)object at loglik)
-
-setMethod('conv.rec','mif2d.pomp',
- function (object, pars, transform = FALSE, ...) {
- conv.rec.internal(object=object,pars=pars,transform=transform,...)
- }
- )
-
-## mif2List class
-setClass(
- 'mif2List',
- contains='list',
- validity=function (object) {
- if (!all(sapply(object,is,'mif2d.pomp'))) {
- retval <- paste0(
- "error in ",sQuote("c"),
- ": dissimilar objects cannot be combined"
- )
- return(retval)
- }
- d <- sapply(object,function(x)dim(x at conv.rec))
- if (!all(apply(d,1,diff)==0)) {
- retval <- paste0(
- "error in ",sQuote("c"),
- ": to be combined, ",sQuote("mif2d.pomp"),
- " objects must equal numbers of iterations"
- )
- return(retval)
- }
- TRUE
- }
- )
-
-setMethod(
- 'c',
- signature=signature(x='mif2d.pomp'),
- definition=function (x, ...) {
- y <- list(...)
- if (length(y)==0) {
- new("mif2List",list(x))
- } else {
- p <- sapply(y,is,'mif2d.pomp')
- pl <- sapply(y,is,'mif2List')
- if (any(!(p||pl)))
- stop("cannot mix ",sQuote("mif2d.pomp"),
- " and non-",sQuote("mif2d.pomp")," objects")
- y[p] <- lapply(y[p],list)
- y[pl] <- lapply(y[pl],as,"list")
- new("mif2List",c(list(x),y,recursive=TRUE))
- }
- }
- )
-
-setMethod(
- 'c',
- signature=signature(x='mif2List'),
- definition=function (x, ...) {
- y <- list(...)
- if (length(y)==0) {
- x
- } else {
- p <- sapply(y,is,'mif2d.pomp')
- pl <- sapply(y,is,'mif2List')
- if (any(!(p||pl)))
- stop("cannot mix ",sQuote("mif2d.pomp"),
- " and non-",sQuote("mif2d.pomp")," objects")
- y[p] <- lapply(y[p],list)
- y[pl] <- lapply(y[pl],as,"list")
- new("mif2List",c(as(x,"list"),y,recursive=TRUE))
- }
- }
- )
-
-setMethod(
- "[",
- signature=signature(x="mif2List"),
- definition=function(x, i, ...) {
- new('mif2List',as(x,"list")[i])
- }
- )
-
-setMethod(
- 'conv.rec',
- signature=signature(object='mif2List'),
- definition=function (object, ...) {
- lapply(object,conv.rec,...)
- }
- )
-
-mif2.diagnostics <- function (z) {
- ## assumes that z is a list of mif2d.pomps with identical structure
- mar.multi <- c(0,5.1,0,2.1)
- oma.multi <- c(6,0,5,0)
- xx <- z[[1]]
- parnames <- names(coef(xx,transform=xx at transform))
- estnames <- parnames
-
- ## plot filter means
- filt.diag <- rbind("eff. sample size"=xx at eff.sample.size,filter.mean(xx))
- filtnames <- rownames(filt.diag)
- plotnames <- filtnames
- lognames <- filtnames[1] # eff. sample size
- nplots <- length(plotnames)
- n.per.page <- min(nplots,10)
- if(n.per.page<=4) nc <- 1 else nc <- 2
- nr <- ceiling(n.per.page/nc)
- oldpar <- par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc),ask=dev.interactive(orNone=TRUE))
- on.exit(par(oldpar))
- low <- 1
- hi <- 0
- time <- time(xx)
- while (hi<nplots) {
- hi <- min(low+n.per.page-1,nplots)
- for (i in seq(from=low,to=hi,by=1)) {
- n <- i-low+1
- logplot <- if (plotnames[i]%in%lognames) "y" else ""
- dat <- sapply(
- z,
- function(po, label) {
- if (label=="eff. sample size")
- po at eff.sample.size
- else
- filter.mean(po,label)
- },
- label=plotnames[i]
- )
- matplot(
- y=dat,
- x=time,
- axes = FALSE,
- xlab = "",
- log=logplot,
- ylab = "",
- type = "l"
- )
- box()
- y.side <- 2
- axis(y.side, xpd = NA)
- mtext(plotnames[i], y.side, line = 3)
- do.xax <- (n%%nr==0||n==n.per.page)
- if (do.xax) axis(1,xpd=NA)
- if (do.xax) mtext("time",side=1,line=3)
- }
- low <- hi+1
- mtext("Filter diagnostics (last iteration)",3,line=2,outer=TRUE)
- }
-
- ## plot mif convergence diagnostics
- other.diagnostics <- c("loglik", "nfail")
- plotnames <- c(other.diagnostics,estnames)
- nplots <- length(plotnames)
- n.per.page <- min(nplots,10)
- nc <- if (n.per.page<=4) 1 else 2
- nr <- ceiling(n.per.page/nc)
- par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc))
- ## on.exit(par(oldpar))
- low <- 1
- hi <- 0
- iteration <- seq(0,xx at Nmif)
- while (hi<nplots) {
- hi <- min(low+n.per.page-1,nplots)
- for (i in seq(from=low,to=hi,by=1)) {
- n <- i-low+1
- dat <- sapply(z,function(po,label) conv.rec(po,label),label=plotnames[i])
- matplot(
- y=dat,
- x=iteration,
- axes = FALSE,
- xlab = "",
- ylab = "",
- type = "l"
- )
- box()
- y.side <- 2
- axis(y.side,xpd=NA)
- mtext(plotnames[i],y.side,line=3)
- do.xax <- (n%%nr==0||n==n.per.page)
- if (do.xax) axis(1,xpd=NA)
- if (do.xax) mtext("MIF iteration",side=1,line=3)
- }
- low <- hi+1
- mtext("MIF convergence diagnostics",3,line=2,outer=TRUE)
- }
- invisible(NULL)
-}
-
-
-setMethod(
- "plot",
- "mif2d.pomp",
- function (x, y, ...) {
- if (!missing(y)) {
- y <- substitute(y)
- warning(sQuote(y)," is ignored")
- }
- mif2.diagnostics(list(x))
- }
- )
-
-setMethod(
- "plot",
- signature=signature(x='mif2List'),
- definition=function (x, y, ...) {
- if (!missing(y)) {
- y <- substitute(y)
- warning(sQuote(y)," is ignored")
- }
- mif2.diagnostics(x)
- }
- )
Modified: pkg/pomp/man/mif.Rd
===================================================================
--- pkg/pomp/man/mif.Rd 2015-06-04 23:29:13 UTC (rev 1181)
+++ pkg/pomp/man/mif.Rd 2015-06-04 23:29:29 UTC (rev 1182)
@@ -13,8 +13,6 @@
\alias{continue-mif}
\alias{mif-class}
\alias{mif-methods}
-\alias{logLik,mif-method}
-\alias{logLik-mif}
\alias{conv.rec}
\alias{conv.rec,mif-method}
\alias{conv.rec-mif}
@@ -31,10 +29,14 @@
\alias{c,mifList-method}
\alias{[-mifList}
\alias{[,mifList-method}
-\alias{compare.mif}
\description{
Iterated filtering algorithms for estimating the parameters of a partially-observed Markov process.
Running \code{mif} causes the iterated filtering algorithm to run for a specified number of iterations.
+ At each iteration, the particle filter is performed on a perturbed version of the model.
+ Specifically, parameters to be estimated are subjected to random perturbations at each observation.
+ This extra variability effectively smooths the likelihood surface and combats particle depletion by introducing diversity into the population of particles.
+ At the iterations progress, the magnitude of the perturbations is diminished according to a user-specified cooling schedule.
+ For most purposes, \code{mif} has been superseded by \code{\link{mif2}}.
}
\usage{
\S4method{mif}{pomp}(object, Nmif = 1, start, ivps = character(0),
@@ -49,7 +51,6 @@
cooling.type, cooling.fraction.50,
method, tol, transform, \dots)
\S4method{continue}{mif}(object, Nmif = 1, \dots)
-\S4method{logLik}{mif}(object, \dots)
\S4method{conv.rec}{mif}(object, pars, transform = FALSE, \dots)
\S4method{conv.rec}{mifList}(object, \dots)
}
@@ -141,9 +142,16 @@
}
\item{pars}{names of parameters.}
}
+\value{
+ Upon successful completion, \code{mif} returns an object of class \code{mif}.
+ The latter inherits from the \code{\link{pfilterd.pomp}} and \code{\link{pomp}} classes.
+}
+\section{IF2}{
+ A more full-featured version of the improved iterated filtering algorithm (IF2) is implemented as \code{\link{mif2}}.
+}
\section{Re-running \code{mif} Iterations}{
To re-run a sequence of \code{mif} iterations, one can use the \code{mif} method on a \code{mif} object.
- By default, the same parameters used for the original \code{mif} run are re-used (except for \code{weighted}, \code{tol}, \code{max.fail}, and \code{verbose}, the defaults of which are shown above).
+ By default, the same parameters used for the original \code{mif} run are re-used (except for \code{tol}, \code{max.fail}, and \code{verbose}, the defaults of which are shown above).
If one does specify additional arguments, these will override the defaults.
}
\section{Continuing \code{mif} Iterations}{
@@ -171,7 +179,7 @@
NB: this is \emph{not} the same as the likelihood of the model at the MLE!
}
\item{c}{
- Concatenates \code{mif} objects into an \code{mifList}.
+ Concatenates \code{mif} objects into a \code{mifList}.
}
\item{plot}{
Plots a series of diagnostic plots when applied to a \code{mif} or \code{mifList} object.
@@ -215,7 +223,7 @@
}
\author{Aaron A. King \email{kingaa at umich dot edu}}
\seealso{
- \code{\link{pomp}}, \code{\link{pfilter}}
+ \code{\link{pomp}}, \code{\link{pfilter}}, \code{\link{mif2}}
}
\keyword{optimize}
\keyword{ts}
Added: pkg/pomp/man/mif2.Rd
===================================================================
--- pkg/pomp/man/mif2.Rd (rev 0)
+++ pkg/pomp/man/mif2.Rd 2015-06-04 23:29:29 UTC (rev 1182)
@@ -0,0 +1,201 @@
+\name{Iterated filtering 2}
+\title{IF2: Maximum likelihood by iterated, perturbed Bayes maps}
+\docType{methods}
+\alias{mif2}
+\alias{mif2,mif2d.pomp-method}
+\alias{mif2-mif2d.pomp}
+\alias{mif2,pfilterd.pomp-method}
+\alias{mif2-pfilterd.pomp}
+\alias{mif2,pomp-method}
+\alias{mif2-pomp}
+\alias{continue,mif2d.pomp-method}
+\alias{continue-mif2d.pomp}
+\alias{mif2d.pomp-class}
+\alias{mif2d.pomp-methods}
+\alias{conv.rec,mif2d.pomp-method}
+\alias{conv.rec-mif2d.pomp}
+\alias{conv.rec,mif2List-method}
+\alias{conv.rec-mif2List}
+\alias{plot-mif2d.pomp}
+\alias{plot,mif2d.pomp-method}
+\alias{plot-mif2List}
+\alias{plot,mif2List-method}
+\alias{mif2List-class}
+\alias{c-mif2d.pomp}
+\alias{c,mif2d.pomp-method}
+\alias{c-mif2List}
+\alias{c,mif2List-method}
+\alias{[-mif2List}
+\alias{[,mif2List-method}
+\alias{rw.sd}
+\description{
+ An iterated filtering algorithm for estimating the parameters of a partially-observed Markov process.
+ Running \code{mif2} causes the improved iterated filtering algorithm (IF2) to run for a specified number of particle-filter iterations.
+ At each iteration, the particle filter is performed on a perturbed version of the model.
+ Specifically, parameters to be estimated are subjected to random perturbations at each observation.
+ This extra variability effectively smooths the likelihood surface and introduces diversity into the population of particles to combat depletion.
+ At the iterations progress, the magnitude of the perturbations is diminished according to a user-specified cooling schedule.
+ The algorithm is presented and justified in Ionides et al. (2015).
+}
+\usage{
+\S4method{mif2}{pomp}(object, Nmif = 1, start, Np, rw.sd, transform = FALSE,
+ cooling.type = c("hyperbolic", "geometric"),
+ cooling.fraction.50, tol = 1e-17, max.fail = Inf,
+ verbose = getOption("verbose"), \dots)
+\S4method{mif2}{pfilterd.pomp}(object, Nmif = 1, Np, tol, \dots)
+\S4method{mif2}{mif2d.pomp}(object, Nmif, start, Np, rw.sd,
+ transform, cooling.type, cooling.fraction.50, tol, \dots)
+\S4method{continue}{mif2d.pomp}(object, Nmif = 1, \dots)
+\S4method{conv.rec}{mif2d.pomp}(object, pars, transform = FALSE, \dots)
+\S4method{conv.rec}{mif2List}(object, \dots)
+rw.sd(\dots)
+}
+\arguments{
+ \item{object}{
+ An object of class \code{pomp}.
+ }
+ \item{Nmif}{
+ The number of filtering iterations to perform.
+ }
+ \item{start}{
+ named numerical vector;
+ the starting guess of the parameters.
+ By default, \code{start=coef(object)}.
+ }
+ \item{Np}{
+ the number of particles to use in filtering.
+ This may be specified as a single positive integer, in which case the same number of particles will be used at each timestep.
+ Alternatively, if one wishes the number of particles to vary across timestep, one may specify \code{Np} either as a vector of positive integers (of length \code{length(time(object,t0=TRUE))}) or as a function taking a positive integer argument.
+ In the latter case, \code{Np(k)} must be a single positive integer, representing the number of particles to be used at the \code{k}-th timestep:
+ \code{Np(0)} is the number of particles to use going from \code{timezero(object)} to \code{time(object)[1]},
+ \code{Np(1)}, from \code{timezero(object)} to \code{time(object)[1]},
+ and so on, while when \code{T=length(time(object,t0=TRUE))},
+ \code{Np(T)} is the number of particles to sample at the end of the time-series.
+ }
+ \item{rw.sd}{
+ specification of the magnitude of the random-walk perturbations that will be applied to some or all model parameters.
+ Parameters that are to be estimated should have positive perturbations specified here.
+ The specification is given using the \code{rw.sd} function (see below), which creates a list of unevaluated expressions.
+ The latter are evaluated in a context where the model time variable is defined (as \code{time}).
+ The expression \code{ivp(s)} can be used in this context as shorthand for \preformatted{ifelse(time==time[1],s,0).}
+ Likewise, \code{ivp(s,lag)} is equivalent to \preformatted{ifelse(time==time[lag],s,0).}
+ See below for some examples.
+ }
+ \item{transform}{
+ logical;
+ if \code{TRUE}, optimization is performed on the estimation scale, as defined by the user-supplied parameter transformations (see \code{\link{pomp}}).
+ }
+ \item{cooling.type, cooling.fraction.50}{
+ specifications for the cooling schedule, i.e., the manner in which the intensity of the parameter perturbations is reduced with successive filtering iterations.
+ \code{cooling.type} specifies the nature of the cooling schedule.
+
+ When \code{cooling.type="geometric"}, on the n-th \code{mif} iteration, the relative perturbation intensity is \code{cooling.fraction.50^(n/50)}.
+
+ When \code{cooling.type="hyperbolic"}, on the n-th \code{mif} iteration, the relative perturbation intensity is \code{(s+1)/(s+n)}, where \code{(s+1)/(s+50)=cooling.fraction.50}.
+ \code{cooling.fraction.50} is the relative magnitude of the parameter perturbations after 50 \code{mif} iterations.
+ }
+ \item{tol, max.fail}{
+ passed to the particle filter.
+ See the descriptions under \code{\link{pfilter}}.
+ }
+ \item{verbose}{
+ logical; if TRUE, print progress reports.
+ }
+ \item{\dots}{
+ additional arguments that override the defaults.
+ }
+ \item{pars}{names of parameters.}
+}
+\value{
+ Upon successful completion, \code{mif2} returns an object of class \code{mif2d.pomp}.
+ This class inherits from the \code{\link{pfilterd.pomp}} and \code{\link{pomp}} classes.
+}
+\section{The \code{rw.sd} function}{
+ This function simply returns a list containing its arguments as unevaluated expressions.
+ These are then evaluated in a context containing the model \code{time} variable.
+ This allows for easy specification of the structure of the perturbations that are to be applied.
+}
+\section{Re-running \code{mif2} Iterations}{
+ To re-run a sequence of \code{mif2} iterations, one can use the \code{mif2} method on a \code{mif2d.pomp} object.
+ By default, the same parameters used for the original \code{mif2} run are re-used (except for \code{tol}, \code{max.fail}, and \code{verbose}, the defaults of which are shown above).
+ If one does specify additional arguments, these will override the defaults.
+}
+\section{Continuing \code{mif2} Iterations}{
+ One can resume a series of \code{mif2} iterations from where one left off using the \code{continue} method.
+ A call to \code{mif2} to perform \code{Nmif=m} iterations followed by a call to \code{continue} to perform \code{Nmif=n} iterations will produce precisely the same effect as a single call to \code{mif2} to perform \code{Nmif=m+n} iterations.
+ By default, all the algorithmic parameters are the same as used in the original call to \code{mif2}.
+ Additional arguments will override these defaults.
+}
+\section{Methods}{
+ Methods that can be used to manipulate, display, or extract information from a \code{mif2d.pomp} object:
+ \describe{
+ \item{conv.rec}{
+ \code{conv.rec(object, pars = NULL)} returns the columns of the convergence-record matrix corresponding to the names in \code{pars}.
+ By default, all rows are returned.
+ }
+ \item{logLik}{
+ Returns the value in the \code{loglik} slot.
+ NB: this is \emph{not} the same as the likelihood of the model at the MLE!
+ }
+ \item{c}{
+ Concatenates \code{mif2d.pomp} objects into a \code{mif2List}.
+ }
+ \item{plot}{
+ Plots a series of diagnostic plots when applied to a \code{mif2d.pomp} or \code{mif2List} object.
+ }
+ }
+}
+\section{Details}{
+ If \code{particles} is not specified, the default behavior is to draw the particles from a multivariate normal distribution.
+ \strong{It is the user's responsibility to ensure that, if the optional \code{particles} argument is given, that the \code{particles} function satisfies the following conditions:}
+
+ \code{particles} has at least the following arguments:
+ \code{Np}, \code{center}, \code{sd}, and \code{\dots}.
+ \code{Np} may be assumed to be a positive integer;
+ \code{center} and \code{sd} will be named vectors of the same length.
+ Additional arguments may be specified;
+ these will be filled with the elements of the \code{userdata} slot of the underlying \code{pomp} object (see \code{\link{pomp}}).
+
+ \code{particles} returns a \code{length(center)} x \code{Np} matrix with rownames matching the names of \code{center} and \code{sd}.
+ Each column represents a distinct particle.
+
+ The center of the particle distribution returned by \code{particles} should be \code{center}.
+ The width of the particle distribution should vary monotonically with \code{sd}.
+ In particular, when \code{sd=0}, the \code{particles} should return matrices with \code{Np} identical columns, each given by the parameters specified in \code{center}.
+}
+\examples{
+\dontrun{
+pompExample(ou2)
+
+guess1 <- guess2 <- coef(ou2)
+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 <- 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 <- 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))
+
+rbind(mle1=c(coef(m1),loglik=logLik(pfilter(m1,Np=1000))),
+ mle2=c(coef(m2),loglik=logLik(pfilter(m1,Np=1000))),
+ truth=c(coef(ou2),loglik=logLik(pfilter(m1,Np=1000))))
+}
+}
+\references{
+ E. L. Ionides, D. Nguyen, Y. Atchad{\\'e}, S. Stoev, and A. A. King.
+ Inference for dynamic and latent variable models via iterated, perturbed Bayes maps.
+ Proc. Natl. Acad. Sci. U.S.A., 112:719--724, 2015.
+}
+\author{Aaron A. King \email{kingaa at umich dot edu}, Edward L. Ionides, and Dao Nguyen}
+\seealso{
+ \code{\link{pomp}}, \code{\link{pfilter}}, \code{\link{mif}}, and the tutorials on the \href{http://pomp.r-forge.r-project.org}{package website}.
+}
+\keyword{optimize}
+\keyword{ts}
Modified: pkg/pomp/man/pfilter.Rd
===================================================================
--- pkg/pomp/man/pfilter.Rd 2015-06-04 23:29:13 UTC (rev 1181)
+++ pkg/pomp/man/pfilter.Rd 2015-06-04 23:29:29 UTC (rev 1182)
@@ -9,6 +9,7 @@
\alias{pfilter,pfilterd.pomp-method}
\alias{pfilter-pfilterd.pomp}
\alias{pfilterd.pomp-class}
+\alias{pfilterd.pomp}
\alias{logLik,pfilterd.pomp-method}
\alias{logLik-pfilterd.pomp}
\alias{$,pfilterd.pomp-method}
Modified: pkg/pomp/man/pmcmc.Rd
===================================================================
--- pkg/pomp/man/pmcmc.Rd 2015-06-04 23:29:13 UTC (rev 1181)
+++ pkg/pomp/man/pmcmc.Rd 2015-06-04 23:29:29 UTC (rev 1182)
@@ -32,7 +32,6 @@
\description{
The Particle MCMC algorithm for estimating the parameters of a partially-observed Markov process.
Running \code{pmcmc} causes a particle random-walk Metropolis-Hastings Markov chain algorithm to run for the specified number of proposals.
-
}
\usage{
\S4method{pmcmc}{pomp}(object, Nmcmc = 1, start, proposal, pars, rw.sd, Np,
Modified: pkg/pomp/tests/ou2-mif2.R
===================================================================
--- pkg/pomp/tests/ou2-mif2.R 2015-06-04 23:29:13 UTC (rev 1181)
+++ pkg/pomp/tests/ou2-mif2.R 2015-06-04 23:29:29 UTC (rev 1182)
@@ -113,11 +113,12 @@
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,
+m2 <- mif2(ou2,Nmif=50,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))
+m2 <- continue(m2,Nmif=50)
plot(c(m1,m2))
Modified: pkg/pomp/tests/ou2-mif2.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-mif2.Rout.save 2015-06-04 23:29:13 UTC (rev 1181)
+++ pkg/pomp/tests/ou2-mif2.Rout.save 2015-06-04 23:29:29 UTC (rev 1182)
@@ -144,11 +144,12 @@
+ 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,
+> m2 <- mif2(ou2,Nmif=50,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))
+> m2 <- continue(m2,Nmif=50)
>
> plot(c(m1,m2))
>
@@ -181,4 +182,4 @@
>
> proc.time()
user system elapsed
- 69.328 0.052 69.748
+ 68.756 0.088 69.231
More information about the pomp-commits
mailing list