Using GenericSSMs.jl

This section discusses how to use GenericSSMs.jl for SSMs defined as discussed in Section Defining SSMs using GenericSSMs.jl.

Most of the methods exported by GenericSSMs.jl do not (heap) allocate memory, but instead operate on storage objects, such that a typical method call (in this case to a method named operation) takes the form:

operation!(storage, ... additional arguments ...)

The object storage is simply an object that collects all pre-allocated memory necessary for making a particular computation (such as particle filtering or predicting). There are also `allocating versions` (lacking the exclamation mark) of some of the functionality, but the non-allocating versions should be used for performance especially in situations where a particular method needs to called multiple times (see Julia performance tips: pre-allocating outputs).

There are three storage objects used in GenericSSMs.jl:

PFStorage (for particle filtering)
CPFStorage (for conditional particle filtering)
PredictStorage (for prediction)

These are discussed below with their respective use cases. In what follows, we will in general use N to refer to the number of particles, and n to the number of time points (or length of the observed time series).

Particle filtering

The particle filter implementation of GenericSSMs.jl can be invoked by calling pf_forward_pass!:

GenericSSMs.pf_forward_pass! — Method

pf_forward_pass!(st::PFStorage, model::GenericSSM, resampling[, rng = Random.GLOBAL_RNG])

Run the standard particle filter for model model using storage st and resampling resampling that implements resample! (see ?GenericSSMs.resample!). An RNG may be specified as the last argument.

The return value is the logarithm of the normalising constant estimate of the particle filter.

pf_forward_pass! populates its first argument, a PFStorage object, which has the following structure:

GenericSSMs.PFStorage — Type

A storage object used for standard particle filtering. The object contains all memory needed for particle filtering.

Type parameters:

P: Particle type.
F: Floating point type.

Fields:

X::Matrix{P}, the N x n particle system of particles
W::Matrix{F}, the N x n unnormalised logweight matrix,
A::Matrix{Int}, the N x n - 1 ancestor index matrix.

The additional field _w is used internally for normalised weights (do not modify). The ancestor of the particle X[i, k] is X[A[i, k - 1], k - 1].

The easiest way to construct a PFStorage object is via the constructor:

GenericSSMs.PFStorage — Method

PFStorage(model::GenericSSM, N::Integer, n::Integer; F::Type{<: AbstractFloat} = Float64)

Construct a storage object for standard particle filtering model model with N particles and time series length n. F may be used to set the floating point type.

ASSUMPTIONS:

N >= 2, n >= 1 (checked, throws error)

The other arguments to pf_forward_pass! above are

the SSM model (expected to be defined as in Defining SSMs using GenericSSMs.jl)
resampling object, that implements resample! from the resampling API of GenericSSMs.jl

For the latter, it is possible to use any resampling provided by Resamplings.jl (see Section Examples for examples of this).

Finally, for convenience, pf_forward_pass! can also be called with a CPFStorage storage object (discussed in Conditional particle filtering), as follows:

GenericSSMs.pf_forward_pass! — Method

pf_forward_pass!(st::CPFStorage, model::GenericSSM, resampling[, rng = Random.GLOBAL_RNG])

Additional method for pf_forward_pass! that allows running the standard particle filter with the CPFStorage storage object. This method does not alter the field st.ref since it is not needed in standard particle filtering.

Conditional particle filtering

The conditional particle filter (CPF) is otherwise similar to the (standard) particle filter, but features a `reference trajectory` that is never mutated as the algorithm proceeds. Due to this, the CPF requires its own storage object, CPFStorage:

GenericSSMs.CPFStorage — Type

A storage object used for conditional particle filtering. The object contains all memory needed for conditional particle filtering.

Type parameters:

P: Particle type.
F: Floating point type.

Fields:

pfst::PFStorage{P, F}, the storage used for standard particle filtering (see ?PFStorage),
ref::Vector{Int}, a vector of length n that records which rows (elements) in each column of pfst.X currently hold the reference trajectory.

The construction of a CPFStorage object is similar to that of PFStorage:

GenericSSMs.CPFStorage — Method

CPFStorage(model::GenericSSM, N::Integer, n::Integer; F::Type{<: AbstractFloat} = Float64)

Construct a storage object for conditional particle filtering model model with N particles and time series length n. F may be used to set the floating point type.

