实现一维卡尔曼滤波的Python平滑化

Question

我想测试Kalmnan滤波器来平滑我的一组数据。请注意，x轴的区间是不相等的。

x = [1,10,22,35,40,51,59,72,85,90,100]
y = [0.2,0.23,0.3,0.4,0.5,0.2,0.65,0.67,0.62,0.5,0.4]
plt.plot(x,y, 'go-');

其中每个点是一个观测值。很明显，x=50的点是噪声。因此，我希望卡尔曼滤波的结果是这样的。

我不是数学专家，所以我不确定这是否重要，但我的数据不是速度或位置(我找到的所有卡尔曼的例子都是指这种情况).问题是我不知道如何在Python中把这个相当简单的问题实现为卡尔曼滤波器。我看到很多人使用 pykalman 包裹

我的第一个问题是--卡尔曼滤波可以处理不相等的时间间隔吗？如果答案是否定的，那么我还是想得到一个答案，假设我的数据中的时间间隔是相等的.我也在例子中看到，数据应该是一个特定的方式，而不是像我的例子中的 "简单 "的两个列表。所以我的第二个问题是，我如何在Python中应用Kalman filtersmooth，盯着我的 "简单 "的两个列表（你可以改变x区间为相等，以防这是一个问题）。

Answer 1

卡尔曼滤波可以处理不相等的时间间隔吗？

是的，你需要注意两点--随着时间间隔的变化，你需要考虑这将对过渡矩阵（描述系统动态--这些通常会有一个delta-t依赖性）和协方差矩阵的影响--特别是过渡协方差（观测之间的时间越长，系统演变的不确定性就越大。

我不知道这是否重要，但我的数据不是速度或位置（我找到的所有卡尔曼的例子都是指这种情况）。

你可以按照自己的意愿应用卡尔曼滤波。但是请记住，卡尔曼滤波器实际上是一个状态估计器。特别是对于具有线性动力学和瓜森噪声的系统，它是一个最佳状态估计器。术语 "滤波器 "可能有点误导。如果你没有一个系统，你想代表它的动态，你需要 "编造 "一些动态来捕捉你对产生数据的物理过程的直觉理解。

很明显，x=50的点是噪声。

对我来说并不明显，因为我不知道你的数据是什么，也不知道它是如何收集的。所有的测量都会有噪声，而卡尔曼滤波器在拒绝噪声方面非常出色。在这个例子中，你似乎想做的是完全拒绝异常值。

下面是一些代码，可能有助于实现这个目标。基本上，它对每一个数据点进行多次掩蔽（忽略）的KF训练，然后通过评估其对观测协方差的影响来确定存在离群值的可能性。请注意，很可能有更好的方法来进行离群值剔除。

from pykalman import KalmanFilter
import numpy as np
import matplotlib.pyplot as plt
import copy

outlier_thresh = 0.95

# Treat y as position, and that y-dot is
# an unobserved state - the velocity,
# which is modelled as changing slowly (inertia)

# state vector [y,
#               y_dot]

# transition_matrix =  [[1, dt],
#                       [0, 1]]

observation_matrix = np.asarray([[1, 0]])

# observations:
t = [1,10,22,35,40,51,59,72,85,90,100]

# dt betweeen observations:
dt = [np.mean(np.diff(t))] + list(np.diff(t))
transition_matrices = np.asarray([[[1, each_dt],[0, 1]]
                                    for each_dt in dt])

# observations
y = np.transpose(np.asarray([[0.2,0.23,0.3,0.4,0.5,0.2,
                              0.65,0.67,0.62,0.5,0.4]]))

y = np.ma.array(y)


leave_1_out_cov = []

for i in range(len(y)):
    y_masked = np.ma.array(copy.deepcopy(y))
    y_masked[i] = np.ma.masked

    kf1 = KalmanFilter(transition_matrices = transition_matrices,
                   observation_matrices = observation_matrix)

    kf1 = kf1.em(y_masked)

    leave_1_out_cov.append(kf1.observation_covariance[0,0])

# Find indexes that contributed excessively to observation covariance
outliers = (leave_1_out_cov / np.mean(leave_1_out_cov)) < outlier_thresh

for i in range(len(outliers)):
    if outliers[i]:
        y[i] = np.ma.masked


kf1 = KalmanFilter(transition_matrices = transition_matrices,
                   observation_matrices = observation_matrix)

kf1 = kf1.em(y)

(smoothed_state_means, smoothed_state_covariances) = kf1.smooth(y)


plt.figure()
plt.plot(t, y, 'go-', label="Observations")
plt.plot(t, smoothed_state_means[:,0], 'b--', label="Value Estimate" )
plt.legend(loc="upper left")
plt.xlabel("Time (s)")
plt.ylabel("Value (unit)")

plt.show()

这就产生了下面的图。