python - 如何使用 PyKalman 获得标准差？

Question

拥有一维测量数据，我想知道使用卡尔曼滤波器在每个点的状态标准偏差。我的程序如下：

from pykalman import KalmanFilter
import numpy as np

measurements = np.asarray([2, 1, 3, 6, 3, 2, 7, 3, 4, 4, 5, 1, 10, 3, 1, 5])
kf = KalmanFilter(transition_matrices=[1],
                  observation_matrices=[1],
                  initial_state_mean=measurements[0],
                  initial_state_covariance=1,
                  observation_covariance=1,
                  transition_covariance=0.01)
state_means, state_covariances = kf.filter(measurements)
state_std = np.sqrt(state_covariances[:,0])
print state_std

这导致了以下奇怪的结果：

[[ 0.70710678]
 [ 0.5811612 ]
 [ 0.50795838]
 [ 0.4597499 ]
 [ 0.42573145]
 [ 0.40067908]
 [ 0.38170166]
 [ 0.36704314]
 [ 0.35556214]
 [ 0.34647811]
 [ 0.33923608]
 [ 0.33342945]
 [ 0.32875331]
 [ 0.32497478]
 [ 0.32191347]
 [ 0.31942809]]

我预计最后一个数据点的方差会增加。我究竟做错了什么？

Answer 1

由于您提供的所有协方差矩阵（测量、转换）都很小（这意味着您不会期望观察中有太多的不确定性），因此状态协方差不能反映您渐近增加的观察离散度，因此卡尔曼滤波器输出非常顺利。但是，如果您认为测量、转换等存在更多不确定性，我认为您可以提供更高的协方差，结果您将获得 KF 输出不是很平滑（几乎遵循测量），但渐近增加将反映在KF 输出协方差也是如此，如下所示。

from pykalman import KalmanFilter
import numpy as np

measurements = np.asarray([2, 1, 3, 6, 3, 2, 7, 3, 4, 4, 5, 1, 10, 3, 1, 5])
kf = KalmanFilter(transition_matrices=[1],
                  observation_matrices=[1],
                  initial_state_mean=measurements[0],
                  initial_state_covariance=1,
                  observation_covariance=5,
                  transition_covariance=9) #0.01)
state_means, state_covariances = kf.filter(measurements)
state_std = np.sqrt(state_covariances[:,0])
print state_std
print state_means   
print state_covariances
import matplotlib.pyplot as plt
plt.plot(measurements, '-r', label='measurment')
plt.plot(state_means, '-g', label='kalman-filter output')
plt.legend(loc='upper left')
plt.show()

measurement_std = [np.std(measurements[:i]) for i in range(len(measurements))]
plt.plot(measurement_std, '-r', label='measurment std')
plt.plot(state_std, '-g', label='kalman-filter output std')
plt.legend(loc='upper left')
plt.show()