ASSUMPTIONS:

N >= 2, n >= 1 (checked, throws error)

The main function for running the (forward) conditional particle filter is called cpf_forward_pass!:

GenericSSMs.cpf_forward_pass! — Function

cpf_forward_pass!(st::CPFStorage, model::GenericSSM, resampling[, rng = Random.GLOBAL_RNG])

Run the forward pass of the conditional particle filter using storage st for model model with resampling resampling that implements conditional_resample! (see ?GenericSSMs.conditional_resample!). An RNG may be specified as the last argument.

ASSUMPTIONS:

the reference trajectory has been initialised to st (fields st.ref and st.pfst have been populated). Results output by this function WILL BE WRONG, if this has not been done. Use initialise_reference! to properly initialise reference.

Here, the resampling object must implement conditional_resample! from the resampling API of GenericSSMs.jl. Again, any conditional resampling from Resamplings.jl may be used.

Note that care must be taken to ensure that the reference trajectory has been initialised to a sensible value before calling cpf_forward_pass!. For the first run of cpf_forward_pass!, a straightforward way to do this is to call initialise_reference! which uses the standard particle filter to initialise the reference:

GenericSSMs.initialise_reference! — Function

initialise_reference!(st::CPFStorage, model::GenericSSM, resampling[, rng = Random.GLOBAL_RNG])

Initialise the reference trajectory for running the conditional particle filter forward pass (cpf_forward_pass!).

To do the CPF backward pass, the function traceback! can be used, with either AncestorTracing or BackwardSampling. This effectively records a new reference trajectory to CPFStorage.

GenericSSMs.traceback! — Method

traceback!(st::CPFStorage, model::GenericSSM, AncestorTracing[, rng = Random.GLOBAL_RNG])

Populate the reference indices st.ref using ancestor tracing. To obtain the concrete reference trajectory values, use get_reference!.

ASSUMPTIONS:

the forward pass (pf_forward_pass! or cpf_forward_pass!) has been run before invoking this function and the contents of st have not changed after that. The results from this function WILL BE WRONG, if this is not the case.

GenericSSMs.traceback! — Method

traceback!(st::CPFStorage, AncestorTracing[, rng = Random.GLOBAL_RNG])

Convenience method for traceback using ancestor tracing (model not needed by ancestor tracing).

ASSUMPTIONS:

the forward pass (pf_forward_pass! or cpf_forward_pass!) has been run before invoking this function and the contents of st have not changed after that. The results from this function WILL BE WRONG, if this is not the case.

GenericSSMs.traceback! — Method

traceback!(st::CPFStorage, model::GenericSSM, BackwardSampling[, rng = Random.GLOBAL_RNG])

Populate reference indices st.ref using backward sampling after running the forward pass, that is, pf_forward_pass! or cpf_forward_pass!. To obtain the concrete reference trajectory values, see get_reference!.

ASSUMPTIONS:

the forward pass (pf_forward_pass! or cpf_forward_pass!) has been run before invoking this function and the contents of st have not changed after that. The results from this function WILL BE WRONG, if this is not the case.

The function get_reference! may be used to read the concrete value of the reference trajectory to a vector of appropriate type and length:

GenericSSMs.get_reference! — Function

get_reference!(x::AbstractVector{P}, st::CPFStorage{P})

Read the concrete reference trajectory to x from st.

ASSUMPTIONS:

traceback! has been run before invoking this function and the contents of st have not changed after that. The results from this function WILL BE WRONG, if this has not been done.
length(x) == length(st) (checked, throws error)

Unconditional simulation and prediction from SSMs

GenericSSMs.jl provides functionality for unconditional simulation and prediction from generic SSMs at state and/or observation levels. The level at which simulation or prediction should be performed, is indicated by using a parametric type Level:

GenericSSMs.Level — Type

Level{:state}, Level{:observation}, Level{(:state, :observation)}

A struct tag (struct with no fields) representing the level (state level or observation level) at which simulation or prediction should be carried out. This type is used in the functions related to simulating and predicting (functions simulate!/simulate and predict!/predict).

For example, using Level{:state} in the aforementioned functions yields predictions/simulations at the state level, and using Level{:observation} yields them at the observation level. Level{(:state, :observation)} or Level{(:observation, :state)} yields predictions/simulations at both levels simultaneously. Note that :observation may be shortened to :obs.

