【问题标题】:Kalman Filter (pykalman): Value for obs_covariance and model without intercept卡尔曼滤波器(pykalman):obs_covariance 的值和没有截距的模型
【发布时间】:2018-09-10 17:31:18
【问题描述】:

我正在查看示例中显示的来自 pykalman 的 KalmanFilter:

pykalman documentation

Example 1

Example 2

我想知道

observation_covariance=100,

vs

observation_covariance=1,

文档说明

observation_covariance R: e(t)^2 ~ Gaussian (0, R)

这里应该如何正确设置值?

另外,是否可以在上面的模块中应用无截距的卡尔曼滤波器?

【问题讨论】:

    标签: transition covariance kalman-filter pykalman


    【解决方案1】:

    观察协方差显示您假设输入数据中有多少错误。卡尔曼滤波器适用于正态分布的数据。在此假设下,您可以使用 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()
    

    【讨论】:

    • Sigma 和协方差显示了信号的相同属性,但方式不同。在正态分布的情况下,它们都取决于最大误差,因为它非常接近 3*sigma 水平 (99.73%)。知道最大误差可以让您计算 sigma 和协方差(在本例中为方差)。
    • 要估计最大误差,您需要非常了解您的系统。通常,传感器的制造商会告诉您传感器的精确度(您会得到 sigma 值)。如果您没有此类信息,则需要自行定义公差范围。如果您的体重秤告诉您您的体重是 80 公斤,那么您的实际体重几乎不可能是 20 或 140 公斤,但它可以在 75 到 85 公斤之间。在这种情况下,您可以估计最大误差为 5 kg,方差 = (5/3)^2 kg^2。
    • 没错!如果你愿意,你可以从我之前的帖子中获取 python 代码并使用协方差值。该代码绘制了在 pykalman 中估计的车辆轨迹。工作数据集也在那里。 stackoverflow.com/questions/47599886/…
    • 一旦我有一些线性回归问题的代码和图表,我会在这里更新我的答案。
    • 看看我的更新。我想这就是你想要得到的,不是吗?
    猜你喜欢
    • 2021-05-06
    • 2011-04-14
    • 1970-01-01
    • 2011-12-31
    • 1970-01-01
    • 2013-11-23
    • 2017-08-11
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多