nimbleSMC (version 0.11.0)

buildBootstrapFilter: Create a bootstrap particle filter algorithm to estimate log-likelihood.

Description

Create a bootstrap particle filter algorithm for a given NIMBLE state space model.

Usage

buildBootstrapFilter(model, nodes, control = list())

Arguments

model

A nimble model object, typically representing a state space model or a hidden Markov model.

nodes

A character vector specifying the stochastic latent model nodes over which the particle filter will stochastically integrate to estimate the log-likelihood function. All provided nodes must be stochastic. Can be one of three forms: a variable name, in which case all elements in the variable are taken to be latent (e.g., 'x'); an indexed variable, in which case all indexed elements are taken to be latent (e.g., 'x[1:100]' or 'x[1:100, 1:2]'); or a vector of multiple nodes, one per time point, in increasing time order (e.g., c("x[1:2, 1]", "x[1:2, 2]", "x[1:2, 3]", "x[1:2, 4]")).

control

A list specifying different control options for the particle filter. Options are described in the details section below.

<code>returnESS()</code> Method

Calling the returnESS() method of a bootstrap filter after that filter has been run() for a given model will return a vector of ESS (effective sample size) values, one value for each time point.

Author

Daniel Turek and Nicholas Michaud

Details

Each of the control() list options are described in detail here:

thresh

A number between 0 and 1 specifying when to resample: the resampling step will occur when the effective sample size is less than thresh times the number of particles. Defaults to 0.8. Note that at the last time step, resampling will always occur so that the mvEWsamples modelValues contains equally-weighted samples.

resamplingMethod

The type of resampling algorithm to be used within the particle filter. Can choose between 'default' (which uses NIMBLE's rankSample() function), 'systematic', 'stratified', 'residual', and 'multinomial'. Defaults to 'default'. Resampling methods other than 'default' are currently experimental.

saveAll

Indicates whether to save state samples for all time points (TRUE), or only for the most recent time point (FALSE)

smoothing

Decides whether to save smoothed estimates of latent states, i.e., samples from f(x[1:t]|y[1:t]) if smoothing = TRUE, or instead to save filtered samples from f(x[t]|y[1:t]) if smoothing = FALSE. smoothing = TRUE only works if saveAll = TRUE.

timeIndex

An integer used to manually specify which dimension of the latent state variable indexes time. Only needs to be set if the number of time points is less than or equal to the size of the latent state at each time point.

initModel

A logical value indicating whether to initialize the model before running the filtering algorithm. Defaults to TRUE.

The bootstrap filter starts by generating a sample of estimates from the prior distribution of the latent states of a state space model. At each time point, these particles are propagated forward by the model's transition equation. Each particle is then given a weight proportional to the value of the observation equation given that particle. The weights are used to draw an equally-weighted sample of the particles at this time point. The algorithm then proceeds to the next time point. Neither the transition nor the observation equations are required to be normal for the bootstrap filter to work.

The resulting specialized particle filter algorthm will accept a single integer argument (m, default 10,000), which specifies the number of random \'particles\' to use for estimating the log-likelihood. The algorithm returns the estimated log-likelihood value, and saves unequally weighted samples from the posterior distribution of the latent states in the mvWSamples modelValues object, with corresponding logged weights in mvWSamples['wts',]. An equally weighted sample from the posterior can be found in the mvEWSamples modelValues object.

Note that if the thresh argument is set to a value less than 1, resampling may not take place at every time point. At time points for which resampling did not take place, mvEWSamples will not contain equally weighted samples. To ensure equally weighted samples in the case that thresh < 1, we recommend resampling from mvWSamples at each time point after the filter has been run, rather than using mvEWSamples.

References

Gordon, N.J., D.J. Salmond, and A.F.M. Smith. (1993). Novel approach to nonlinear/non-Gaussian Bayesian state estimation. IEEE Proceedings F (Radar and Signal Processing). Vol. 140. No. 2. IET Digital Library, 1993.

See Also

Other particle filtering methods: buildAuxiliaryFilter, buildEnsembleKF, buildIteratedFilter2(), buildLiuWestFilter

Examples

## For illustration only.
exampleCode <- nimbleCode({
  x0 ~ dnorm(0, var = 1)
  x[1] ~ dnorm(.8 * x0, var = 1)
  y[1] ~ dnorm(x[1], var = .5)
  for(t in 2:10){
    x[t] ~ dnorm(.8 * x[t-1], var = 1)
    y[t] ~ dnorm(x[t], var = .5)
  }
})

model <- nimbleModel(code = exampleCode, data = list(y = rnorm(10)),
                     inits = list(x0 = 0, x = rnorm(10)))
my_BootF <- buildBootstrapFilter(model, 'x',
                control = list(saveAll = TRUE, thresh = 1))
## Now compile and run, e.g.,
## Cmodel <- compileNimble(model)
## Cmy_BootF <- compileNimble(my_BootF, project = model)
## logLik <- Cmy_BootF$run(m = 1000)
## ESS <- Cmy_BootF$returnESS()
## boot_X <- as.matrix(Cmy_BootF$mvEWSamples, 'x')