Unconditional simulation

Unconditional simulation means simulating a continuous slice of states $x_{l}, x_{l+1}, \ldots, x_u$ from the prior of $x_{l:u}$, where $l < u$ and $l \geq 1$ such that the initial state $x_l$ is given as an argument. Similarly, a continuous slice of observations $y_{l:u}$ may be drawn. (see below) This functionality is achieved by the following allocating and nonallocating versions of simulate. By default, $l = 1$ and $x_1$ is first drawn from $m_1$ and then conditioned on when simulating $x_2, \ldots, x_u$, but this may be changed with the argument initial:

GenericSSMs.simulate — Function

simulate(n::Integer, model::GenericSSM, ::Type{L <: Level} [, rng = Random.GLOBAL_RNG; initial::Tuple{Integer, P} = (1, m1(model, rng))])

Return a vector of length n, dest, by simulating at level L from model model. L specifies the level at which simulations should be carried out, and may be set to one of Level{:state}, Level{:observation}, Level{(:state, :observation)} (with possible abbreviations, see ?Level) which results in simulations at state, observation or state and observation levels, respectively. In the final case, each simulated value corresponds to a tuple of the form (p, o), where p is a simulated particle value, and o a simulated observation simulated given p.

The value initial = (k, x) specifies that the time index associated with dest[1] should be k, and that the state value associated with dest[1] is x. Indeed, if states are simulated x is placed to dest[1]. The default of initial corresponds to (1, m1(model, rng)) meaning that simulation begins from the first time index. See simulate! for in place version.

REQUIRES:

m1, mk, gk from the GenericSSMs interface (see ?GenericSSM) depending on input L.

GenericSSMs.simulate! — Function

simulate!(dest::AbstractVector{P}, model::GenericSSM, Level{:state}[, rng = Random.GLOBAL_RNG; initial::Tuple{Integer, P} = (1, m1(model, rng))])

Fill dest by simulating states from model model.

The value initial = (k, x) specifies:

that x should be placed to dest[1],
that the remaining elements in indices 2:length(dest) of dest should be simulated from mj where j in (k + 1):(k + length(dest) - 1).

The default of initial corresponds to (1, m1(model, rng)) meaning that the function returns a draw of state values at time indices 1:length(dest). See simulate for allocating version.

REQUIRES:

m1 and mk from the GenericSSMs interface (see ?GenericSSM).

simulate!(dest::AbstractVector{Y}, model::GenericSSM, Level{:observation}[, rng = Random.GLOBAL_RNG; initial::Tuple{Integer, P} = (1, m1(model, rng))])

Fill dest by simulating observations from model model.

The value initial = (k, x) specifies:

that x is the state that is conditioned on when simulating the first observation to dest[1],
that the remaining observations in indices 2:length(dest) of dest should be simulated conditional on states simulated from mj where j in (k + 1):(k + length(dest) - 1).

The default of initial corresponds to (1, m1(model, rng)) meaning that the function returns a draw of observations at time indices 1:length(dest). See simulate for allocating version.

REQUIRES:

m1, mk and gk from the GenericSSMs interface (see ?GenericSSM).

simulate!(dest::AbstractVector{Tuple{P, Y}}, model::GenericSSM, Level{(:state, :observation)}[, rng = Random.GLOBAL_RNG; initial::Tuple{Integer, P} = (1, m1(model, rng))])

Fill dest by simulating states and observations from model model.

The value initial = (k, x) specifies:

that x should be placed to dest[1][1],
that the remaining elements in indices 2:length(destx) of dest should be simulated from mj and gj where j in (k + 1):(k + length(dest) - 1).

The default of initial corresponds to (1, m1(model, rng)) meaning that the function returns a draw of states and observations at time indices 1:length(dest). See simulate for allocating version.

REQUIRES:

m1, mk and gk from the GenericSSMs interface (see ?GenericSSM).

Prediction

Prediction means simulating from the posterior distribution of future states $\pi(x_{n+1:n+h} \mid y_{1:n})$ or observations $\pi(y_{n+1:n+h} \mid y_{1:n})$ for some prediction horizon $h \geq 1$. Prediction requires that the particle filter has been run. The current implementation of prediction draws particles from the final particle approximation returned by the particle filter using stratified resampling and then simulates future state and/or observation trajectories given the chosen particles.

