A plain vanilla sequential Monte Carlo (particle filter) algorithm. Resampling is performed at each observation.
# 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, …)An object of class pomp or inheriting class pomp.
optional named numeric vector containing the parameters at which the filtering should be performed.
By default, params = coef(object).
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).
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.
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.
logical; if TRUE, the prediction means are calculated for the state variables and parameters.
logical; if TRUE, the prediction variances are calculated for the state variables and parameters.
logical; if TRUE, the filtering means are calculated for the state variables and parameters.
logical; if TRUE, a filtered trajectory is returned for the state variables and parameters.
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.
logical; if TRUE, progress information is reported as pfilter works.
Names of parameters.
Names of state variables.
additional arguments that override the defaults.
An object of class pfilterd.pomp.
This class inherits from class pomp.
The following additional slots can be accessed via the $ operator:
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}]\).
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\).
the number of particles used, failure tolerance, and number of filtering failures (see below), respectively.
Extracts the estimated log likelihood.
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\).
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\).
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\).
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\).
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.
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.
pomp, mif, pmcmc, bsmc2,
and the tutorials on the package website.
# 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