kalman

Basic implementation of the Kalman filter (and smoother).

Overview

The Kalman filter/smoother is a well-known algorithm for computing recursively the filtering/smoothing distributions of a linear Gaussian model, i.e. a model of the form:

\[\begin{split}X_0 & \sim N(\mu_0,C_0) \\ X_t & = F X_{t-1} + U_t, \quad U_t \sim N(0, C_X) \\ Y_t & = G X_t + V_t, \quad V_t \sim N(0, C_Y)\end{split}\]

Linear Gaussian models and the Kalman filter are covered in Chapter 7 of the book.

MVLinearGauss class and subclasses

To define a specific linear Gaussian model, we instantiate class MVLinearGauss (or one its subclass) as follows:

import numpy as np
from particles import kalman

ssm = kalman.MVLinearGauss(F=np.ones((1, 2)), G=np.eye(2), covX=np.eye(2),
                           covY=.3)

where the parameters have the same meaning as above. It is also possible to specify mu0`and  `cov0 (the mean and covariance of the initial state $X_0$). (See the documentation of the class for more details.)

Class MVLinearGauss is a sub-class of StateSpaceModel in module state_space_models, so it inherits methods from its parent such as:

true_states, data = ssm.simulate(30)

Class MVLinearGauss implements methods proposal, proposal0 and logeta, which correspond respectively to the optimal proposal distributions and auxiliary function for a guided or auxiliary particle filter; see Chapter 11 and module state_space_models for more details. (That the optimal quantities are tractable is, of course, due to the fact that the model is linear and Gaussian.)

To define a univariate linear Gaussian model, you may want to use instead the more conveniently parametrised class LinearGauss (which is a sub-class of ``MVLinearGauss`):

ssm = LinearGauss(rho=0.3, sigX=1., sigY=.2, sig0=1.)

which corresponds to model:

\[\begin{split}X_0 & \sim N(0, \sigma_0^2) \\ X_t|X_{t-1}=x_{t-1} & \sim N(\rho * X_{t-1},\sigma_X^2) \\ Y_t |X_t=x_t & \sim N(x_t, \sigma_Y^2)\end{split}\]

Another sub-class of MVLinearGauss defined in this module is MVLinearGauss_Guarniero_etal, which implements a particular class of linear Gaussian models often used as a benchmark (after Guarniero et al, 2016).

Kalman class

The Kalman filter is implemented as a class, Kalman, with methods filter and smoother. When instantiating the class, one passes as arguments the data, and an object that represents the considered model (i.e. an instance of MvLinearGauss, see above):

kf = kalman.Kalman(ssm=ssm, data=data)
kf.filter()

The second line implements the forward pass of a Kalman filter. The results are stored as lists of MeanAndCov objects, that is, named tuples with attributes ‘mean’ and ‘cov’ that represent a Gaussian distribution. For instance:

kf.filt[3].mean  # mean of the filtering distribution at time 3
kf.pred[7].cov  # cov matrix of the predictive distribution at time 7

The forward pass also computes the log-likelihood of the data:

kf.logpyt[5] # log-density of Y_t | Y_{0:t-1} at time t=5

Smoothing works along the same lines:

kf.smoother()

then object kf contains a list called smooth, which represents the successive (marginal) smoothing distributions:

kf.smth[8].mean  # mean of the smoothing dist at time 8

It is possible to call method smoother directly (without calling filter first). In that case, the filtering step is automatically performed as a preliminary step.

Kalman objects as iterators

It is possible to perform the forward pass step by step; in fact a Kalman object is an iterator:

kf = kalman.Kalman(ssm=ssm, data=data)
next(kf)  # one step
next(kf)  # one step

If you run the smoother after k such steps, you will obtain the smoothing distribution based on the k first data-points. It is therefore possible to compute recursively the successive smoothing distributions, but (a) at a high CPU cost; and (b) at each time, you must save the results somewhere, as attribute kf.smth gets written over and over.

Functions to perform a single step

The module also defines low-level functions that perform a single step of the forward or backward step. Some of these function makes it possible to perform such steps in parallel (e.g. for N predictive means). The table below lists these functions. Some of the required inputs are MeanAndCov objects, which may be defined as follows:

my_predictive_dist = kalman.MeanAndCov(mean=np.ones(2), cov=np.eye(2))
Function (with signature)
predict_step(F, covX, filt)
filter_step(G, covY, pred, yt)
filter_step_asarray(G, covY, pred, yt)
smoother_step(F, filt, next_pred, next_smth)

Module summary

Kalman Kalman filter/smoother.
predict_step Predictive step of Kalman filter.
filter_step Filtering step of Kalman filter.
filter_step_asarray Filtering step of Kalman filter: array version.
smoother_step Smoothing step of Kalman filter/smoother.