Digital Reactor
モンテカルロ法

粒子フィルタ入門：非線形・非ガウスな状態空間を逐次モンテカルロで推定する

粒子フィルタ入門：非線形・非ガウスな状態空間を逐次モンテカルロで推定する

はじめに

線形でガウスならカルマンフィルタが最適解を閉じた形で与えます。では、状態の遷移が非線形で、観測が状態を歪めて返すような場合はどうするでしょうか。今回扱う成長モデルは、観測が状態の2乗に比例するため、観測からは状態の符号が分かりません。この非線形性の前では、カルマンフィルタを線形近似で拡張した拡張カルマンフィルタ（EKF）は崩れます。

同じ100ステップの系列で試すと差は明快でした。EKFの推定は真の状態に対してRMSE22.5、ところどころで±70も外れることがあります。一方、粒子フィルタは同じデータでRMSE4.6に収まり、真の軌道をなめらかに追いました。この記事では、事後分布を重み付きの粒子で近似する逐次モンテカルロの考え方を、退化とリサンプリングという要所まで含めて実装で確かめます。

対象読者:

  • 時系列の状態推定・トラッキングに関わる方
  • モンテカルロやベイズ更新の基礎があり、カルマンフィルタの適用限界を越えたい方
  • 非線形・非ガウスなモデルを扱う必要のある方

記事のポイント:

  • 状態空間モデルとフィルタリング分布の逐次更新を整理します
  • 重点サンプリングとリサンプリングで事後を粒子近似する仕組みを理解します
  • 粒子の退化と有効サンプルサイズ、リサンプリングの役割を数値で確認します

状態空間モデルとフィルタリング問題

状態が確率的に遷移し、その状態を通してノイズ込みの観測が得られます。状態遷移 xt∼f(xt∣xt−1)x_t \sim f(x_t \mid x_{t-1})、観測 yt∼g(yt∣xt)y_t \sim g(y_t \mid x_t) という枠組みで、知りたいのはフィルタリング分布 p(xt∣y1:t)p(x_t \mid y_{1:t})、つまり「今までの観測を踏まえた、今の状態の分布」です。

これは予測（遷移で1ステップ進める）と更新（観測でベイズ補正する）のベイズ再帰で書けます。線形ガウスならこの再帰が閉形式で解けてカルマンフィルタになりますが、非線形・非ガウスだと積分が解けません。そこを標本で近似するのが逐次モンテカルロ、すなわち粒子フィルタです。

今回のモデルは非線形状態空間の定番で、次の形をとります。

xt=0.5 xt−1+25 xt−11+xt−12+8cos⁡(1.2 t)+wtx_t = 0.5\, x_{t-1} + \frac{25\, x_{t-1}}{1 + x_{t-1}^2} + 8\cos(1.2\,t) + w_t yt=xt220+vty_t = \frac{x_t^2}{20} + v_t

観測が xt2x_t^2 に比例するため、yty_t からは xtx_t の符号が復元できません。事後分布は正負の2つの山を持つ多峰型になりやすく、単峰のガウスで近似するEKFには本質的に不利な設定です。

重点サンプリングから逐次モンテカルロへ

事後分布を、重み付きの粒子の集合 {x(i),w(i)}\{x^{(i)}, w^{(i)}\} で近似します。各粒子は状態空間上の1つの仮説で、重みがその尤もらしさを表します。時刻を進めるたびに、粒子を遷移モデルで動かし（予測）、観測との合致度で重みを付け替えます（更新）。これを逐次重点サンプリング（SIS）と呼びます。

SISには退化という致命的な問題があります。時間が経つと、ほとんどの粒子の重みが0に近づき、1つの粒子だけが全重みを持つようになります。粒子は名目上たくさんあっても、実際に効いているのは1つだけ、という状態です。この「実効的に効いている粒子数」を測るのが有効サンプルサイズ（ESS）で、ESS=1/∑i(w(i))2\mathrm{ESS} = 1/\sum_i (w^{(i)})^2 で定義されます。

ブートストラップ粒子フィルタ（SIR）

退化への対処がリサンプリングです。ESSが閾値（よく使うのは N/2N/2）を下回ったら、重みに比例して粒子を復元抽出し直し、重みを均等に戻します。尤度の低い粒子は淘汰され、高い粒子が複製されます。予測、重み付け、リサンプリングを繰り返すこの構成がブートストラップ粒子フィルタ（SIR）です。

リサンプリングには多項・系統・層化などの方式があり、分散の小さい系統リサンプリングがよく使われます。実装は短いです。

