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

由于您提供的所有协方差矩阵(测量、转换)都很小(这意味着您预计观察中不会有太多不确定性)，因此状态协方差不会反射(reflect)您渐近增加的观察离散度，因此卡尔曼拟合输出非常平滑。但是，如果您认为测量、转换等存在更多不确定性，我认为您可以提供更高的协方差，因此您将得到不太平滑的 KF 输出(几乎遵循测量结果)，但渐近增加将反射(reflect)在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()