Learn R Programming

pomp (version 1.19)

Particle filter: Particle filter

Description

A plain vanilla sequential Monte Carlo (particle filter) algorithm. Resampling is performed at each observation.

Usage

# S4 method for pomp
pfilter(object, params, Np, tol = 1e-17,
    max.fail = Inf, pred.mean = FALSE, pred.var = FALSE,
    filter.mean = FALSE, filter.traj = FALSE, save.states = FALSE,
    save.params = FALSE, verbose = getOption("verbose"), …)
# S4 method for pfilterd.pomp
pfilter(object, params, Np, tol, …)
# S4 method for pfilterd.pomp
logLik(object, …)
# S4 method for pfilterd.pomp
cond.logLik(object, …)
# S4 method for pfilterd.pomp
eff.sample.size(object, …)
# S4 method for pfilterd.pomp
pred.mean(object, pars, …)
# S4 method for pfilterd.pomp
pred.var(object, pars, …)
# S4 method for pfilterd.pomp
filter.mean(object, pars, …)
# S4 method for pfilterd.pomp
filter.traj(object, vars, …)

Arguments

object

An object of class pomp or inheriting class pomp.

params

optional named numeric vector containing the parameters at which the filtering should be performed. By default, params = coef(object).

Np

the number of particles to use. 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 timesteps, 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. When object is of class mif, this is by default the same number of particles used in the mif iterations.

One should omit Np if params is a matrix of parameters, with one column for each particle. In this case, obviously, the number of particles is ncol(params).

tol

positive numeric scalar; particles with likelihood less than tol are considered to be incompatible with the data. See the section on Filtering Failures below for more information.

max.fail

integer; the maximum number of filtering failures allowed (see below). If the number of filtering failures exceeds this number, execution will terminate with an error. By default, max.fail is set to infinity, so no error can be triggered.

pred.mean

logical; if TRUE, the prediction means are calculated for the state variables and parameters.

pred.var

logical; if TRUE, the prediction variances are calculated for the state variables and parameters.

filter.mean

logical; if TRUE, the filtering means are calculated for the state variables and parameters.

filter.traj

logical; if TRUE, a filtered trajectory is returned for the state variables and parameters.

save.states, save.params

logical. If save.states=TRUE, the state-vector for each particle at each time is saved in the saved.states slot of the returned pfilterd.pomp object. If save.params=TRUE, the parameter-vector for each particle at each time is saved in the saved.params slot of the returned pfilterd.pomp object.

verbose

logical; if TRUE, progress information is reported as pfilter works.

pars

Names of parameters.

vars

Names of state variables.

additional arguments that override the defaults.

Value

An object of class pfilterd.pomp. This class inherits from class pomp. The following additional slots can be accessed via the $ operator:

saved.states

If pfilter was called with save.states=TRUE, this is the list of state-vectors at each time point, for each particle. It is a length-ntimes list of nvars-by-Np arrays. In particular, saved.states[[t]][,i] can be considered a sample from \(f[X_t|y_{1:t}]\).

saved.params

If pfilter was called with save.params=TRUE, this is the list of parameter-vectors at each time point, for each particle. It is a length-ntimes list of npars-by-Np arrays. In particular, saved.params[[t]][,i] is the parameter portion of the i-th particle at time \(t\).

Np, tol, nfail

the number of particles used, failure tolerance, and number of filtering failures (see below), respectively.

Methods

logLik

Extracts the estimated log likelihood.

cond.logLik

Extracts the estimated conditional log likelihood $$\ell_t(\theta) = \mathrm{Prob}[y_t \vert y_1, \dots, y_{t-1}],$$ where \(y_t\) are the data, at time \(t\).

eff.sample.size

Extracts the (time-dependent) estimated effective sample size, computed as $$\left(\sum_i\!w_{it}^2\right)^{-1},$$ where \(w_{it}\) is the normalized weight of particle \(i\) at time \(t\).

pred.mean, pred.var

Extract the mean and variance of the approximate prediction distribution. This prediction distribution is that of $$X_t \vert y_1,\dots,y_{t-1},$$ where \(X_t\), \(y_t\) are the state vector and data, respectively, at time \(t\).

filter.mean

Extract the mean of the filtering distribution, which is that of $$X_t \vert y_1,\dots,y_t,$$ where \(X_t\), \(y_t\) are the state vector and data, respectively, at time \(t\).

Filtering failures

If the degree of disagreement between model and data becomes sufficiently large, a “filtering failure” results. A filtering failure occurs when, at some time point, none of the Np particles is compatible with the data. In particular, if the conditional likelihood of a particle at any time is below the tolerance value tol, then that particle is considered to be uninformative and its likelihood is taken to be zero. A filtering failure occurs when this is the case for all particles. A warning is generated when this occurs unless the cumulative number of failures exceeds max.fail, in which case an error is generated.

References

M. S. Arulampalam, S. Maskell, N. Gordon, & T. Clapp. A Tutorial on Particle Filters for Online Nonlinear, Non-Gaussian Bayesian Tracking. IEEE Trans. Sig. Proc. 50:174--188, 2002.

See Also

pomp, mif, pmcmc, bsmc2, and the tutorials on the package website.

Examples

Run this code
# NOT RUN {
pompExample(gompertz)
pf <- pfilter(gompertz,Np=1000)	## use 1000 particles
plot(pf)
logLik(pf)
cond.logLik(pf)			## conditional log-likelihoods
eff.sample.size(pf)             ## effective sample size
logLik(pfilter(pf))      	## run it again with 1000 particles
## run it again with 2000 particles
pf <- pfilter(pf,Np=2000,filter.mean=TRUE)
fm <- filter.mean(pf)    	## extract the filtering means
# }

Run the code above in your browser using DataLab