⚠️ このページは機械翻訳されています

この日本語訳は機械翻訳によるもので、誤りや不自然な表現が含まれている場合があります。正確な内容については英語版をご参照ください。

This page was machine-translated and may contain errors. Please refer to the English version for accuracy.

誤りを見つけられた場合は、joseph.austerweil@gmail.com までご連絡いただくか、GitHub リポジトリにプルリクエストをお送りください。

パーティクルフィルタリング：昨日の事後分布は今日の事前分布

この章の位置づけ（および課題についての注記）

第16章では、サンプルの集合から固定された量（平均弁当、裾確率）を推定しました。しかしデータは1つずつ到着することが多く、推定したい対象は変化します。この章では、そのストリーミングのケースに重点サンプリングを適応させます。モンテカルロ課題には必要ありませんが、第16章のモンテカルロと第18章のマルコフ連鎖モンテカルロをつなぐ自然な橋となります。そして、サンプリングが心がリアルタイムで変化する世界を追跡する方法のモデルとなるのはここです。

Chibanyの廊下を追いかける

ChibanyはラボへとつながるC長い廊下を歩いており、Jamalは自席に座りながら、数秒ごとにおおよその距離を返すドア上の雑音センサーを使ってどこにいるかを把握しようとしています。

Jamal：「ピング値がばらばらです。一つは1.3メートル、次は1.8、そして3.4へと飛ぶ。どれか1つを信頼することはできません。」
Alyssa：「信頼する必要はありません。センサーが知らない2つのことを知っています：Chibanyは前方に一定のペースで歩き続けており、各ピングは真の位置に少しノイズが加わったものです。両方を使いましょう。」
Jamal：「つまり、ピングのたびに推定を更新して、少し前の推定が次の出発点になると？」
Alyssa：「そのとおり。昨日の答えが今日の問いです。」

Alyssaが説明したのは逐次推論であり、そのためのツールは移動する目標の後ろをはい進む重み付きモンテカルロサンプルの群れです。まずパーツを整理し、それを構築しましょう。

1つずつ到着するデータ

「雑音センサー」をモデルに変えるには2つの要素が必要です。それらは状態空間モデルの標準的な道具立てです——ノイズを通じてのみ観測される隠れた変化するものです。

隠れ状態 $x_t$ — 実際に知りたいもの（時刻$t$におけるChibanyの真の位置）。直接観測することはできません。
運動モデル $p(x_t \mid x_{t-1})$ — 状態がどのように変化するか。ここでは、Chibanyはティックごとに約1メートル前進し、少しのぶれがあります：$x_t = x_{t-1} + 1 + \text{noise}$。
観測モデル $p(z_t \mid x_t)$ — 雑音のある読み取り値$z_t$が真の状態とどのように関係するか。ここでは、ピングは真の位置にセンサーノイズを加えたものです：$z_t \sim \mathcal{N}(x_t, \sigma_{\text{obs}}^2)$。

どちらのモデルも選択です——世界について我々が立てる仮定です。これらを前提として、目標はこれまでのすべてのピングを条件とした現在の状態の事後分布 $p(x_t \mid z_1, \dots, z_t)$ です。これをストリーム上で扱いやすくする鍵となるのが再帰です：

$$p(x_t \mid z_1, \dots, z_t) \propto p(z_t \mid x_t) \underbrace{p(x_t \mid z_1, \dots, z_{t-1})}_{\text{昨日から予測}}。$$

右から左に読んでください。昨日の事後分布を取り、運動モデルを通じて前進させ、Chibanyが今どこにいるかを予測します——その予測が今日の事前分布です。次に、各可能性が今日のピングをどれだけうまく説明するか、尤度 $p(z_t \mid x_t)$ で重み付けします。その積が（定数倍を除いて）今日の事後分布です。昨日の事後分布は今日の事前分布——この一文がアルゴリズム全体です。

逐次重点サンプリング

その事後分布を閉じた形式で書き下すことはできませんが、第16章ではできないときの対処法を学びました：分布をサンプルの雲で表現し、重み付けし直すことです。ここでは雲はパーティクルの集合です——$M$個の隠れ状態の候補値であり、Chibanyがどこにいるかに関するちょっとした推測のそれぞれです——そして重点サンプリングをティックごとに1回適用します。その繰り返しのストリーミング重点サンプリングは逐次重点サンプリング（SIS）と呼ばれ、もう1つのステップを加えるとパーティクルフィルタになります。

