mif: Maximum likelihood by iterated filtering

In pomp: Statistical Inference for Partially Observed Markov Processes

Description

Iterated filtering algorithms for estimating the parameters of a partially-observed Markov process. Running 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, mif has been superseded by mif2.

Usage

## S4 method for signature 'pomp'
mif(object, Nmif = 1, start, ivps = character(0),
    particles, rw.sd, Np, ic.lag, var.factor = 1,
    cooling.type, cooling.fraction.50,
    method = c("mif","unweighted","fp","mif2"),
    tol = 1e-17, max.fail = Inf,
    verbose = getOption("verbose"), transform = FALSE, ...)
## S4 method for signature 'pfilterd.pomp'
mif(object, Nmif = 1, Np, tol, ...)
## S4 method for signature 'mif'
mif(object, Nmif, start, ivps,
    particles, rw.sd, Np, ic.lag, var.factor,
    cooling.type, cooling.fraction.50,
    method, tol, transform, ...)
## S4 method for signature 'mif'
continue(object, Nmif = 1, ...)
## S4 method for signature 'mif'
conv.rec(object, pars, transform = FALSE, ...)
## S4 method for signature 'mifList'
conv.rec(object, ...)

Arguments

object	An object of class pomp.
Nmif	The number of filtering iterations to perform.
start	named numerical vector; the starting guess of the parameters.
ivps	optional character vector naming the initial-value parameters (IVPs) to be estimated. Every parameter named in ivps must have a positive random-walk standard deviation specified in rw.sd. If there are no non-IVP parameters with positive rw.sd, i.e., only IVPs are to be estimated, see below “"Using mif to estimate initial-value parameters only"”.
particles	Function of prototype particles(Np,center,sd,...) which sets up the starting particle matrix by drawing a sample of size Np from the starting particle distribution centered at center and of width sd. If particles is not supplied by the user, the default behavior is to draw the particles from a multivariate normal distribution with mean center and standard deviation sd.
rw.sd	numeric vector with names; the intensity of the random walk to be applied to parameters. The random walk is only applied to parameters named in pars (i.e., not to those named in ivps). The algorithm requires that the random walk be nontrivial, so each element in rw.sd[pars] must be positive. rw.sd is also used to scale the initial-value parameters (via the particles function). Therefore, each element of rw.sd[ivps] must be positive. The following must be satisfied: names(rw.sd) must be a subset of names(start), rw.sd must be non-negative (zeros are simply ignored), the name of every positive element of rw.sd must be in either pars or ivps.
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 Np either as a vector of positive integers (of length length(time(object,t0=TRUE))) or as a function taking a positive integer argument. In the latter case, Np(k) must be a single positive integer, representing the number of particles to be used at the k-th timestep: Np(0) is the number of particles to use going from timezero(object) to time(object)[1], Np(1), from timezero(object) to time(object)[1], and so on, while when T=length(time(object,t0=TRUE)), Np(T) is the number of particles to sample at the end of the time-series.
ic.lag	a positive integer; the timepoint for fixed-lag smoothing of initial-value parameters. The mif update for initial-value parameters consists of replacing them by their filtering mean at time times[ic.lag], where times=time(object). It makes no sense to set ic.lag>length(times); if it is so set, ic.lag is set to length(times) with a warning. For method="mif2", the default is ic.lag=length(times).
var.factor	optional positive scalar; the scaling coefficient relating the width of the starting particle distribution to rw.sd. In particular, the width of the distribution of particles at the start of the first mif iteration will be random.walk.sd*var.factor. By default, var.factor=1.
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. cooling.type specifies the nature of the cooling schedule. When cooling.type="geometric", on the n-th mif iteration, the relative perturbation intensity is cooling.fraction.50^(n/50). When cooling.type="hyperbolic", on the n-th mif iteration, the relative perturbation intensity is (s+1)/(s+n), where (s+1)/(s+50)=cooling.fraction.50. cooling.fraction.50 is the relative magnitude of the parameter perturbations after 50 mif iterations.
method	method sets the update rule used in the algorithm. method="mif" uses the iterated filtering update rule (Ionides 2006, 2011); method="unweighted" updates the parameter to the unweighted average of the filtering means of the parameters at each time; method="fp" updates the parameter to the filtering mean at the end of the time series. method="mif2" implements an incomplete version of the iterated Bayes map method of Ionides (2015). The latter method is, by every indication, both more efficient and more stable.
tol, max.fail	See the description under pfilter.
verbose	logical; if TRUE, print progress reports.
transform	logical; if TRUE, optimization is performed on the transformed scale, as defined by the user-supplied parameter transformations (see pomp).
...	additional arguments that override the defaults.
pars	names of parameters.

