【カルマンフィルタ#1】Pythonで時系列データのノイズ除去

やりたいことが色々あるんですが、全部中途半端になってます。てことで、今日は時系列データのノイズ除去を目的として、カルマンフィルタの実装をしてみたいと思います。今後オーソドックスな線形カルマンフィルタに加えて非線形カルマンフィルタについても実装方法を記事にできたらなと思ってます。

カルマンフィルタとは

カルマンフィルタは、ルドマン・カルマンによって提唱されたアルゴリズムで、時系列を表現するモデルである、状態空間モデルを逐次的に計算するものみたいです。詳しくは調べればいくらでも出てくるのでここでは、実際の実装をメインに進めていきたいと思います。

状態空間モデル

xは真値、yは観測値、wとvはノイズ、残りのG,H,Fは係数の行列で変数の書き方は教科書によって違いますが、つまりは、観測値yが得られた状態で、ノイズに紛れた本当に知りたい値である真値xを効率よく推定するアルゴリズムです。

今回はシステムモデルに1次 or 2次のトレンドモデルを切り替えて使えるようにしました。

一応1期先予測とフィルタ、平滑化の計算方法を以下に書いておきます。真値の共分散行列がVとなっておりますが、コード内ではCになっています。

実装の流れ

実装の流れを簡単に示します。ベースはこちらのサイトです。ていうかほぼここを参考にしました。それに加えてパラメータの推定や欠損値がある場合の処理を行いより実践に近いものを作ればと思います。

  1. データの作成
  2. カルマンフィルタ(1期先予測・フィルタ・平滑化)の作成
  3. パラメータ推定(最尤推定法)
  4. カルマフィルタの実行
  5. 図化+出力

プログラム作成

1. データの作成

まずは今回使うデータを作成します。実際に使用する場合はcsvとかから読み込むことになると思います。jupyter notebookで書いたままになります。適宜書き換えて使ってください。1次、2次のトレンドモデルの使い分けは9行目と10行目の一方をコメントアウトだけで大丈夫です。

必要なモジュールのインポートと初期値の設定とデータの作成になります。今回は時間1000stepのランダムウォークでデータを作成しました。以下のグラフが出力されます。真値が黒色、観測値が灰色です。ランダムのシード値は固定しました。

データの赤い点線の部分(310~360step)を見るとわかると思いますが、観測が行われなかったと仮定して欠損値としています(灰色が直線になっていますが、この間に観測値はありません)。

カルマンフィルタを使って観測値(灰色)から真値(黒色)を推定していきたいと思います。

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

この部分はまるまる上で紹介したこちらのサイトを参考にしています。

上から1期先予測、フィルタ、平滑化の関数です。これだけでカルマンフィルタを実行できるなんて簡単ですね!サイトと異なるところはフィルタに尤度を計算する部分を追加したくらいです。

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

状態空間モデルの中のシステムノイズ観測ノイズははじめの状態だと未知です。ちょうどいい値を探したいのですが、このとき使うのが最尤推定法です。詳しくはこちらのサイトが良いと思います。

関数calcLogLikeで尤度を算出します。そしてscipyモジュールの中にある最適化アルゴリズムを使ってこの尤度を最大化(最適化アルゴリズムは最小化かので関数の戻り値はマイナスにしています)するパラメータを推定します。ポイントは欠損値がある場合には1期先予測を実行しているところです。これでおそらく欠損値を考慮した尤度が算出できると思います。

算出する期間はtimesで決めています。ながければより精度は上がると思いますが、計算に時間がかかります。このくらいだと数秒から数分だと思います。

尤度計算の方法ですが、調べてると準ニュートン法の中のL-BFGS-B法ってのが多いみたいなので、上記ではそれを使っています。ただパラメータが増えてくると局所的な解に落ちることがあったので、差分進化のコードも最後に追記しています。差分進化の方法は準ニュートン法に比べて時間はかかりますが、大域的な解にたどり着けると思います。

上記の方法で1次トレンドで算出した場合、システムノイズの分散0.0114133観測ノイズの分散は0.0408188尤度-46.9124となりました。はじめデータを作った際のシステムノイズの分散0.01観測ノイズの分散0.04だったので、かなり近い値になっていると思います。



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

3でパラメータが求まりましたので、カルマンフィルタを実行します。と言っても最尤推定のときにもうやってるので、同じコードになっちゃいます。

以上でカルマンフィルタとカルマン平滑化の結果がm、Cとs、Sに格納されていると思います。

5. 図化+出力

最後の結果の図化と出力をします。

図化したグラフがこれ

これは1次のトレンドモデルの場合です。上が観測値とフィルタ+予測中が観測値と平滑化下が観測値とフィルタ+予測と真値です。tが310~360の間は上を見ると分かる通り、フィルタでは同じ値が続いていますが、中に見られる通り、平滑化を行うことでなめらかに結合されています。

下を見れば分かる通りだいぶカルマンフィルタを使うことで真値に近づいたのではないでしょうか。

まとめ

今回はPythonを使って実装メインで最尤推定法によるパラメータの推定とカルマンフィルタを実施してみました。ご質問、ご指摘ございましたらコメントいただければ幸いです。

また非線形の方法についてもやっていけたらと思いますので、よろしくお願いします。

参考文献

コメント

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