【カルマンフィルタ#3】【Python】階段関数を使った非ガウス型カルマンフィルタへの対応

前回からまた大分時間が経ってしまいましたが、今回は非ガウスカルマンフィルタの実装ということで、やってみたいと思います。

非線形・非ガウス型カルマンフィルタ

非線形・非ガウス型カルマンフィルタについては、こちらの動画こちらのサイトこちらの資料を見られると良いかもしれません。まあ難しい。私はこちらの本を買って読んで見たのですが、なかなか難しく、ちゃんと勉強する方には良いかもしれませんが、私のような新参ものが理解するには少々難儀します。

といっても実装したい意志はかわらないので、いろいろな資料を紐解きつつ実装してみました。

非線形・非ガウス型カルマンフィルタはこれまでの線形カルマンフィルタの手順では解けないもので、イズの分布が正規分布でなかったり、システム方程式が非線形であったりする場合に使用するものです。

私がやりたかったのは、システムノイズがコーシー分布である状態方程式で、これは通常の線形カルマンフィルタでは解けません。

非線形・非ガウス型カルマンフィルタを解く方法は色々あるみたいで、一番一般的なものは拡張カルマンフィルタという非線形の部分を線形に近似して解く方法です。ただこちらの方法で私のやりたかったことを実現することはできませんでした(私の実力不足の可能性もあります)。

その他にもアンサンブルカルマンフィルタパーティクルフィルタなどがありました。ただこの2つは私の知る限りでは乱数を発生させるため再現性がありません。また、私がやりたかったことは、何度やっても同じ結果を再現したいということでしたので、別の方法を模索しました。

そして今回行き着いたのが、階段関数によるカルマンフィルタです。詳しい解説はこちらの資料の46p以降に書いてあります。分布を細かく区分して、システムノイズ分布や観測ノイズ分布と畳み込み積分を使って事象の確率分布を計算するものです。

計算手順

簡単のため1次のトレンドモデルの場合について、手順を記載します。上記のサイトとほぼ同じなので、必要ない方は飛ばしてください。

フィルタと1期先予測は以下のとおりです。

規格化については、 各ステップの後で密度関数の全積分 が1になるように行います。

以上を繰り返します、平滑化については、今回も割愛します。

コードの作成

以上の計算をPythonで実装していきます。

1. データの作成

今回は予め作っておいたcsvを読み込みます。

今回は途中(11/3 03:00)でデータが急上昇するモデルにしました。こういったモデルの場合は正規分布では対応しきれないことがあるため、コーシー分布を使うほうが良い場合もあります。

コードは以下の通り。

これまでと違う部分は一番下の部分です。意味的には以下の図の通りで、階段関数で近似する幅と一つ一つのステップ幅を定義しています。

2. カルマンフィルタ(1期先予測・フィルタ)の作成

1期先予測とフィルタのコードです。

予測の関数内でそのまま計算する方法を記載していますが、計算を高速化するために、以下の通り予めシステムノイズの行列Qを作成しておいて、それを利用して計算しています。以下を計算することでn期先予測も1回の積で計算できます。Qの計算自体は次の過程で行っています。

その他予測の計算を高速化する方法として、畳み込み積分に高速フーリエ変換を使って計算する方法等もありますが、上記の方法のほうが早かったため今回は記載していません。詳しくはこちらのサイトを御覧ください。

3. パラメータ推定(最尤推定法)

最尤推定法自体はこれまでの方法と同じです。step配列はフィルタ分布等のx座標を格納する配列です。この方法では階段関数の幅(今回は8m)が1回の計算のたびに移動することを想定しているため、stepの値が移動していきます。calcLogLike関数の最後がその部分です。フィルタの平均値を中央としてその前後(今回は前後4m)をstep幅として移動させています。移動させるにあたって端の部分に関しては線形近似を行い値が負になるところは0にしています。このあたりの値のとり方によっては値が大きく変化するような変化がある場合には影響を受けてしまうところだと思います。色々試したのですが、この方法が安定していましたので、これにしました。本来は計算メモリが足りるのであれば、予想される範囲すべてをもとからモデル化しstepが移動しないようにするのが正確であるように思います。

