粒子フィルターを用いたサンプリング | About SMC Samplers
粒子フィルターは 30 年前に「万能」非線型フィルタリング手法として開発されたが，それは粒子系を輸送するメカニズムとしての万能性も意味するのであり，汎用サンプラーとしても「万能」であるのかもしれないのである．近年，最適化や最適輸送の理論と結びつき，その真の力がますます明らかになりつつある．本稿では現在までのサンプラーとしての SMC 手法に対する理解をまとめる．
A Blog Entry on Bayesian Computation by an Applied Mathematician
$$
$$
1 フィルタリング問題の歴史
1.1 フィルタリング（濾波）問題とは何か？
フィルタ（濾波器）の第一義は，液体から不純物を取り除くための装置である．そのアナロジーで「フィルタリング問題」と言った場合は，信号処理の意味で電圧や電波の信号を「濾過」してノイズを除去し本当に注目したい部分を純粋化する営みのことを指す．凡ゆる通信機器において装置の熱運動によるノイズが入ることは避けられぬ自然の摂理であり，フィルタリング問題は普遍的な課題である．1
Gauss は天体観測の経験から「誤差論」と最小二乗法を発明し，これが現代の統計的推定理論の先駆けとなった．これと同じように，通信と制御の分野では「時々刻々と受信するデータから時々刻々と変化する信号をどのようにうまく濾波するか」という独自の課題から，独自の理論が発展していった．2
特に，デジタル回路がない時代では，「どのような電気回路のシステムとして濾波機をデザインすれば良いか？」という電気工学的な回路設計の問題としての側面も大きかった．
1.2 最初のフィルタリング理論
このアナログフィルタの時代で，フィルタリングの問題を統計的技術で解くための理論3 が，まず離散時間の場合が (Kolmogorov, 1941)，続いて連続時間の場合が (Wiener, 1949) によって模索された．4
しかし，この Kolmogorov と Wiener 理論では「信号とノイズの過程が定常である」という仮定の下で展開されており，この定常性の制約が Kolmogorov-Wiener 理論の広い応用を阻んでいた．
だがこれは，「当時の技術（抵抗器やコンデンサーなど）で実装出来る範囲」という制約がある上で考えられた理論としては，仕方ないことでもあった．更なる理論的発展も，デジタル技術の登場を待つ必要があった．
1.3 デジタルフィルタリングの登場
トランジスタというデジタル技術が使われるようになり，集積回路の製造技術が発達すると，「フィルタ（濾波器）」はアナログデジタル変換器，レジスタ，メモリ，マイクロプロセッサから構成されたデジタルフィルタが主に使われるようになった．アナログフィルタから物理的な姿は全く変わり，もはや肉眼では見えない装置になってしまったのである．5
その中で (Kalman, 1960) が，定常性の仮定が満たされない場合でも使えるアルゴリズムである カルマンフィルター を提案すると，すぐに Apollo 計画に導入されてスペースシャトルの制御へ実用化され，更には海軍の潜水艦などにも応用されていった．6
あらゆるシステムがデジタル化されていく中で，アナログフィルタは徹底的にデジタルフィルタに代替されるようになった．これに伴い，現代でフィルタリング問題と言った場合，現在 \(t\) までの観測 \(Y_1,\cdots,Y_t\) からノイズを除去してメッセージ部分 \(X_t\) をなるべく正確に推定するアルゴリズム という完全に数学的で抽象的な存在として研究が進められていくことになる．
1.4 状態空間モデルという語彙
Kalman の論文 (Kalman, 1960) の新規性はアルゴリズムだけではなく，状態空間モデル という枠組みを導入し，どのようなシステムに適用可能かに関する共通言語を提供したことが，即時的な応用に寄与した面もあるだろう．7 例えばこの語彙に従えば，「Kalman filter は状態空間モデルが線型正規であれば最適なフィルタリング手法」ということになる．
この状態空間モデルの枠組みと同様なものが (Baum and Petrie, 1966) によって隠れ Markov モデルとして展開され，こちらは (Baker, 1975) により音声認識に応用された．現代でも多くの音声認識システムは隠れ Markov モデルに基づく．8
状態空間モデルはその後，(Akaike, 1974) で時系列モデリングに応用され，再び統計学の分野と深く関わるようになる．9
1.5 カルマンフィルターの限界
しかし Kalman filter にも大きな制約があった．それは モデルに線型かつ正規であるという大きな制約があった ということである．一方で現実のシステムは殆どが非線型性を持つ．
そこで NASA の Ames 研究センター ではすぐに 拡張 Kalman フィルタ と共分散行列の二乗根を保持する実装が考案され，これが本当の意味で Kalman フィルターを実用に耐えるものにした (McGee and Schmidt, 1985)．10
だが，「拡張」の名前の通り本質的な解決とは言えず，システムの非線型性が強い場合は性能が伸びない．11 こうして，計算機や情報通信技術の発展と共に複雑化していくシステムに併せて，様々なフィルターが考案されていく必要があるのである．統計計算の時代の黎明である．
1.6 非線型性・非正規性という強敵
実は，1960年に開発された Kalman filter を真に超克する手法は，1990年代に入るのを待つ必要がある．他のシミュレーションに基づく ベイズ計算手法 と同じく，計算機の十分な発達を待つ必要があったのである．
というのも，Kalman filter は，線型正規状態空間モデルの下でフィルタリング分布 \[ \mathcal{L}[X_t|Y_1,\cdots,Y_t] \tag{1}\] は再び正規分布であり，平均と共分散という２つの量のみで完全に特徴付けられるため，これらのみを考慮し，微分方程式を解いて更新規則を事前に得ておくことで，計算量を大幅に削減することが出来る，というトリックに基づく．
このように線型性と正規性が同在する状況，または状態空間が有限であるなどの限られた状況でない限り，フィルタリング分布 式 1 は有限次元の十分統計量を持たない．12
従って，一般の状況に対応出来るフィルタには，上述のような計算量削減のトリックは絶対に存在せず，正面から分布 式 1 を近似する方法を考える必要がある．これには新しいアイデアが必要であると同時に，一定の性能を持つ計算機の出現を待つ必要もあったのである．
1.7 種々の近似戦略
フィルタリング分布 式 1 \[\mathbb{P}_t(X_t\in dx_t):=\mathcal{L}[X_t|Y_{0:t}]\] は Bayes の定理を通じて再帰的な関係 \[ \begin{align*} &\mathbb{P}_{t-1}(X_t\in dx_t|Y_{0:t-1})\\ &=\int_{x_{t-1}\in E}\mathbb{P}_{t-1}(X_{t-1}\in dx_{t-1}|\\ &\qquad Y_{0:t-1}=y_{0:t-1})P_t(x_{t-1},dx_t),\\ &\mathbb{P}_t(X_t\in dx_t|Y_{0:t}=y_{0:t})\\ &=\frac{1}{p_t(y_t|y_{0:t-1})}f_t(x_t|y_t)\\ &\qquad\mathbb{P}_{t-1}(X_t\in dx_t|Y_{0:t-1}=y_{0:t-1}). \end{align*} \] を満たすが，この積分を計算する必要がある．13
粒子フィルターでは相関を持った粒子系により，これらの値を逐次的に近似していくが，このアルゴリズムが出来る前に提案された近似手法を総覧する．
1.7.1 解析的な方法
節 1.5 で紹介した拡張 Kalman filter は，モデルが線型に近い場合には有効な近似手法であるが，一致推定量にはならない．この点を修正するために，線型近似の誤差を修正する IEKF (Iterated EKF) などが開発された．
1.7.2 数値積分による方法
状態空間 \(E\) の次元が3以下である場合，積分は数値積分法によっても効率的に近似できる (Kitagawa, 1987)．14
特に状態空間が有限である場合は積分は加法に退化するため，正確に実行することができる．15 この場合が隠れMarkovモデルに当たる．
隠れMarkovモデルに於て MAP 推定量を動的計画法に基づいて探索する Viterbi 算譜 (Forney, 1973), (Viterbi, 1982) とパラメータ推定のための EM アルゴリズムの変種 Baum-Welch 算譜 (Baum and Eagon, 1967), (Gopalakrishnan et al., 1989) とは音声認識の分野で広く使われており，また状態空間が近似的に離散である場合にも近似手法として採用される．16
1.7.3 Gauss混合で近似する方法
フィルタリング分布 式 1 を正規分布の有限混合で近似する方法は Gaussian sum filter と呼ばれる (Sorenson and Alspach, 1971)．これは多峰性を帯びる事後分布のモデリングに強く，物体追跡の分野で広く使われることとなったが，一般の設定に使える普遍的な手法ではなかった．17
1.7.4 アンサンブルによる近似
フィルタリング分布 式 1 を正規分布を表現する決定論的に設計された粒子系によって近似し，これをシステムモデルに従って伝播させる．これにより，フィルタリング分布の2次以下の積率までは性格に近似できていることになる．この手法を 無香 Kalman filter (Unscented Kalman filter) という (S. Julier et al., 2000), (S. J. Julier and Uhlmann, 2004)．
この手法は，宇宙から大気圏に再突入する弾道物体の追跡などの場面で，粒子フィルターよりやや正確性が劣るが計算が速く，実用的であることも知られている (Ristic et al., 2004, p. 101)．
この手法は本質的に (Evensen, 1994) の EnKF (Ensemble Kalman filter) によって拡張された（ 節 3.5 も参照）．
2 粒子フィルターの発明
線型性や正規性の仮定を 全く 必要とせず，あらゆる状態空間モデルに使えて，加えて高次元でも適用可能な夢の新フィルタリング手法は，確率的シミュレーションを利用する Monte Carlo 法の一種 であった．(Gordon et al., 1993) はこれを bootstrap filter という名前で発表し，角度観測のみを用いた物体追跡の問題 (bearings-only tracking) への応用も付した．
実は，この bearings-only tracking は極めて非線型性が強い問題として古典的なものであり，従来の拡張 Kalman filter の方法では精度が全く伸びなかった．これに比べて bootstrap filter では圧倒的な性能改善が見られたのであった．18
それだけでなく，この手法はフィルタリング問題の範疇を超えて広範な応用先を見つけつつあり，MCMC と並ぶ ベイズ計算法 となっている（ 節 2.3 ）．
なお，北川源四郎も同年（1993年）のカンファレンスにて，Monte Carlo filter の名前で同様のアルゴリズムを発表している．そのジャーナル版は (北川源四郎, 1996a)．19 日本語文献 (北川源四郎, 1996b) はウェブ上からも読める．
2.1 逐次重点サンプリングの修正としての粒子フィルター
この手法は 重点サンプリング を繰り返すという逐次重点サンプリングの改良として開発された．逐次重点サンプリングのアイデアは古く，(Hammersley and Morton, 1954) で提案され，(Mayne, 1966), (Handschin and Mayne, 1969) で逐次推定に応用された．
しかし，これには荷重の分散が指数増大するという致命的な欠点があった (Chopin, 2004)．特に何度か反復を経ると，１つの粒子を除いて他の粒子は全て荷重を殆ど持たなくなってしまうという現象が起こる (Del Moral and Doucet, 2003)．このように，荷重の分散が大きくなり，少数の粒子しか推定に関与しなくなる現象を 荷重の縮退 という．20
そこで，(Rubin, 1987) のアイデアを基に，21 リサンプリング という新たな機構を取り入れることを考える．これは 遺伝的変異・選択機構 とも呼ばれ，22 尤度の高い粒子を複製する一方で尤度の低い粒子は削除するというものである
これにより定期的に荷重をリセットすることで分散の指数増大を抑えることができ，が保たれるのである．しかしながらこの仕組みにより粒子の間に相関が生じるために，理論的解析を困難にする．この点から粒子フィルターは （平均場）相関粒子法 ともいう．23
リサンプリング機構の分計算負荷は上がるが，1990年代では計算機の性能はすでにこれを補って余りある段階に達していたのである．なお，(Gordon et al., 1993) はリサンプリング手法としては多項リサンプリングを採用しており，極めて実装が簡単という点も多くの応用を生んだ理由である．24
リサンプリングの実際の実装については，粒子フィルターの実装の稿 も参照．
2.2 粒子フィルターの応用
(Gordon et al., 1993) による粒子フィルターの考案は，物体追跡（と防衛目的）への応用が念頭にあり，この分野では粒子フィルターが極めて有効である (Ristic et al., 2004)．これは非線型性をものともしない性質に加えて，事前情報を柔軟に取り入れやすいという粒子フィルターの性格も大きく貢献している．25
コンピュータビジョンへの応用も早期から取り組まれており (Isard and Andrew, 1998)，これに続いてロボティクス，HCI (Human-Computer Interaction) 分野への応用もなされている (岡兼司, 2005), (Wills and Schön, 2023)．
一方で，(北川源四郎, 1996a) は季節調整モデルなど非定常時系列への応用が念頭にあった．確率的ボラティリティモデルなど，ファイナンスで扱う時系列は非線型性・非正規性を示すと同時にデータ数も多い．逐次推定のステップ数が増えようとも誤差が蓄積しない粒子フィルターが見事に推定を実行する．26
加えて，マクロ経済学の分野で 動学的確率的一般均衡モデル (DSGE) の推定にも応用されている．DSGE は非線型なミクロ経済学的モデルの上に構築された大規模なモデルで， 従来は MCMC を用いたベイズ推論が実行されていたが，粒子フィルターに焼戻し法や並列計算を組み合わせることでこの問題を回避でき，さらに事後分布が多峰性を持つ場合でも有用である (Herbst and Schorfheide, 2013)．27
近年では他の社会科学分野でもベイジアンモデリングが用いられ（ベイズ計算の稿 参照），粒子フィルターも エージェント・ベースド・モデル への応用などが試みられている (Lux, 2018)．
2.3 粒子フィルターの一般の推定問題への応用
こうして，粒子フィルターは非正規・非線型フィルタリング問題の解決のために開発されたアルゴリズムであったが，非線型フィルタリングに限らず極めて広い問題へと応用出来ることが徐々に明らかになった．
特にサンプラーとしても極めて有効であり（ 節 2.3.2 ），粒子フィルターは MCMC と併せて ベイズ計算 の主要トピックの１つに躍り出た (Martin et al., 2023, p. 11)．この意味で，粒子フィルターは広く 逐次 Monte Carlo 法（Sequential Monte Carlo methods 略して SMC ）とも呼ばれる．28
2.3.1 ベイズ学習
逐次的でない「静的」な設定の下での Bayes 推論に SMC を用いる方法は (Chopin, 2002) が草分け的な仕事をした．この枠組みでは，Bayes 事後分布 \(\pi(\theta|y_1,\cdots,y_N)\) の近似において，途中の \(\pi(\theta|y_1,\cdots,y_n)\) を経由して逐次的に近似することが，自然な計算コスト削減法として理解できる．
2.3.2 サンプラーとしてのSMC
複雑な分布からのサンプリングや，その正規化定数の計算という MCMC と同様の用途に SMC を使うこともできる (Del Moral et al., 2006)．
SMC は 調音 (tempering) を通じてサンプリング問題に応用される．これは目標の分布 \(\pi_p\in\mathcal{P}(E)\) に対して，これに至る \(\mathcal{P}(E)\) 上の道 \[ [p]\ni n\mapsto\pi_n\in\mathcal{P}(E) \] を通じて，より簡単な分布 \(\pi_1,\pi_2,\cdots\) から逐次的にサンプリングをするというアイデアである．
この媒介的な分布 \(\pi_n\) を焼き戻し分布 (tempered distribution) または架橋分布 (bridging distribution) などとも呼ぶ．29
詳しくは，SMC サンプラーの稿 を参照．
2.3.3 最適化
SMC を最適化へ応用することができる．
まず，任意の確率的最適化アルゴリズムに対して，これを並列して実行し，うまくいっているものとうまくいっていないものの間に遺伝的変異・選択機構を導入することでより性能の良い発見的アルゴリズムを導出出来ることを (Aldous and Vazirani, 1994) が指摘しており，この機構を “go with the winners” と呼んでいる．
古典的な大域的最適化法に 焼きなまし法 (simulated annealing) (Kirkpartick et al., 1983) があるが，(Schäfer, 2013) は特定の目的関数に対してこれを一般化して粒子法に基づく最適化法を提案し，特に多峰性を持つ場合に大きく性能を改善した．
(Johansen et al., 2008) では潜在変数モデルのパラメータの最尤推定に，EM アルゴリズム (Dempster et al., 1977) の代わりに simulated annealing に基づいた手法を用いている．分布の台が最尤推定量に収束するような分布の列 \[ \overline{\pi}_{\gamma_n}(\theta)\,\propto\,p(\theta)p(y|\theta)^{\gamma_n} \] を構成し，これから逐次的にサンプリングをするのである．最終的なアルゴリズムは，EM アルゴリズムや MCMC を用いる場合より，局所解に囚われることが少なく，初期値の設定に殆ど左右されないという利点がある．
この枠組みは一般の非凸最適化アルゴリズムになる可能性がある．30
2.3.4 稀現象シミュレーション
SMC によるサンプリングは，直接のシミュレーションが困難な分布 \(\pi_p\) に対しても，\(\pi_1,\cdots,\pi_p\) と逐次的に近似することでこれを可能にするというアイデアであった．
特に，シミュレーションが困難である分布 \(\pi_p\) の例として，極めて稀な事象 \(A_p\) に関する条件付き分布などがあり得る．これに対して事象列 \(A_0\supset\cdots\supset A_p\) を取り， \[ \pi_n(d\theta)=\frac{1}{L_n}1_{A_n}(\theta)\nu(d\theta) \] という仲介分布の列を取るのである．この問題を 稀現象シミュレーション という．
2.3.5 確率的グラフィカルモデルの推論手法として
変数間の統計的依存関係が無向グラフによって与えられるモデルを 確率的グラフィカルモデル という．ニューラルネットワーク (Rumelhart et al., 1987) も 状態空間モデル もその例である．
確率的グラフィカルモデル自体多くの応用先をもち，疫学における疾病マッピング (Green and Richardson, 2002)，画像解析 (Carbonetto and Freitas, 2003) などにも応用されている．31
このようなモデルの尤度は極めて複雑になるが，これへ至る道を自然な方法で構成することができる (Hamze and de Freitas, 2005)．複雑なモデルでは MCMC は尤度の正規化定数を評価する方法を持たないが，SMC ではこれを自然に評価することができる．
2.3.6 ポリマーシミュレーション
ポリマーシミュレーションにおいても重点サンプリング法と同じ発想が提案されたのは極めて早い段階であった (Rosenbluth and Rosenbluth, 1955)．加えてリサンプリングにあたる enrichment の考え方もあった (Wall and Erpenbeck, 1959)．
なお，このような物理・化学的文脈では，状態空間を探索する主体としての意味を強調し，「粒子」の代わりに walker と呼ぶ．32
これら２つを組み合わせることで，長いポリマー鎖のシミュレーションを可能にする方法として，相関粒子法に基づくアルゴリズム PERM (pruned-enriched Rosenbluth algorithm) が提案されている (Grassberger, 1997)．
この手法はたんぱく質の折り畳み問題にも応用されている (Hansmann and Okamoto, 1999)．
2.3.7 量子系シミュレーション
量子多体系の基底状態のシミュレーションにおけるモンテカルロ法である QMC (Quantum Monte Carlo)33 でも， branching というリサンプリング機構を取り入れた相関粒子法が用いられている (Assaraf et al., 2000)．
これは大規模な疎行列の最小固有値・固有ベクトルを近似する手法として，量子系に限らない幅広い応用がある．34
2.4 粒子フィルターの弱点
計算量
粒子フィルターは高い汎用性の代償として，多数の粒子による高精度な推論のためには多くの計算量を必要とすることは欠点に挙げられる．が，CPU や並列計算の発展により十分な量の粒子を用意できる場面も増えたため，その問題点も形骸化してきてると言える．35
荷重の縮退
また粒子フィルターは，観測 \(Y_t\) の次元が大きいなど，観測から得られる情報量が多く，尤度（ポテンシャル）の尖度が高いとき，リサンプリング機構があってもやはり縮退を起こしてしまう．36 このような場合は，観測の情報を柔軟に取り入れた提案核を構築し，誘導粒子フィルターをうまく設計する必要がある．
粒子フィルタを適用する際の課題の一つは，各粒子に割り当てられる重みが１粒子に集中する，いわゆる退化の問題を限られた数の粒子でいかに克服するかである．(上野玄太, 2019)
高次元
地球科学や天気予報の分野では \(Y_t\) は大きく（\(10^7\) を超えることもある），このような場合は粒子フィルターは実行可能でなくなる．加えて Kalman フィルターも逆行列の計算が不安定になり，アンサンブル Kalman フィルタという粒子法が用いられる（ 節 3.5 ）．
3 今後の研究
粒子フィルターの更なる応用には次の点の研究が肝要である．
3.1 Feynman-Kac 理論と粒子法の統一的理解
粒子フィルターは状態空間モデルのフィルタリングに使えるだけでなく，Feynman-Kac 測度の確率的近似に使える汎用手法である というのがより新しい数学的理解である．37 この文脈では 相関粒子法 (interacting particle methods) と呼ばれる．38
the mathematical concepts and models are now at a point where they provide a very natural and unifying mathematical basis for a large class of Monte Carlo algorithms. (Del Moral and Doucet, 2014, p. 2)
この枠組みからならば，粒子フィルターに限らず，物理学・化学・工学で用いられている多くの 発見的手法 について，理論的な解析が可能になる可能性がある．
3.2 提案分布の取り方
節 2.4 で触れた通り，特に観測の情報量が大きい場合，提案分布の選び方が粒子フィルターの精度を大きく左右する．これは粒子を高確率領域に始めから誘導するように設計する 誘導粒子フィルター によってある程度対処できる．39
参照 Markov 核 \(M_t\) がより良い攪拌性を持ち，ポテンシャル \(G_t\) が平坦であるほど性能は良い傾向があるが，40 このときの提案分布の取り方について普遍的な指針というものが得られていない．
The key factors for a successful application of particle filters in practice are therefore a good choice of the importance density and Rao-Blackwellization if possible. (Ristic et al., 2004) Epilogue
(Guarniero et al., 2017) と (Heng et al., 2020) は提案分布にパラメトリックモデルを用意し，粒子推定量の分散を最小化するようにそのモデル内で逐次的に最適化していく機構を提案している．
(Naesseth et al., 2015) が提唱する nested SMC は，は各時刻での提案分布を近似するために，もう一つのSMCを内部に走らせる．当然計算量は二倍になるが，それでも単純な bootstrap filter から大きく性能が改善する場合が多い．
加えて機械学習の観点からも，真の事後分布との KL 距離を適応的に最小化する汎用手法に，提案分布のパラメトリックモデルにニューラルネットワークによる大型のパラメトリックモデルを併せる手法が提案されている (Gu et al., 2015)．
3.3 高次元性への対応
状態空間が高次元になることと，既存のモデルを状態空間モデルに定式化することに困難が伴うことが，朧げながら共通課題のように思われるが，その現れ方と解決法は個々の事例で異なる．
3.3.1 部分的な線型構造の利用
周辺化粒子フィルター，または Rao-Blackwellized particle filter とも呼ばれる方法である．これは多くの場面で，仮定されている状態空間モデルが部分的に線型である場合に，線型の部分をKalmanフィルターによって正確に解き，残った部分のみを粒子フィルターで解くことで，精度の向上とアルゴリズムの効率化を図る方法である．
(Ristic et al., 2004, p. 287) では，探知前追跡 (track-before-detect) の問題41において，周辺化粒子フィルタが，必要な粒子数を大きく削減してくれることを紹介している．
その他にも，問題毎の特有の構造を利用して計算量を削減・パフォーマンスを最適化することは重要な営みである．
3.4 漸近論
3.4.1 粒子数に関する漸近論
目前の問題を，所与の精度で解くために必要な粒子数は幾ばくか？という問題は実用上も有益だと思われる．42
3.4.2 時間離散化に関する漸近論
3.5 EnKFの数学的解析
状態空間が高次元である場合，（どうなる？）
モデルが正規性を持つならば，EnKF (Ensemble Kalman Filter) (Evensen, 1994) が有効であり，データ同化の分野で広く使われている．43
一方で，その数学的な振る舞いはほとんど解明されておらず，1次元の場合から調べている状況である (Del Moral and Horton, 2023)．
3.6 粒子フィルターの応用
リアルタイム性と ベイズ推定の意思決定への応用との相性の良さ が，今後多くの重要な応用を見る可能性がある．
3.6.1 属人化医療への応用
筆者は属人化医療への応用が大きなモチベーションになっている．44
病気の進行（あるいは健康）のモニタリングのために，健康診断やウェアラブルデバイス，フォローアップから得られるデータは，治療の見直しや異常の早期発見のために即時処理されることが望ましい．これに逐次モンテカルロ法でBayes的に迫る研究がある (Alvares et al., 2021) ！
さらに，個々人の日常生活のレベルではSMCを用いているものはどうやらまだなく，運動と睡眠時間の間の関係と処置効果をMonte Carlo法を用いて推定している研究はある (Daza and Schneider, 2022)．これを逐次化することで，よりリアルタイムで自分に合った生活習慣への示唆が得られるアプリを開発できるかもしれない．
また，属人化医療においてはシステム生物学的なモデルに基づいた薬効推定が欠かせない．小規模な患者群と測定時間（服薬時間）に関する不確定性を考慮した手法が (Krengel et al., 2013) で考察されている．
また，腫瘍サンプルに含まれる体細胞の突然変異に関するデータから，SMCを用いて腫瘍の発達と進行の状態を理解する手法も提案されている (Ogundijo et al., 2019)．
3.6.2 ゲノム解析への応用
次世代DNAシークエンサーでは，DNAの各塩基ごとに異なる蛍光物質を結合させ，蛍光の波長と強度により塩基を読み取る仕組みであり，蛍光強度の生データからDNA配列データへ変換するベースコールと呼ばれる段階で粒子フィルターを使うことも提案されている (Shen and Vikalo, 2012)．
3.6.3 疫学への応用
Covid-19のようなパンデミックにおいて，疫学モデルを通じて時間変動する再生産数をリアルタイムでモニタリングをして意思決定に繋げるためのSMC手法も考えられており，実際にノルウェーで使用され有効性が実証された (Strovik et al., 2023)．
References
Akaike, H. (1974). Markovian representation of stochastic processes and its application to the analysis of autoregressive moving average processes. Annals of the Institute of Statistical Mathematics, 26, 363–387.
Aldous, D., and Vazirani, U. (1994). "Go with the winners" algorithms. In Proceedings 35th annual symposium on foundations of computer science, pages 492–501.
Alvares, D., Armero, C., Forte, A., and Chopin, N. (2021). Sequential monte carlo methods in bayesian joint models for longitudinal and time-to-event data. Statistical Modelling, 21(1-2), 161–181.
Anderson, B. D. O., and Moore, J. B. (1979). Optimal filtering. Dover Publications.
Assaraf, R., Caffarel, M., and Khelif, A. (2000). Diffusion monte carlo methods with a fixed number of walkers. Physical Review E, 61, 4566.
Bain, A., and Crisan, D. (2009). Fundamentals of stochsatic filtering. Springer New York.
Baker, J. (1975). The DRAGON system–an overview. IEEE Transaction on Acoustics, Speech, and Signal Processing, 23(1), 24–29.
Baum, L. E., and Eagon, J. A. (1967). An inequality with applications to statistical estimation for probabilistic functions of markov processes and to a model for ecology. Bulletin of the American Mathematical Society, 73(3), 360–363.
Baum, L. E., and Petrie, T. (1966). Statistical inference for probabilistic functions of finite state markov chains. The Annals of Mathematical Statistics, 37(6), 1554–1563.
Behrens, G., Friel, N., and Hurn, M. (2012). Tuning tempered transitions. Statistics and Computing, 22, 65–78.
Cappé, O., Moulines, E., and Rydén, T. (2005). Inference in hidden markov models. Springer New York.
Carbonetto, P., and Freitas, N. de. (2003). Why can’t josé read? The problem of learning semantic associations in a robot environment. In Proceedings of the HLT-NAACL 2003 workshop on learning word meaning from non-linguistic data, pages 54–61.
Chopin, N. (2002). A sequential particle filter method for static models. Biometrika, 89(3).
Chopin, N. (2004). Central limit theorem for sequential monte carlo methods and its application to bayesian inference. The Annals of Statistics, 32(6), 2385–2411.
Chopin, N., and Papaspiliopoulos, O. (2020). An introduction to sequential monte carlo. Springer Cham.
Chopin, N., Singh, S. S., Soto, T., and Vihola, M. (2022). On Resampling Schemes for Particle Filter with Weakly Informative Observations. The Annals of Statistics, 50(6), 3197–3222.
Creal, D. (2011). A survey of sequential monte carlo methods for economics and finance. Econometric Reviews, 31(3), 245–296.
Crisan, D., and Doucet, A. (2002). A survey of convergence results on particle filtering methods for practitioners. IEEE Transactions on Signal Processing, 50(3), 736–746.
Daza, E. J., and Schneider, L. (2022). Model-twin randomization (MoTR): A monte carlo method for estimating the within-individual average treatment effect using wearable sensors.
Del Moral, P. (2013). Mean field simulation for monte carlo integration,Vol. 126. Chapman; Hall/CRC.
Del Moral, P., and Doucet, A. (2003). Séminaire de probabilités XXXVII. In J. Azéma, M. Émery, M. Ledoux, and M. Yor, editors,Vol. 1832, pages 415–446. Springer Berlin, Heidelberg.
Del Moral, P., and Doucet, A. (2014). Particle methods: An introduction with applications. ESAIM: Proceedings and Surveys, 44, 1–46.
Del Moral, P., Doucet, A., and Jasra, A. (2006). Sequential Monte Carlo Samplers. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(3), 411–436.
Del Moral, P., and Horton, E. (2023). A theoretical analysis of one-dimensional discrete generation ensemble kalman particle filters. The Annals of Applied Probability, 33(2), 1327–1372.
Del Moral, P., and Penev, S. (2014). Stochastic processes. Chapman; Hall/CRC.
Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society. Series B (Methodological), 39(1), 1–22.
Doucet, A., Freitas, N., and Gordon, N. (Eds.). (2001). Sequential monte carlo methods in practice. Springer New York.
Evensen, G. (1994). Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics. Journal of Geophysical Research, 99(C5), 10143–10162.
Forney, G. D. (1973). The vitebi algorithm. Proceedings of the IEEE, 61(3), 268–278.
Gales, M., and Young, S. (2007). The application of hidden markov models in speech recognition. Foundations and Trends in Signal Processing, 1(3).
Gopalakrishnan, P. S., Kanevsky, D., Nadas, A., and Nahamoo, D. (1989). A generalization of the baum algorithm to rational objective functions. In IEEE conference on acoustics, speech, and signal processing,Vol. 1, pages 631–634.
Gordon, N. G., Salmond, D. J., and Smith, A. F. M. (1993). Novel approach to nonlinear/non-gaussian bayesian state estimation. IEE Proceedings-F, 140(2), 107–113.
Grassberger, P. (1997). Pruned-enriched rosenbluth method: Simulations of \(\theta\)-polymers of chain length up to 1 000 000. Physical Review E, 56, 3682–3693.
Green, P. J., and Richardson, S. (2002). Hidden markov models and disease mapping. Journal of the American Statistical Association, 97(460), 1055–1070.
Gu, S., Ghahramani, Z., and Turner, R. E. (2015). Neural adaptive sequential monte carlo. In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett, editors, Advances in neural information processing systems,Vol. 28. Curran Associates, Inc.
Guarniero, P., Johansen, A. M., and Lee, A. (2017). The iterated auxiliary particle filter. Journal of the American Statistical Association, 112(520).
Hammersley, J. M., and Morton, K. W. (1954). Poor man’s monte carlo. Journal of the Royal Statistical Society. Series B (Methodological), 16(1), 23–38.
Hamze, F., and de Freitas, N. (2005). Hot coupling: A particle approach to inference and normalization on pairwise undirected graphs. Advances in Neural Information Processing Systems, 18.
Handschin, J. E., and Mayne, D. Q. (1969). Monte carlo techniques to estimate the conditional expectation in multi-stage non-linear filtering. International Journal of Control, 9(5), 547–559.
Hansmann, U. H. E., and Okamoto, Y. (1999). New monte carlo algorithms for protein folding. Current Opinion in Structural Biology, 9(2), 177–183.
Heng, J., Bishop, A. N., Deligiannidis, G., and Doucet, A. (2020). Controlled sequential monte carlo. The Annals of Statistics, 48(5), 2904–2929.
Herbst, E. P., and Schorfheide, F. (2013). Sequential monte carlo sampling for DSGE models. NBER Working Paper Series, 19152.
Iba, Y. (2001). Population monte carlo algorithms. 人工知能学会論文誌, 16(2), 279–286.
Isard, M., and Andrew. (1998). CONDENSATION–conditional density propagation for visual tracking. International Journal of Computer Vision, 29, 5–28.
Johansen, A. M., Doucet, A., and Davy, M. (2008). Particle methods for maximum likelihood estimation in latent variable models. Statistics and Computing, 18, 47–57.
Julier, S. J., and Uhlmann, J. K. (2004). Unscented filtering and nonlinear estimation. Proceedings of the IEEE, 92(3), 401–422.
Julier, S., Uhlmann, J., and Durrant-Whyte, H. F. (2000). A new method for the nonlinear transformation of means and covariances in filters and estimators. IEEE Transactions on Automatic Control, 45(3), 477–482.
Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1), 35–45.
Kirkpartick, S., Gelatt, C. D., and Vecchi, M. P. (1983). Optimization by simulated annealing. Science, 220(4598), 671–680.
Kitagawa, G. (1987). Non-gaussian state-space modeling of nonstationary time series. Journal of the American Statistical Association, 82(400), 1032–1041.
Kitagawa, G. (1998). A self-organizing state-space model. Journal of the American Statistical Association, 93(443), 1203–1215.
Kolmogorov, A. N. (1941). Interpolation and extrapolation of stationary sequences. Bulletin de l-Acad´emie Des Sciences de U.S.S.R., Ser. Math., 5, 3–14.
Kremer, K., and Binder, K. (1988). Monte carlo simulation of lattice models for macromolecules. Computer Physics Reports, 7(6), 259–310.
Krengel, A., Hauth, J., Taskinen, M., Adiels, M., and Jirstrand, M. (2013). A continuous-time adaptive particle filter for estimations under measurement time uncertainties with an application to a plasma-leucine mixed effects model. BMC Systems Biology, 7(8).
Kummert, J., Schulz, A., Redick, T., Ayoub, N., Modabber, A., Abel, D., and Hammer, B. (2021). Efficient reject options for particle filter object tracking in medical applications. Sensors, 21(6), 2114.
Lux, T. (2018). Estimation of agent-based model using sequential monte carlo methods. Journal of Economic Dynamics and Control, 91, 391–408.
Martin, G. M., Fraizier, D. T., and Robert, C. P. (2023). Computing bayes: From then ‘til now. Statistical Science, Advanced Publication, 1–17.
Mayne, D. Q. (1966). A solution of the smoothing problem for linear dynamic systems. Automatica, 4(2), 73–92.
McGee, L. A., and Schmidt, S. F. (1985). Discovery of the kalman filter as a practical tool for aerospace and industry. NASA Technical Memorandum, 86847.
Naesseth, C., Lindsten, F., and Schon, T. (2015). Nested sequential monte carlo methods. Proceedings of the 32nd International Conference on Machine Learning Research, 37, 1292–1301.
Ogundijo, O. E., Zhu, K., Wang, X., and Anastassiou, D. (2019). A sequential monte carlo algorithm for inference of subclonal structure in cancer. PLoS One, 14(1).
Rabiner, L. R. (1989). A tutorial on hidden markov models and selected applications in speech recognition. Proceedings of the IEEE, 77(2), 257–286.
Ristic, B., Arulampalam, S., and Gordon, N. (2004). Beyond the kalman filter: Particle filters for tracking applications. Artech House Publishers: Boston, London.
Robert, C. P., and Casella, G. (2004). Monte carlo statistical methods. Springer New York.
Rosenbluth, M. N., and Rosenbluth, A. W. (1955). Monte carlo calculation of the average extension of molecular chains. The Journal of Chemical Physics, 23(2), 356–359.
Rubin, D. B. (1987). A noniterative sampling/importance resampling alternative to the data augmentation algorithm for creating a few imputations when the fraction of missing information is modest: The SIR algorithm (discussion of tanner and wong). Journal of the American Statistical Association, 82(398), 543–546.
Rumelhart, D. E., Hinton, G. E., and Williams, R. J. (1987). Parallel distributed processing: Explorations in the microstructure of cognition: foundations. In D. E. Rumelhart and J. L. McClelland, editors, pages 318–362. MIT Press.
Schäfer, C. (2013). Particle algorithms for optimization on binary spaces. AMC Transactions on Modeling and Computer Simulation, 12(1), 1–25.
Shen, X., and Vikalo, H. (2012). ParticleCall: A particle filter for base calling in next-generation sequencing systems. BMC Bioinformatics, 13(160).
Sorenson, H. W., and Alspach, D. L. (1971). Recursive bayesian estimation using gaussian sums. Automatica, 7(4), 465–479.
Strovik, G., Palomares, A. D.-L., Engebretsen, S., Rø, G. Ø. I., Engø-Monsen, K., Kristoffersen, A.-B., … Frigessi, A. (2023). A sequential monte carlo approach to estimate a time varying reproduction number in infectious disease models: The covid-19 case. Journal of the Royal Statistical Society Series A: Statistics in Society, qnad043.
Viterbi, A. J. (1982). Error bounds for convolutional codes and an asymptotically optimum decoding algorithm. IEEE Transactions on Information Theory, 13, 260–269.
Wall, F. T., and Erpenbeck, J. J. (1959). Statistical computation of radii of gyration and mean internal dimensions of polymer molecules. The Journal of Chemical Physics, 30(3), 637–640.
Wiener, N. (1949). Extrapolation, interpolation & smoothing of stationary time series. The MIT Press, Cambridge, Mass.
Wills, A. G., and Schön, T. B. (2023). Sequential monte carlo: A unified review. Annual Review of Control, Robotics, and Autonomous Systems, 6, 159–182.
Yang, J., Yao, Y., and Yang, D. (2023). Particle filter based on harris hawks optimization algorithm for underwater visual tracking. Journal of Marine Science and Engineering, 11(7), 1456.
上野玄太. (2019). 粒子フィルタとデータ同化. 統計数理, 67(2), 241–253.
北川源四郎. (1996a). Monte carlo filter and smoother for non-gaussian nonlinear state space models. Journal of Computational and Graphical Statistics, 5(1), 1–25.
北川源四郎. (1996b). モンテカルロ・フィルタおよび平滑化について. 統計数理, 44(1), 31–48.
岡兼司. (2005). 適応的拡散制御を伴うパーティクルフィルタを用いた頭部姿勢推定システム. 電子情報通信学会論文誌, 88(8), 1601–1613.
有本卓. (1970). カルマン・フィルタ. 計測と制御, 9(8), 590–598.
矢野浩一. (2014). 粒子フィルタの基礎と応用. 日本統計学会誌, 44(1), 189–216.
青山友紀. (1986). ディジタルフィルタ. 計測と制御, 6(25).
Footnotes
(Anderson and Moore, 1979) 第1.1節．↩︎
このフィルタリング問題の統計的な側面を，前述の電気工学的な側面から区別して，stochastic filteringと呼んだりもする．↩︎
(Bain and Crisan, 2009) 1.3節．↩︎
(青山友紀, 1986) 「当初は大型コンピュータを用いたシミュレーションの技法であったディジタルフィルタが，今では超LSI 1チップで実現されるまでになった．」他にも，発表時1986年では，音声信号を中心とする低周波帯域ではデジタルフィルタがアナログフィルタを駆逐しつつあること，通信システムのデジタル化に伴ってこの勢いは完全に代替するまで進むだろうとの筆者の考えが述べられている．↩︎
(McGee and Schmidt, 1985) がNASAからの資料．また，(Del Moral and Penev, 2014) も参照．加えて，(Kalman, 1960) からちょうど10周年の文献 (有本卓, 1970) に当時の雰囲気も感じさせる良いサーベイがある．↩︎
(Ristic et al., 2004, p. 3) “The state-space approach is convenient for handling multivariate data and nonlinear/non-Gaussian processes and it provides a significant advantage over traditional time-series techniques for these problems.”↩︎
(Rabiner, 1989)．(Gales and Young, 2007) には “almost all present day large vocabulary continuous speech recognition (LVCSR) systems are based on HMMs” とある．↩︎
(Kitagawa, 1998) で指摘されている．↩︎
また，同様にアポロ計画の中で，飛行体に載積出来るような小規模な計算資源と短い単語長でも安定してKalman filterが動くように，共分散行列の二乗根を保持するという “square-root” formulation of the filter が 1972年に考案されたことが，summary でも触れられている (McGee and Schmidt, 1985)．↩︎
(S. J. Julier and Uhlmann, 2004, p. 402) など．拡張 Kalman filter は遷移関数を線型近似することに基づく．よって Jacobi 行列の計算が必要であり，このため 解析的方法 とも呼ばれる (Ristic et al., 2004, p. 21)．よって遷移関数が可微分でない場合も実行不可能である．とはいっても，拡張 Kalman filter はナビゲーションシステムや GPS のデファクトスタンダードである Wikipedia．(Ristic et al., 2004) 第7章では距離のみでの追跡 (range-only tracking) では拡張 Kalman filter の性能と大差なく，計算量の問題から EKF の選択を推奨している．第8章でも弾道物体追跡の問題で，限られた設定では粒子フィルターと同等の性能を見せている．↩︎
この結果は (Chopin and Papaspiliopoulos, 2020) 第5章 など．記法 \(Y_{0:t}\) は 本サイトの数学記法一覧 を参照↩︎
(Crisan and Doucet, 2002) に3の数字が例示されている．↩︎
(Chopin and Papaspiliopoulos, 2020) 第6章，(Ristic et al., 2004) 2.2節．↩︎
(Gordon et al., 1993) の結果であると同時に，(Ristic et al., 2004) 第6章でも種々の手法と比較した数値実験がなされている．一方第7章にて，距離のみでの追跡 (range-only tracking) では拡張 Kalman filter の性能と大差なく，計算量の問題から EKF の選択を推奨している．↩︎
Del MoralのWebサイトも歴史的背景に詳しい．↩︎
weight degeneracy (Creal, 2011) p.253，(Cappé et al., 2005) 第7章，(Robert and Casella, 2004, p. 551) 14.3.3節 など．↩︎
SIR (Sampling/Importance Resampling) Algorithm と呼ばれるものであった．これの逐次化が粒子フィルターだとみなせる (Robert and Casella, 2004, p. 552) 14.3.4節．↩︎
(Del Moral and Doucet, 2014) などの用語である．↩︎
(Del Moral and Horton, 2023) では mean-field type interacting particle methods と呼んでいる．呼び方については (Iba, 2001) と (Del Moral, 2013) と 紹介記事 も参照．↩︎
(Creal, 2011) p.256．↩︎
事前情報というのは，自動運転の文脈では自動車が動き得る領域というのは極めて限られている，というような事前に判明しているが，うまく取り入れにくい情報のことをいう．(Ristic et al., 2004) 第6章や第9章で繰り返し種々の設定で実証されている．(Yang et al., 2023) は水中での物体追跡が縮退により従来は粒子フィルタが使えなかった問題の解決を試みている．(Kummert et al., 2021) はロボット支援を用いた手術中に物体追跡を利用する際に，尤度が低すぎるなどの要素から追跡対象を失った状態を検出してアラートを出す機構を開発している．↩︎
ファイナンスにおける確率的ボラティリティモデルなどの例が挙げられている (Creal, 2011) p.256．↩︎
粒子フィルターの経済学での応用が増えたきっかけが Fernández-Villaverde and Rubio-Ramírez (2005, 2007) による（小規模な）DSGEモデルの推定への応用だった (Creal, 2011, p. 246)．(矢野浩一, 2014) は実物景気循環モデル (Real Business Cycles Model) への応用を解説している．↩︎
(Crisan and Doucet, 2002) ではすでに SMC とも呼ばれることが記されている．(Chopin and Papaspiliopoulos, 2020) 第3章にSMCのフィルタリング以外の多くの応用が紹介されている．↩︎
(Behrens et al., 2012) はいずれも用いている．(Hamze and de Freitas, 2005) に tempered distribution の語用法がある，↩︎
疾病マッピングは疾病の空間的な分布を把握すること，画像解析とは画像データから特徴を抽出してその意味論を理解することを指す．↩︎
(Iba, 2001), (Kremer and Binder, 1988) などが良いレビューを提供している．(Assaraf et al., 2000) が量子系のシミュレーションの文脈で．↩︎
diffusion Monte Carlo, projector Monte Carlo などと呼ばれる手法に等しい．Green’s function Monte Carlo もポテンシャル \(G\) の取り方が違うのみである (Assaraf et al., 2000)．また (Iba, 2001) 3.1節 p.282 にも言及がある．↩︎
(矢野浩一, 2014) p.190，(Ristic et al., 2004) 前文 p.xi．↩︎
(Chopin and Papaspiliopoulos, 2020, pp. 19.1節 p.371), (Creal, 2011, pp. 2.5.1節 p.258)．↩︎
そのアイデアは (Doucet et al., 2001, pp. 79–95) 第4章 から．↩︎
探知前追跡とは，信号が弱い，またはノイズが強い環境下において，信頼のおける初期信号を頼りにせずとも，物体追跡を実行するための手法．↩︎
(Ristic et al., 2004, p. 288) に示唆されている．↩︎
(Del Moral and Horton, 2023) は ensemble Kalman particle filtering methodology と呼んでいる．↩︎