Value

Upon successful completion, mif returns an object of class mif. The latter inherits from the pfilterd.pomp and pomp classes.

IF2

A more full-featured version of the improved iterated filtering algorithm (IF2) is implemented as mif2.

Re-running mif Iterations

To re-run a sequence of mif iterations, one can use the mif method on a mif object. By default, the same parameters used for the original mif run are re-used (except for tol, max.fail, and verbose, the defaults of which are shown above). If one does specify additional arguments, these will override the defaults.

Continuing mif Iterations

One can resume a series of mif iterations from where one left off using the continue method. A call to mif to perform Nmif=m iterations followed by a call to continue to perform Nmif=n iterations will produce precisely the same effect as a single call to mif to perform Nmif=m+n iterations. By default, all the algorithmic parameters are the same as used in the original call to mif. Additional arguments will override the defaults.

Using mif to estimate initial-value parameters only

One can use mif's fixed-lag smoothing to estimate only initial value parameters (IVPs). In this case, pars is left empty and the IVPs to be estimated are named in ivps. If theta is the current parameter vector, then at each mif iteration, Np particles are drawn from a distribution centered at theta and with width proportional to var.factor*rw.sd, a particle filtering operation is performed, and theta is replaced by the filtering mean at time(object)[ic.lag]. Note the implication that, when mif is used in this way on a time series any longer than ic.lag, unnecessary work is done. If the time series in object is longer than ic.lag, consider replacing object with window(object,end=ic.lag).

Methods

Methods that can be used to manipulate, display, or extract information from a mif object:

conv.rec

conv.rec(object, pars = NULL) returns the columns of the convergence-record matrix corresponding to the names in pars. By default, all rows are returned.

logLik

Returns the value in the loglik slot. NB: this is not the same as the likelihood of the model at the MLE!

c

Concatenates mif objects into a mifList.

plot

Plots a series of diagnostic plots when applied to a mif or mifList object.

Details

If particles is not specified, the default behavior is to draw the particles from a multivariate normal distribution. It is the user's responsibility to ensure that, if the optional particles argument is given, that the particles function satisfies the following conditions:

particles has at least the following arguments: Np, center, sd, and .... Np may be assumed to be a positive integer; center and sd will be named vectors of the same length. Additional arguments may be specified; these will be filled with the elements of the userdata slot of the underlying pomp object (see pomp).

particles returns a length(center) x Np matrix with rownames matching the names of center and sd. Each column represents a distinct particle.

The center of the particle distribution returned by particles should be center. The width of the particle distribution should vary monotonically with sd. In particular, when sd=0, the particles should return matrices with Np identical columns, each given by the parameters specified in center.

Author(s)

Aaron A. King kingaa at umich dot edu

References

E. L. Ionides, C. Breto, & A. A. King, Inference for nonlinear dynamical systems, Proc. Natl. Acad. Sci. U.S.A., 103:18438–18443, 2006.

E. L. Ionides, A. Bhadra, Y. Atchad\'e, & A. A. King, Iterated filtering, Annals of Statistics, 39:1776–1802, 2011.

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.

A. A. King, E. L. Ionides, M. Pascual, and M. J. Bouma, Inapparent infections and cholera dynamics, Nature, 454:877–880, 2008.

See Also

pomp, pfilter, mif2

pomp documentation built on May 2, 2019, 4:09 p.m.