估计卡尔曼滤波器周围的置信区间

Question

我一直在努力实现卡尔曼滤波器来搜索二维数据集中的异常。与我在这里找到的优秀帖子非常相似。作为下一步，我想预测下一个值将落入的置信区间（例如下限和上限值的 95% 置信度）。所以除了下面的行之外，我还想成为能够生成两条额外的线，表示下一个值将高于地板或低于天花板的置信度为 95%。

我假设我想使用卡尔曼滤波器生成的每个预测返回的不确定性协方差矩阵 (P)，但我不确定它是否正确。任何有关如何执行此操作的指导或参考将不胜感激！

kalman 2d filter in python

上面帖子中的代码随时间生成一组测量值，并使用卡尔曼滤波器对结果进行平滑处理。

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

1-sigma interval 的 2D 泛化是置信椭圆，其特征在于方程 (x-mx).T P^{-1}.(x-mx)==1，x 是参数 2D-Vector，mx 是 2D 均值或椭圆中心，P^{-1} 是逆协方差矩阵。请参阅此answer，了解如何绘制一个。与 sigma 区间一样，椭圆区域对应于真值所在的固定概率。通过使用因子n（缩放间隔长度或椭圆半径）进行缩放，可以获得更高的置信度。请注意，因子n 在一维和二维中具有不同的概率：

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

在 2D 中计算这些值有点复杂，遗憾的是我没有公开引用它。

Answer 2

如果您想要一个 95% 的区间来预测下一个值将落入，那么您需要一个预测区间而不是置信区间 (http://en.wikipedia.org/wiki/Prediction_interval)。

对于 2-D (3-D) 数据，可以通过计算数据的协方差矩阵的特征值并调整半轴的大小来计算椭圆（椭球）的半轴必要的预测概率。

有关计算 95% 预测椭圆或椭球的 Python 代码，请参阅 Prediction ellipse and prediction ellipsoid。 这可能会帮助您计算数据的预测椭圆。

Answer 3

因为您的统计数据当然来自样本，所以总体统计数据大于 2 sigma 标准差的概率为 0.5。因此，如果您没有应用 2 倍标准差的上置信度因子，我会考虑考虑是否对您期望下一个度量值低于概率 0.95 的值有良好预测的重要性。该因素的大小将取决于用于得出 0.5 总体概率的样本量。用于得出协方差矩阵的样本量越小，得出总体 0.95 统计量小于分解后样本统计量的 0.95 概率的因子就越大。