Estimating confidence intervals around Kalman filter

Question

I have been working to implement a Kalman filter to search for anomalies in a two dimensional data set. Very similar to the excellent post that I found here. As a next step, I'd like to predict confidence intervals (for example 95% confidence for floor and ceiling values) for what I predict the next values will fall in. So in addition to the line below, I'd like to be able to generate two additional lines which represent a 95% confidence that the next value will be above the floor or below the ceiling.

I assume that I'll want to use the uncertainty covariance matrix (P) that is returned with each prediction generated by the Kalman filter but I'm not sure if it's right. Any guidance or reference to how to do this would be much appreciated!

kalman 2d filter in python

The code in the post above generates a set of measurements over time and uses a Kalman filter to smooth the results.

import numpy as np
import matplotlib.pyplot as plt

def kalman_xy(x, P, measurement, R,
              motion = np.matrix('0. 0. 0. 0.').T,
              Q = np.matrix(np.eye(4))):
    """
Parameters:    
x: initial state 4-tuple of location and velocity: (x0, x1, x0_dot, x1_dot)
P: initial uncertainty convariance matrix
measurement: observed position
R: measurement noise 
motion: external motion added to state vector x
Q: motion noise (same shape as P)
"""
return kalman(x, P, measurement, R, motion, Q,
              F = np.matrix('''
                  1. 0. 1. 0.;
                  0. 1. 0. 1.;
                  0. 0. 1. 0.;
                  0. 0. 0. 1.
                  '''),
              H = np.matrix('''
                  1. 0. 0. 0.;
                  0. 1. 0. 0.'''))

def kalman(x, P, measurement, R, motion, Q, F, H):
    '''
    Parameters:
    x: initial state
    P: initial uncertainty convariance matrix
    measurement: observed position (same shape as H*x)
    R: measurement noise (same shape as H)
    motion: external motion added to state vector x
    Q: motion noise (same shape as P)
    F: next state function: x_prime = F*x
    H: measurement function: position = H*x

    Return: the updated and predicted new values for (x, P)

    See also http://en.wikipedia.org/wiki/Kalman_filter

    This version of kalman can be applied to many different situations by
    appropriately defining F and H 
    '''
    # UPDATE x, P based on measurement m    
    # distance between measured and current position-belief
    y = np.matrix(measurement).T - H * x
    S = H * P * H.T + R  # residual convariance
    K = P * H.T * S.I    # Kalman gain
    x = x + K*y
    I = np.matrix(np.eye(F.shape[0])) # identity matrix
    P = (I - K*H)*P

    # PREDICT x, P based on motion
    x = F*x + motion
    P = F*P*F.T + Q

    return x, P

def demo_kalman_xy():
    x = np.matrix('0. 0. 0. 0.').T 
    P = np.matrix(np.eye(4))*1000 # initial uncertainty

    N = 20
    true_x = np.linspace(0.0, 10.0, N)
    true_y = true_x**2
    observed_x = true_x + 0.05*np.random.random(N)*true_x
    observed_y = true_y + 0.05*np.random.random(N)*true_y
    plt.plot(observed_x, observed_y, 'ro')
    result = []
    R = 0.01**2
    for meas in zip(observed_x, observed_y):
        x, P = kalman_xy(x, P, meas, R)
        result.append((x[:2]).tolist())
    kalman_x, kalman_y = zip(*result)
    plt.plot(kalman_x, kalman_y, 'g-')
    plt.show()

demo_kalman_xy()

Answer 1

The 2D generalization of the 1-sigma interval is the confidence ellipse which is characterized by the equation (x-mx).T P^{-1}.(x-mx)==1, with x being the parameter 2D-Vector, mx the 2D mean or ellipse center and P^{-1} the inverse covariance matrix. See this answer on how to draw one. Like the sigma-intervals the ellipses area corresponds to a fixed probability that the true value lies within. By scaling with the factor n (scaling the interval length or the ellipse radii) a higher confidence can be reached. Note that the Factors n have different probabilities in one and two dimensions:

|`n` | 1D-Intverval | 2D Ellipse |
==================================
  1  |  68.27%      |  39.35%    
  2  |  95.5%       |  86.47%
  3  |  99.73%      |  98.89%

Calculating these values in 2D is a bit involved and unfortunately I don't have a public reference to it.

Answer 2

If you want a 95% interval to predict the next values will fall in, then you want a prediction interval and not a confidence interval (http://en.wikipedia.org/wiki/Prediction_interval).

For 2-D (3-D) data, the semi-axes of the ellipse (ellipsoid) can be found by calculating the eigenvalues of the covariance matrix of the data and adjusting the size of the semi-axes to account for the necessary prediction probability.

See Prediction ellipse and prediction ellipsoid for a Python code to calculate the 95% prediction ellipse or ellipsoid. This might help you to calculate the prediction ellipse for your data.