graph LR
    A["Mパーティクル<br/>（x_tの推測）"] --> B["重み付け<br/>w ∝ p(z_t | x_i)"]
    B --> C["リサンプリング<br/>重いものを残し、軽いものを捨てる"]
    C --> D["伝播<br/>x ~ p(x_t+1 | x_t)"]
    D --> A
    classDef node fill:none,stroke:#9bbcff,stroke-width:2px,color:#fff
    class A,B,C,D node
    linkStyle default stroke:#9bbcff,stroke-width:2px,color:#fff

ループには3つの動きがあり、それぞれはすでに持っているアイデアに対応しています：

重み付け。 新しいピングの尤度によって各パーティクルをスコア付けします：$w_i \propto p(z_t \mid x_i)$。真の位置の近くにあるパーティクルはピングをうまく説明して大きな重みを得ます；遠いパーティクルはほぼ何も得ません。（これはまさに第16章の重要度重みを1つの新しい観測に適用したものです。）
リサンプリング。 ここが新しい動きです。リサンプリングとは、古いパーティクルから重みに比例した確率で$M$個の新しいパーティクルを抽出することを意味します：重み$0.4$のパーティクルは平均して2回コピーされます；重み$0.001$のパーティクルはほぼ確実に消えます。具体的には、$M$個のインデックスの$\text{Categorical}(\text{正規化された重み})$による1回の抽出です。リサンプリング後、すべてのパーティクルは等しい重みを持ちます——重いものが単純に複製され、軽いものが除去されたので、群れは活動のある場所に集中します。
伝播。 各生存パーティクルを運動モデルを通じて前進させます：$x \sim p(x_{t+1} \mid x_t)$。これで雲は次のティックに向けた新鮮な推測の集合になります——昨日の事後分布が今日の事前分布になります。ループします。

以下は5つのパーティクルでループの1ティックを完全に描いたもので、すべての動きが見えます：

さらに良いのは、ループを自分で動かすことです。以下のウィジェットは、この後に出てくる作業例の廊下エピソードを1クリックごとに1フェーズずつ実行し、今行った動きを説明します：

（M = 200で完全なエピソードをステップ実行してください。次にM = 20に切り替えて何度か再実行してください——パーティクルが少ないと推定がどれほど不安定になるかを観察しましょう。その揺れは最後のセクションの「限られたメモリ」テーマを可視化したものです。）

作業例：直線上の追跡

Chibanyを実際に追ってみましょう。真の位置は$1, 2, 3, 4, 5$（ティックごとに1メートル）；センサーは雑音のあるピング$1.3, 1.8, 3.4, 3.9, 5.2$を報告します。$M = 2000$個のパーティクルを実行し、開始位置近くで初期化して、ピングごとに重み付け–リサンプリング–伝播のループを1回ステップし、毎回の重み付き平均推定値を記録します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
import jax
import jax.numpy as jnp
import jax.random as jr
from jax.scipy.stats import norm

MOTION_SD, OBS_SD = 0.3, 0.7                # wobble in the step; noise in the ping
TRUE = [1.0, 2.0, 3.0, 4.0, 5.0]
OBS  = jnp.array([1.3, 1.8, 3.4, 3.9, 5.2])

def filter_step(key, particles, z):
    kr, kp = jr.split(key)
    logw = norm.logpdf(z, particles, OBS_SD)              # 1. WEIGHT by fit to the ping
    w = jnp.exp(logw - jnp.max(logw)); w = w / jnp.sum(w)
    est = jnp.sum(w * particles)                          #    estimate = weighted mean
    idx = jr.categorical(kr, jnp.log(w), shape=(particles.shape[0],))  # 2. RESAMPLE
    survivors = particles[idx]
    moved = survivors + 1.0 + MOTION_SD * jr.normal(kp, survivors.shape)  # 3. PROPAGATE
    return moved, est

def run_filter(key, M):
    parts = 1.0 + MOTION_SD * jr.normal(key, (M,))        # init near the start (x0 ~ 1)
    ests, k = [], key
    for t in range(len(OBS)):
        k = jr.fold_in(k, t)
        parts, est = filter_step(k, parts, OBS[t])
        ests.append(round(float(est), 2))
    return ests

print("true positions :", TRUE)
print("noisy pings    :", [1.3, 1.8, 3.4, 3.9, 5.2])
print("filter estimate:", run_filter(jr.key(0), 2000))

出力：

