pykalman を使ってみた

はじめに

自分で作ったKalman Filter, Kalman Smoother が合っているのか検証したいという目的で、pykalman というOSSを動かしてみました。pykalman はこちら：

普通にKalman filter, Kalman soother による状態推定の計算を実行できます。pykalmanはそれだけでなく、EMアルゴリズムでモデルの要素を推定する機能もあり、面白いと思いました。ただし、実行してみると、十分な学習データがないと使うのは難しいかもと思いました。

おいらのコードはこちら：

内容

使い方メモ

バッチ処理するにはとても良く設定されています。カルマンフィルタの場合、

• 状態空間の設定(F)、システムノイズ、観測ノイズを与える
• 観測値を与えて状態の推定値、共分散行列を計算する

という使い方をします。

import numpy as np
from pykalman import KalmanFilter

F = np.array([\
        [1.0, dt],
        [0.0, 1.0]])
H = np.array([1, 0]).reshape(1, 2)
Q = np.array([\
        [1.0 / 3.0 * np.power(sig0, 2) * np.power(dt, 3), 1.0 / 2.0 * np.power(sig0, 2) * np.power(dt, 2)], \
        [1.0 / 2.0 * np.power(sig0, 2) * np.power(dt, 2), np.power(sig0, 2)* dt]])
R = np.array([sig1 * sig1]).reshape(1,1)

print(f'F {F.shape} H {H.shape}  Q {Q.shape} R {R.shape} x0 {x_init.shape} x0_cov {x_var_init.shape}')

kf = KalmanFilter(transition_matrices = F, 
                  observation_matrices = H, 
                  transition_covariance = Q, 
                  observation_covariance = R,
                  initial_state_covariance = x_var_init,
                  initial_state_mean = x_init)

# Run Kalman filter
x_est, P_est = kf.filter(pos_obs)

Kalman smoother も同じように使えます。

kf = KalmanFilter(transition_matrices = F, 
                  observation_matrices = H, 
                  transition_covariance = Q, 
                  observation_covariance = R,
                  initial_state_covariance = x_var_init,
                  initial_state_mean = x_init)
x_smooth, P_smooth = kf.smooth(pos_obs)

EM algorithm による学習

イメージですが、カルマンスムーザを実行すると、状態変数の確率分布が求められます。すなわち、ガウス分布を仮定して平均、共分散行列が求められています。これを用いて、他のパラメータを学習することができます。

試しに、初期位置をずらしてからEM学習を行い、正しい状態に収束するか見てみました。

n_iter = 200
x_init = x_init + np.array([-10,0])
kf = KalmanFilter(transition_matrices = F, 
                  observation_matrices = H, 
                  transition_covariance = Q, 
                  observation_covariance = R,
                  initial_state_covariance = x_var_init,
                  initial_state_mean = x_init,
                  em_vars = ['initial_state_mean'])


print('before learning: x_init=', kf.initial_state_mean)
x_est, P_est = kf.smooth(pos_obs)
x_em, P_em = kf.em(pos_obs, n_iter = n_iter).smooth(pos_obs)
print('after learning: x_init', kf.initial_state_mean)

動かしてみると、確かに指定したパラメータは推定できます。しかし、一つの観測系列を与えて学習させると、生成時に指定した値には収束せず、推定される状態も可適応した感じになりました。学習サンプル数を増やしてみましたが、あまり改善をしました。

In [49]: run my_test_kf/kf_pv_em_pykalman.py
before learning: x_init= [-10.   0.]
after learning: x_init [ 0.35570077 -0.13390231]

下記は一例のプロットですが、初期状態の推定値が変わる以外、ほとんど変わっていません。内部のパラメータをみると値は変わっていますが。

学習出来たらいろいろ使えるはずです。何故正解の値に収束しないかは検討が必要です。

はまりポイント

KalmanFilter のパラメータの与え方について、ドキュメントに書いてある観測行列のサイズのところに誤植がありました。pykalman.KalmanFilter? とすると表示されるメッセージの抜粋は下記のようになります。

Init signature:
pykalman.KalmanFilter(
    transition_matrices=None,
    observation_matrices=None,
    transition_covariance=None,
    observation_covariance=None,
    transition_offsets=None,
    observation_offsets=None,
    initial_state_mean=None,
    initial_state_covariance=None,
    random_state=None,
    em_vars=['transition_covariance', 'observation_covariance', 'initial_state_mean', 'initial_state_covariance'],
    n_dim_state=None,
    n_dim_obs=None,
)
Docstring:
Implements the Kalman Filter, Kalman Smoother, and EM algorithm.

...

Parameters
----------
transition_matrices : [n_timesteps-1, n_dim_state, n_dim_state] or     [n_dim_state,n_dim_state] array-like
    Also known as :math:`A`.  state transition matrix between times t and
    t+1 for t in [0...n_timesteps-2]
observation_matrices : [n_timesteps, n_dim_obs, n_dim_obs] or [n_dim_obs,     n_dim_obs] array-like
    Also known as :math:`C`.  observation matrix for times
    [0...n_timesteps-1]
transition_covariance : [n_dim_state, n_dim_state] array-like
    Also known as :math:`Q`.  state transition covariance matrix for times
    [0...n_timesteps-2]
observation_covariance : [n_dim_obs, n_dim_obs] array-like
    Also known as :math:`R`.  observation covariance matrix for times
    [0...n_timesteps-1]
transition_offsets : [n_timesteps-1, n_dim_state] or [n_dim_state]     array-like
    Also known as :math:`b`.  state offsets for times [0...n_timesteps-2]
observation_offsets : [n_timesteps, n_dim_obs] or [n_dim_obs] array-like
    Also known as :math:`d`.  observation offset for times
    [0...n_timesteps-1]
initial_state_mean : [n_dim_state] array-like
    Also known as :math:`\mu_0`. mean of initial state distribution
initial_state_covariance : [n_dim_state, n_dim_state] array-like
    Also known as :math:`\Sigma_0`.  covariance of initial state
    distribution
random_state : optional, numpy random state
    random number generator used in sampling
em_vars : optional, subset of ['transition_matrices',     'observation_matrices', 'transition_offsets', 'observation_offsets',     'transition_covariance', 'observation_covariance', 'initial_state_mean',     'initial_state_covariance'] or 'all'
    if `em_vars` is an iterable of strings only variables in `em_vars`
    will be estimated using EM.  if `em_vars` == 'all', then all
    variables will be estimated.

観測行列は、状態空間から観測ベクトルへの射影なので、[n_dim_obs, n_dim_state] となるはずですが、observation_matrix の表記は [n_dim_state, n_dim_state] となっています。当初戸惑ったのですが、こちら間違えのようで、実際に n_dim_obs x n_dim_state できちんと動作します。実際にこちらの掲示板でも話題になっていました。

尚、時系列情報として与えるのは、拡張カルマンフィルタ（EKF）のように、各時刻ごとに観測行列が変わる場合を想定した使い方をしているのだと思います。ただし、EKFの場合、事前に観測行列は分からない状態、すなわち、オンラインに推定する各状態に依存した感想行列になるので、この枠組みでは使えない気がしました。

まとめ

Kalman filter を手軽に実行するOSSとしてある pykalman を使ってみました。Kalman filter, Kalman smoother で同じ結果を得ることができました。また、パラメータ(状態遷移行列、プロセスノイズ、観測ノイズ、初期状態)をEMアルゴリズムで学習する機能の動作確認しました。とりあえず、自分は自作のを使うかもしれませんが、勉強なりました。
(2023/7/16)