The allocating version of predict takes the form:

GenericSSMs.predict — Function

predict(st::PFStorage, newmodel::GenericSSM[, L = Level{:state}, rng = Random.GLOBAL_RNG]; nahead::Integer, nsim::Integer)::PredictStorage

Predict nahead steps ahead using nsim simulations at level L (see ?Level). st should be a PFStorage object that contains results of particle filtering n time steps (see below). newmodel should be a model struct that contains the data (if any) necessary in making the future predictions. In particular, if newmodel contains data indexed by time in the definition of mk, this function requires that said data must be indexable at the future indices n + 1, n + 2, ..., n + nahead, where n = length(st), that is, the number of time points that have been filtered. (otherwise a bounds error is thrown if bounds checking is not disabled in the definition of mk by the user) To this end, the Julia package OffsetArrays.jl offers convenient ways of defining arrays whose indexing may be offset by n (upon construction of newmodel).

The return value is an object of type PredictStorage (see ?PredictStorage), containing the simulations. See predict! for non-allocating version of this function.

REQUIRES:

mk and gk from the GenericSSMs interface (see ?GenericSSM), depending on the value of L.

ASSUMPTIONS:

pf_forward_pass! has been run with st as the first argument. Results WILL BE WRONG if this is not the case.

Prediction may also be carried out in a non-allocating fashion, by first explicitly constructing a PredictStorage object:

GenericSSMs.PredictStorage — Type

PredictStorage{S, P, L, F}

A storage object needed for performing predictions at state, observation or state and observation level.

Type parameters:

S is the type of predicted values (either the particle type, observation type or tuple consisting of the particle type and observation type).
P is the particle type.
L is the level associated with the prediction storage.
F is a floating point type.

Fields:

X::Matrix{S}: An nahead x nsim matrix for storing nsim simulated trajectories nahead steps ahead from n.
x::Vector{P}: The initial particles at time n.
w::Vector{F}: The normalised weights of the initial particles x.
n::Int: The time index associated with the weighted particles (x, w).
_u::Vector{F}, _ind::Vector{Int}: Temporaries used in predicting (best left untouched).

The easiest way to do this is via the constructor:

GenericSSMs.PredictStorage — Method

PredictStorage(st::PFStorage, model::GenericSSM, L::Type{<: Level}; nahead::Integer, nsim::Integer, n::Integer = length(st))

Construct a PredictStorage object from st and model at level L. Specifically reserve memory for predicting nahead steps ahead using nsim simulations and starting at timepoint n of the storage st.

ASSUMPTIONS:

n is less than or equal to the length of st (checked, throws error)
pf_forward_pass! has been run for st beforehand.

Note that changing the argument n above to a value less than length(st) implicitly assumes that

\[ \pi(x_{1:n} \mid y_{1:n}) \propto M_1(x_1)G_1(x_1) \prod_{k = 2}^{n} M_k(x_k \mid x_{k-1}) G_k(x_{k-1}, x_k) \ \text{ for all } x_{1:n}.\]

holds also for this new value of n (and not only for n = length(st) as discussed in Feynman-Kac representation of an SSM).

Then, predict! can be called with the storage object as the first argument:

GenericSSMs.predict! — Function

predict!(st::PredictStorage, newmodel::GenericSSM[, rng = Random.GLOBAL_RNG])

Predict in place to storage st using model newmodel. The prediction is done at the level at which st was constructed. (see ?PredictStorage)

newmodel should be a model struct that contains the data (if any) necessary in making the future predictions. In particular, if newmodel contains data indexed by time in the definition of mk, this function requires that said data must be indexable at the future indices n + 1, n + 2, ..., n + size(st, 1) at which predictions are requested. (otherwise a bounds error is thrown if bounds checking is not disabled in the definition of mk by the user)

To this end, the Julia package OffsetArrays.jl offers convenient ways of defining arrays whose indexing may be offset by n (upon construction of newmodel). See predict for allocating versions of this function.

REQUIRES: mk and gk from the GenericSSMs interface (see ?GenericSSM), depending on the value of L.

The prediction results can then be obtained from the PredictStorage object written to.