true positions : [1.0, 2.0, 3.0, 4.0, 5.0]
noisy pings    : [1.3, 1.8, 3.4, 3.9, 5.2]
filter estimate: [1.04, 1.98, 3.1, 4.03, 5.08]

フィルタの推定値は生のピングよりもはるかに密着して真の直線$[1, 2, 3, 4, 5]$に沿います——各ピングを一定歩行の事前分布と組み合わせることでノイズを平滑化します。それがセンサーの知らない両方のモデルを使うことの成果です：

リサンプリングと退化

なぜリサンプリング（ステップ2）は省略できないのでしょうか？省略すると、プレーンな逐次重点サンプリングになります：重み付けして伝播するだけで、除去しません。問題は重みが時間とともに乗算されることです。多くのティック後、1つの幸運なパーティクルがほぼすべての重みを蓄積し、他のすべてのパーティクルは本質的にゼロを持つことになります——群れは分布のふりをする1点に崩壊しました。これは重み退化と呼ばれ、第16章の有効サンプルサイズがこれを正確に測定します。リサンプリングなしでトラックが長くなるにつれてそれが崩壊するのを観察してください：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
def make_track(key, T):
    _, kobs = jr.split(key)
    true = jnp.cumsum(jnp.ones(T))                        # 1, 2, ..., T
    return true + OBS_SD * jr.normal(kobs, (T,))          # noisy pings

def ess_no_resample(key, M, T):
    ktrack, kfilt = jr.split(key)
    obs = make_track(ktrack, T)
    parts = 1.0 + MOTION_SD * jr.normal(kfilt, (M,))
    logw, k = jnp.zeros(M), kfilt
    for t in range(T):
        k = jr.fold_in(k, t)
        logw = logw + norm.logpdf(obs[t], parts, OBS_SD)  # accumulate weights, never resample
        parts = parts + 1.0 + MOTION_SD * jr.normal(k, parts.shape)
    w = jnp.exp(logw - jnp.max(logw)); w = w / jnp.sum(w)
    return float(1.0 / jnp.sum(w**2))                     # effective sample size

for T in [5, 10, 20]:
    print(f"T = {T:2d} steps, no resampling:  ESS = {ess_no_resample(jr.key(0), 2000, T):7.1f}  / 2000")

出力：

T =  5 steps, no resampling:  ESS =   957.5  / 2000
T = 10 steps, no resampling:  ESS =   421.4  / 2000
T = 20 steps, no resampling:  ESS =    61.2  / 2000

20ステップで、2000個のパーティクルは約61個分の価値しかありません——残りは無駄な重荷です。トラックをさらに延ばすと崩壊は続きます：

リサンプリングが解決策です：ティックごとに見込みのないパーティクルを除去し、有望なものを複製することで、群れ全体を有用な状態に保ち、フィルタは崩壊することなく無限に実行できます。

プロセスモデルとしてのパーティクル

ここで、これはアルゴリズムであることを超えて心理学になります。パーティクルフィルタは完全な事後分布を保持しません——有限の少数の推測を保持し、データを左から右への厳密なパスの下で更新します。それは制限のように聞こえます。しかしそれはまた、人間が変化する世界を追跡する方法についての驚くほど良い記述でもあり、その対応こそが要点です。

小さなパーティクルフィルタが不完全な推論者について予測することを考えてみましょう：

限られたメモリ。 $M$個のパーティクルだけでは、群れは一度に限られた不確実性しか表現できません——全ての仮説の確率分布ではなく、いくつかの生きた仮説を抱えている人のように。$M$はデータにフィットできるノブになります：人は実質的にいくつの推測を行っているのか？
順序効果。 早期のデータは後のデータが正当化したであろうパーティクルを消し去ることができるため、パーティクルフィルタは証拠が届く順序に敏感です——そして人々もそうです。
試行間の変動性。 リサンプリングは確率的であるため、同じデータでの2回の実行は異なる結果になる可能性があります——人間の判断における真の変動性を反映しています。

これは合理的プロセスモデルのアイデアです：現実的なリソース予算の下で実行された同じベイズ計算が、正しい答えだけでなく、人々がそこから特徴的に逸脱する方法も予測します。パーティクルフィルタは特に、人間のカテゴリ学習、連合学習、およびオンライン文処理のモデル化に使われてきました（Sanborn, Griffiths & Navarro 2010; Austerweil & Griffiths 2013 など）——数値的便宜としてではなく、心についての仮説としてのサンプリングです。