4.カルマンフィルタの実行

未知数が決まったらあとはカルマンフィルタを実行するだけです。

特段変わったことはないです。

5. 図化+出力

図化と出力です。平滑化を省きました。

図化したグラフが以下です。下は急激に上昇している部分を拡大したものです。

正規分布とコーシー分布の比較

これまでの計算はシステムノイズの分布を正規分布としていましたが、コーシー分布として計算し、図化したものが以下の通り、コードで変更したのは尤度計算のqをシステムノイズ分布がコーシー分布の場合に変更したのとカルマンフィルタ実行の同様qの部分をシステムノイズ分布がコーシー分布の場合に変更しただけです。

先程の正規分布と比較すると、急激に上昇する部分でより真値に近い値にすり寄っている様子がわかると思います。これがコーシー分布にした利点です。

timestep500までの尤度は正規分布=-83.12623、コーシー分布=-58.94126でコーシー分布の方が尤度が大きい結果となりました。

今回算出したシステムノイズの正規分布(赤色)とコーシー分布(緑色)を比較したものが以下の図です。下はy軸を拡大したものです。

コーシー分布は正規分布に比べて裾野が広がっている様子がわかると思います。これのおかげて急激な上昇にも対応できるといったわけです。

まとめ

今回は非ガウス型カルマンフィルタということで、階段関数を用いたPythonでの実装をしてみました。メリットとしてはシステムノイズの分布がコーシー分布の場合でも計算できるところですね。デメリットとしては、ある程度幅を大きく取らないと急激な変化に対応できないため、先の観測値が予測できない場合には使用するのが難しいところと分布の細分化について、コーシー分布の中心ではかなり尖った分布形状になります。そのためステップ幅を小さく取らないと計算に誤差が出てしまう恐れがあると思います。ただし細かくしすぎるとメモリが大きくなって計算にも時間がかかってしまうと思います。このあたりはとても難しい問題だと思います。また、今回のように幅が観測値によって移動する場合には端の取り扱いも難しい問題になると思います。

完全に独学で実装したので、間違っているところがあると思いますが、意見をいただければ幸いです。今後はまた非線形カルマンフィルタの別の方法についても検討したいと思っておりますので、形になったら記事にしたいと思います。

参考サイト

コメント

  1. たつや より:

    カルマンフィルタのPythonコーディングに関する記事を拝読させていただき,感銘を受けています.
    初学者のため,とりあえずコードをjupyter notebookにコピペして実行してみている程度なのですが,最後の
    5. 図化+出力
    のパートで下記のエラーとなってしまいます:

    NameError Traceback (most recent call last)
    in
    11 sns.lineplot(x=x, y=y_obs, color=”gray”, ax=axes[0], label=”observation”)
    12 sns.lineplot(x=x, y=y_true, color=”black”, ax=axes[0], label=”true value”)
    —> 13 sns.lineplot(x=x, y=m[:].flatten(), color=”red”, ax=axes[0], label=”kalman filter + prediction”)
    14 #axes.set_ylim(lower, upper)
    15 #axes.set_xlim(x_min, x_max)

    NameError: name ‘m’ is not defined

    確かにそれ以前のパートで変数mは定義されていないようですので,正しい変数名をお教えいただければ幸いです.

    • ashikapengin より:

      たつやさん
      はじめましてあしかぺんぎんと申します。

      ご返信がおそくなって申し訳ございません。
      まさかこんな記事を見てくれている方がいるとは思っておりませんでした。ありがとうございます。
      ご指摘誠にありがとうございます。
      ご指摘の通りmを計算するコードが抜けておりました。
      “4.カルマンフィルタ実行”の最後にm = d_t*np.sum(f * step,1)を追加しました。
      フィルタ分布の平均値をとっているだけです。
      ご確認いただければ幸いです。
      よろしくお願いいたします。

タイトルとURLをコピーしました