Фильтр Калмана (pykalman): значение obs_covariance и модель без перехвата

Ковариация наблюдения показывает, сколько ошибок вы предполагаете в своих входных данных. Фильтр Калмана отлично работает с нормально распределенными данными. При этом предположении вы можете использовать правило 3 сигм для расчета ковариации (в данном случае дисперсии) вашего наблюдения на основе максимальной ошибки наблюдения.

Значения в вашем вопросе можно интерпретировать следующим образом:

Пример 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()

Я не мог получить данные с помощью панд, поэтому я скачал их и прочитал из файла .

Здесь вы можете увидеть предполагаемый наклон и точку пересечения:

Чтобы проверить оценочные данные, я восстановил значение EWC из EWA, используя оценочные параметры:

О значении ковариации наблюдения

Изменяя значение ковариации наблюдения, вы сообщаете Фильтру, насколько точны входные данные (обычно вы просто описываете свою уверенность в наблюдении, используя некоторые таблицы данных или свои знания о системе).

Вот оценочные параметры и восстановленные значения 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()