def systematic_resample(w, rng):
    n = len(w)
    positions = (rng.random() + np.arange(n)) / n
    return np.clip(np.searchsorted(np.cumsum(w), positions), 0, n - 1)


def particle_filter(obs, N, resample=True, seed=0):
    rng = np.random.default_rng(seed)
    particles = rng.normal(0, np.sqrt(5), N)
    w = np.ones(N) / N
    est, ess = np.zeros(T), np.zeros(T)
    for t in range(T):
        particles = f_state(particles, t) + rng.normal(0, np.sqrt(Q), N)  # 予測
        loglik = -0.5 * (obs[t] - h_obs(particles)) ** 2 / R              # 尤度
        w = w * np.exp(loglik - loglik.max()); w /= w.sum()              # 更新
        ess[t] = 1.0 / np.sum(w ** 2)
        est[t] = np.sum(w * particles)
        if resample and ess[t] < N / 2:                                  # リサンプリング
            particles = particles[systematic_resample(w, rng)]
            w = np.ones(N) / N
    return est, ess

非線形モデルで試す

真の状態軌道を既知として合成し、EKFと粒子フィルタで推定します。プロセスノイズ分散10、観測ノイズ分散1、系列長100。乱数シードは42に固定しました。

推定軌道を重ねると、EKFが観測の2乗による符号の曖昧さに翻弄されて大きく振れるのに対し、粒子フィルタは真の状態を安定して追っています。

真の状態・EKF・粒子フィルタの推定

退化とリサンプリングの効果は、ESSの推移に表れます。リサンプリングを切ると、ESSは数ステップで1近くまで落ち、以後ずっと1のまま張り付きます。全重みが1粒子に集中した退化状態です。リサンプリングを入れると、ESSは閾値を割るたびに回復し、健全な水準を保ちます。

有効サンプルサイズの推移

粒子数を変えてRMSEを見ると、増やすほど精度は上がりますが、ある水準で頭打ちになります。

手法粒子数RMSE
EKF—22.45
粒子フィルタ207.94
粒子フィルタ1004.97
粒子フィルタ5004.65
粒子フィルタ20004.65

粒子数に対するRMSE

粒子フィルタの勘所と限界

EKFがなぜ崩れたかから見ます。EKFは各ステップで関数を1次近似し、事後を単峰のガウスで表します。観測が x2x^2 で符号を潰す今回のモデルでは事後が2峰になり、単峰近似では表現できません。線形化のヤコビアンも状態に強く依存して不安定になります。非線形性が強い、あるいは事後が多峰になる問題では、EKFの前提そのものが成り立ちません。

リサンプリングは粒子フィルタの必須部品です。ESSの図が示すとおり、SISのままでは事後が1粒子に潰れ、推定は使い物になりません。ただしリサンプリングにも副作用があり、同じ粒子を何度も複製するため粒子の多様性が失われます（標本の枯渇）。ESSが閾値を割ったときだけ実施する、というのは、この多様性の消耗と退化のバランスを取るための実務的な妥協点です。

粒子数は多ければよいわけではありません。RMSEはN=500あたりで4.6程度に飽和し、それ以上増やしても下がりませんでした。この床はプロセス・観測ノイズと符号の曖昧さに由来する、モデル本来の不確実性です。計算コストは粒子数に比例するので、精度が飽和する手前を狙うのが費用対効果の分かれ目になります。

苦手な領域もはっきりしています。状態が高次元になると、必要な粒子数が次元とともに爆発し、素朴な粒子フィルタは効かなくなります。提案分布を遷移そのものに任せるブートストラップ版は実装が楽な半面、観測が鋭いと効率が落ちます。パラメータ（今回でいうノイズ分散）まで同時に推定したい場合は、粒子フィルタ単体では扱えず、粒子MCMCのような枠組みが必要です。金融の潜在ボラティリティ推定のように、状態が低次元で強く非線形な問題は、粒子フィルタが最も力を発揮する領域です。

まとめ

非線形・非ガウスな状態推定は、事後分布を重み付きの粒子で近似すれば解けます。今回のモデルでは、EKFのRMSE22.5に対し、粒子フィルタは4.6と大きく上回り、リサンプリングがなければ退化して機能しないことも数値で確認できました。手元の問題に適用するなら、まず粒子数を数百から始めてESSの推移を監視し、退化していないか、精度がどこで飽和するかを見てから本番の粒子数を決めるとよいでしょう。状態が高次元なら、素朴な粒子フィルタの前にRao-Blackwell化や次元削減を検討します。

関連記事

← 技術ブログ一覧へ