GenJAXで

手作りのフィルタはすべての3つの動きを詳しく説明しました。GenJAXでは伝播ステップを小さな生成モデルとして書くことができます——運動モデルは@gen関数です——そしてその周りの重み付けとリサンプリングのロジックを再利用します。構造は同一です；変わるのは伝播の行だけです。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
from genjax import gen, normal as gnormal

@gen
def motion(x_prev):
    return gnormal(x_prev + 1.0, MOTION_SD) @ "x"         # the propagate step, as a model

def genjax_filter(key, M):
    parts = 1.0 + MOTION_SD * jr.normal(key, (M,))
    ests, k = [], key
    for t in range(len(OBS)):
        k = jr.fold_in(k, t); kr, kp = jr.split(k)
        logw = norm.logpdf(OBS[t], parts, OBS_SD)         # weight
        w = jnp.exp(logw - jnp.max(logw)); w = w / jnp.sum(w)
        ests.append(round(float(jnp.sum(w * parts)), 2))
        idx = jr.categorical(kr, jnp.log(w), shape=(M,))  # resample
        survivors = parts[idx]
        keys = jr.split(kp, M)
        parts = jax.vmap(lambda kk, xp: motion.simulate(kk, (xp,)).get_retval())(keys, survivors)  # propagate
    return ests

print("genjax filter estimate:", genjax_filter(jr.key(0), 2000))

出力：

genjax filter estimate: [1.04, 1.98, 3.1, 4.03, 5.08]

同じトラック、同じ答えです。（GenJAXのsmcモジュールはこのループ全体——重み付け、リサンプリング、伝播——を再利用可能な逐次モンテカルロルーティンにパッケージ化しています；ここでの手作りバージョンはまさにそれが自動化するものです。）

今できること

状態空間モデル——運動モデルと観測モデルを持つ隠れ状態——を設定し、ストリーミング推論を機能させる再帰を認識できます：昨日の事後分布は今日の事前分布。パーティクルフィルタを実行できます：重み付きパーティクルの群れで事後分布を表現し、重み付け → リサンプリング → 伝播のループを実行します。ここでリサンプリングは重いパーティクルを複製して軽いものを除去する Categorical抽出であり、重み退化を防ぎます。そして、小さなパーティクルフィルタが人間の推論の候補プロセスモデルである理由——限られたメモリ、順序効果、およびリソース予算から生じる変動性——を理解できます。

次の第18章では、目標を固定したままサンプル自体をマルコフ連鎖にします——第13章を逆向きに実行して、定常分布が求める事後分布になる連鎖を設計します。

用語集： パーティクルフィルタ, リサンプリング, 重点サンプリング, 有効サンプルサイズ。

演習

自分で試してみよう

群れを縮小する。 追跡フィルタを2000個ではなく$M = 20$個のパーティクルで再実行してください。それでも直線を追跡できますか？異なるキーで何度か実行してください——試行ごとに推定値はどれほど揺れますか？見えるものを「限られたメモリ」と「試行間の変動性」に関連付けてください。
リサンプリングをオフにする。 リサンプリングの行をスキップするようにfilter_stepを修正してください（重み付きパーティクルを直接伝播し、重みを持ち越す）。5ステップのトラックでは問題ないかもしれません；30ステップのトラックに延ばしてリサンプリングありバージョンと推定値の精度を比較してください。何が悪くなり、なぜですか？
より雑音の多いセンサー。 OBS_SDを1.4に倍増させてください（ピングの雑音がはるかに多い）。まず予測してください：フィルタの推定値は生ピングに近づくべきですか、それとも一定歩行の事前分布に近いままですか？実行して確認してください。

コンパニオンノートブックがすべてをインタラクティブに進めます：

📓 Colabで開く：17_particle_filtering.ipynb

JPPCAがこのチュートリアルシリーズを寛大にご支援いただいたことに深く感謝します。

参考文献

Austerweil, J. L., & Griffiths, T. L. (2013). A nonparametric Bayesian framework for constructing flexible feature representations. Psychological Review, 120(4), 817–851. https://doi.org/10.1037/a0034194
Sanborn, A. N., Griffiths, T. L., & Navarro, D. J. (2010). Rational approximations to rational models: Alternative algorithms for category learning. Psychological Review, 117(4), 1144–1167. https://doi.org/10.1037/a0020511

📽️ 講義より: 第7週 — Monte Carlo Methods · PDF