卡尔曼滤波器 (pykalman)：obs_covariance 的值和无截距的模型
Kalman Filter (pykalman): Value for obs_covariance and model without intercept
我正在查看示例中显示的来自 pykalman 的 KalmanFilter：
我想知道
observation_covariance=100,
vs
observation_covariance=1,
文档说明
observation_covariance R: e(t)^2 ~ Gaussian (0, R)
这里的值应该如何正确设置？
另外，是否可以在上述模块中应用无截距的卡尔曼滤波器？
观察协方差显示您假设输入数据中有多少误差。卡尔曼滤波器适用于正态分布的数据。在此假设下，您可以使用 3-Sigma 规则根据观察中的最大误差计算观察的协方差（在本例中为方差）。
您问题中的值可以解释如下：
示例 1
observation_covariance = 100
sigma = sqrt(observation_covariance) = 10
max_error = 3*sigma = 30
示例 2
observation_covariance = 1
sigma = sqrt(observation_covariance) = 1
max_error = 3*sigma = 3
所以你需要根据你的观测数据来选择数值。 观测越准确，观测协方差越小
另一点：您可以通过操纵协方差来调整您的过滤器，但我认为这不是一个好主意。观察协方差值越高，新观察对过滤器状态的影响越弱。
抱歉，我没看懂你问题的第二部分（关于无截距的卡尔曼滤波器）。你能解释一下你的意思吗？
您正在尝试使用回归模型，截距和斜率都属于它。
--------------------------
更新
我准备了一些代码和图来详细回答你的问题。我使用了EWC和EWA的历史数据来贴近原文
首先是代码（与上面示例中的代码完全相同，但符号不同）
from pykalman import KalmanFilter
import numpy as np
import matplotlib.pyplot as plt
# reading data (quick and dirty)
Datum=[]
EWA=[]
EWC=[]
for line in open('data/dataset.csv'):
f1, f2, f3 = line.split(';')
Datum.append(f1)
EWA.append(float(f2))
EWC.append(float(f3))
n = len(Datum)
# Filter Configuration
# both slope and intercept have to be estimated
# transition_matrix
F = np.eye(2) # identity matrix because x_(k+1) = x_(k) + noise
# observation_matrix
# H_k = [EWA_k 1]
H = np.vstack([np.matrix(EWA), np.ones((1, n))]).T[:, np.newaxis]
# transition_covariance
Q = [[1e-4, 0],
[ 0, 1e-4]]
# observation_covariance
R = 1 # max error = 3
# initial_state_mean
X0 = [0,
0]
# initial_state_covariance
P0 = [[ 1, 0],
[ 0, 1]]
# Kalman-Filter initialization
kf = KalmanFilter(n_dim_obs=1, n_dim_state=2,
transition_matrices = F,
observation_matrices = H,
transition_covariance = Q,
observation_covariance = R,
initial_state_mean = X0,
initial_state_covariance = P0)
# Filtering
state_means, state_covs = kf.filter(EWC)
# Restore EWC based on EWA and estimated parameters
EWC_restored = np.multiply(EWA, state_means[:, 0]) + state_means[:, 1]
# Plots
plt.figure(1)
ax1 = plt.subplot(211)
plt.plot(state_means[:, 0], label="Slope")
plt.grid()
plt.legend(loc="upper left")
ax2 = plt.subplot(212)
plt.plot(state_means[:, 1], label="Intercept")
plt.grid()
plt.legend(loc="upper left")
# check the result
plt.figure(2)
plt.plot(EWC, label="EWC original")
plt.plot(EWC_restored, label="EWC restored")
plt.grid()
plt.legend(loc="upper left")
plt.show()
我无法使用 pandas 检索数据，所以我下载了它们并从 file 中读取。
这里可以看到估计的斜率和截距：
为了测试估计数据，我使用估计参数从 EWA 恢复了 EWC 值：
关于观察协方差值
通过改变观察协方差值，您可以告诉过滤器输入数据的准确性（通常您只需使用一些数据表或您对系统的了解来描述您对观察结果的信心）。
这里是使用不同观测协方差值估计的参数和恢复的 EWC 值：
您可以看到滤波器更好地遵循原始函数，对观察的置信度更高（R 更小）。如果置信度较低（R 较大），则滤波器离开初始估计值（斜率 = 0，截距 = 0）非常缓慢，并且恢复的函数与原始函数相去甚远。
关于冻结拦截
如果出于某种原因想要冻结截距，则需要更改整个模型和所有滤波器参数。
在正常情况下我们有：
x = [slope; intercept] #estimation state
H = [EWA 1] #observation matrix
z = [EWC] #observation
现在我们有：
x = [slope] #estimation state
H = [EWA] #observation matrix
z = [EWC-const_intercept] #observation
结果：
代码如下：
from pykalman import KalmanFilter
import numpy as np
import matplotlib.pyplot as plt
# only slope has to be estimated (it will be manipulated by the constant intercept) - mathematically incorrect!
const_intercept = 10
# reading data (quick and dirty)
Datum=[]
EWA=[]
EWC=[]
for line in open('data/dataset.csv'):
f1, f2, f3 = line.split(';')
Datum.append(f1)
EWA.append(float(f2))
EWC.append(float(f3))
n = len(Datum)
# Filter Configuration
# transition_matrix
F = 1 # identity matrix because x_(k+1) = x_(k) + noise
# observation_matrix
# H_k = [EWA_k]
H = np.matrix(EWA).T[:, np.newaxis]
# transition_covariance
Q = 1e-4
# observation_covariance
R = 1 # max error = 3
# initial_state_mean
X0 = 0
# initial_state_covariance
P0 = 1
# Kalman-Filter initialization
kf = KalmanFilter(n_dim_obs=1, n_dim_state=1,
transition_matrices = F,
observation_matrices = H,
transition_covariance = Q,
observation_covariance = R,
initial_state_mean = X0,
initial_state_covariance = P0)
# Creating the observation based on EWC and the constant intercept
z = EWC[:] # copy the list (not just assign the reference!)
z[:] = [x - const_intercept for x in z]
# Filtering
state_means, state_covs = kf.filter(z) # the estimation for the EWC data minus constant intercept
# Restore EWC based on EWA and estimated parameters
EWC_restored = np.multiply(EWA, state_means[:, 0]) + const_intercept
# Plots
plt.figure(1)
ax1 = plt.subplot(211)
plt.plot(state_means[:, 0], label="Slope")
plt.grid()
plt.legend(loc="upper left")
ax2 = plt.subplot(212)
plt.plot(const_intercept*np.ones((n, 1)), label="Intercept")
plt.grid()
plt.legend(loc="upper left")
# check the result
plt.figure(2)
plt.plot(EWC, label="EWC original")
plt.plot(EWC_restored, label="EWC restored")
plt.grid()
plt.legend(loc="upper left")
plt.show()
我正在查看示例中显示的来自 pykalman 的 KalmanFilter：
我想知道
observation_covariance=100,
vs
observation_covariance=1,
文档说明
observation_covariance R: e(t)^2 ~ Gaussian (0, R)
这里的值应该如何正确设置？
另外，是否可以在上述模块中应用无截距的卡尔曼滤波器？
观察协方差显示您假设输入数据中有多少误差。卡尔曼滤波器适用于正态分布的数据。在此假设下，您可以使用 3-Sigma 规则根据观察中的最大误差计算观察的协方差（在本例中为方差）。
您问题中的值可以解释如下：
示例 1
observation_covariance = 100
sigma = sqrt(observation_covariance) = 10
max_error = 3*sigma = 30
示例 2
observation_covariance = 1
sigma = sqrt(observation_covariance) = 1
max_error = 3*sigma = 3
所以你需要根据你的观测数据来选择数值。 观测越准确，观测协方差越小
另一点：您可以通过操纵协方差来调整您的过滤器，但我认为这不是一个好主意。观察协方差值越高，新观察对过滤器状态的影响越弱。
抱歉，我没看懂你问题的第二部分（关于无截距的卡尔曼滤波器）。你能解释一下你的意思吗？ 您正在尝试使用回归模型，截距和斜率都属于它。
--------------------------
更新
我准备了一些代码和图来详细回答你的问题。我使用了EWC和EWA的历史数据来贴近原文
首先是代码（与上面示例中的代码完全相同，但符号不同）
from pykalman import KalmanFilter
import numpy as np
import matplotlib.pyplot as plt
# reading data (quick and dirty)
Datum=[]
EWA=[]
EWC=[]
for line in open('data/dataset.csv'):
f1, f2, f3 = line.split(';')
Datum.append(f1)
EWA.append(float(f2))
EWC.append(float(f3))
n = len(Datum)
# Filter Configuration
# both slope and intercept have to be estimated
# transition_matrix
F = np.eye(2) # identity matrix because x_(k+1) = x_(k) + noise
# observation_matrix
# H_k = [EWA_k 1]
H = np.vstack([np.matrix(EWA), np.ones((1, n))]).T[:, np.newaxis]
# transition_covariance
Q = [[1e-4, 0],
[ 0, 1e-4]]
# observation_covariance
R = 1 # max error = 3
# initial_state_mean
X0 = [0,
0]
# initial_state_covariance
P0 = [[ 1, 0],
[ 0, 1]]
# Kalman-Filter initialization
kf = KalmanFilter(n_dim_obs=1, n_dim_state=2,
transition_matrices = F,
observation_matrices = H,
transition_covariance = Q,
observation_covariance = R,
initial_state_mean = X0,
initial_state_covariance = P0)
# Filtering
state_means, state_covs = kf.filter(EWC)
# Restore EWC based on EWA and estimated parameters
EWC_restored = np.multiply(EWA, state_means[:, 0]) + state_means[:, 1]
# Plots
plt.figure(1)
ax1 = plt.subplot(211)
plt.plot(state_means[:, 0], label="Slope")
plt.grid()
plt.legend(loc="upper left")
ax2 = plt.subplot(212)
plt.plot(state_means[:, 1], label="Intercept")
plt.grid()
plt.legend(loc="upper left")
# check the result
plt.figure(2)
plt.plot(EWC, label="EWC original")
plt.plot(EWC_restored, label="EWC restored")
plt.grid()
plt.legend(loc="upper left")
plt.show()
我无法使用 pandas 检索数据，所以我下载了它们并从 file 中读取。
这里可以看到估计的斜率和截距：
为了测试估计数据，我使用估计参数从 EWA 恢复了 EWC 值：
关于观察协方差值
通过改变观察协方差值，您可以告诉过滤器输入数据的准确性（通常您只需使用一些数据表或您对系统的了解来描述您对观察结果的信心）。
这里是使用不同观测协方差值估计的参数和恢复的 EWC 值：
您可以看到滤波器更好地遵循原始函数，对观察的置信度更高（R 更小）。如果置信度较低（R 较大），则滤波器离开初始估计值（斜率 = 0，截距 = 0）非常缓慢，并且恢复的函数与原始函数相去甚远。
关于冻结拦截
如果出于某种原因想要冻结截距，则需要更改整个模型和所有滤波器参数。
在正常情况下我们有：
x = [slope; intercept] #estimation state
H = [EWA 1] #observation matrix
z = [EWC] #observation
现在我们有：
x = [slope] #estimation state
H = [EWA] #observation matrix
z = [EWC-const_intercept] #observation
结果：
代码如下：
from pykalman import KalmanFilter
import numpy as np
import matplotlib.pyplot as plt
# only slope has to be estimated (it will be manipulated by the constant intercept) - mathematically incorrect!
const_intercept = 10
# reading data (quick and dirty)
Datum=[]
EWA=[]
EWC=[]
for line in open('data/dataset.csv'):
f1, f2, f3 = line.split(';')
Datum.append(f1)
EWA.append(float(f2))
EWC.append(float(f3))
n = len(Datum)
# Filter Configuration
# transition_matrix
F = 1 # identity matrix because x_(k+1) = x_(k) + noise
# observation_matrix
# H_k = [EWA_k]
H = np.matrix(EWA).T[:, np.newaxis]
# transition_covariance
Q = 1e-4
# observation_covariance
R = 1 # max error = 3
# initial_state_mean
X0 = 0
# initial_state_covariance
P0 = 1
# Kalman-Filter initialization
kf = KalmanFilter(n_dim_obs=1, n_dim_state=1,
transition_matrices = F,
observation_matrices = H,
transition_covariance = Q,
observation_covariance = R,
initial_state_mean = X0,
initial_state_covariance = P0)
# Creating the observation based on EWC and the constant intercept
z = EWC[:] # copy the list (not just assign the reference!)
z[:] = [x - const_intercept for x in z]
# Filtering
state_means, state_covs = kf.filter(z) # the estimation for the EWC data minus constant intercept
# Restore EWC based on EWA and estimated parameters
EWC_restored = np.multiply(EWA, state_means[:, 0]) + const_intercept
# Plots
plt.figure(1)
ax1 = plt.subplot(211)
plt.plot(state_means[:, 0], label="Slope")
plt.grid()
plt.legend(loc="upper left")
ax2 = plt.subplot(212)
plt.plot(const_intercept*np.ones((n, 1)), label="Intercept")
plt.grid()
plt.legend(loc="upper left")
# check the result
plt.figure(2)
plt.plot(EWC, label="EWC original")
plt.plot(EWC_restored, label="EWC restored")
plt.grid()
plt.legend(loc="upper left")
plt.show()