粒子フィルタのソースを表示

{{翻訳中途|en:Particle Filter 15:26, 20 September 2007|date=2007年10月}}
'''粒子フィルタ'''（りゅうしフィルタ、{{lang-en-short|particle filter}}）や'''逐次[[モンテカルロ法]]'''（ちくじモンテカルロほう、{{lang-en-short|sequential Monte Carlo; SMC}}）とは、シミュレーションに基づく複雑なモデルの[[推定]]法である。[[1993年]][[1月]]に[[北川源四郎]]が'''モンテカルロフィルタ'''の名称で<ref name="Kitagawa1993">{{cite journal |last=Kitagawa |first= G. |date=1993-01 |title=A Monte Carlo filtering and smoothing method for non-Gaussian nonlinear state space models |url=https://www.ism.ac.jp/~kitagawa/1993_US-Japan.pdf |journal=Proceedings of the 2nd U.S.-Japan Joint Seminar on Statistical Time Series |pages=110-131 |accessdate=2019-11-20 }}</ref>、1993年[[4月]]にN.J. Gordonらが'''ブートストラップフィルタ'''の名称で<ref name="Gordon1993">{{cite journal |last= Gordon |first= N.J. |last2=Salmond |first2=D.J. |last3= Smith |first3= A.F.M. |date= 1993-04 |title= Novel approach to nonlinear/non-Gaussian Bayesian state estimation |journal=IEE Proceedings F - Radar and Signal Processing |volume=140 |issue=2 |pages= 107-113 |doi=10.1049/ip-f-2.1993.0015 }}</ref>それぞれ同時期に同様のものを発表した。

この手法はふつう[[ベイズ統計学|ベイズ]]モデルを推定するのに用いられ、バッチ処理である[[マルコフ連鎖モンテカルロ法]] (MCMC) の逐次 (オンライン) 版である。またこの手法は[[重点サンプリング]]法にも似たところがある。
うまく設計すると、粒子フィルタはMCMCよりも高速である。[[カルマンフィルター#拡張カルマンフィルター|拡張カルマンフィルタ]]や[[カルマンフィルター#Unscented カルマンフィルター|無香カルマンフィルタ]]（Unscented カルマンフィルタ）に比べて、サンプル点が十分多くなるとベイズ最適推定に近付くことからより高い精度の解が得られるので、これらの代わりに用いられることがある。また手法を組み合わせて、カルマンフィルタを粒子フィルタの提案分布として使うこともできる。

==目的==
粒子フィルタの目的は、観測不可能な状態列 <math>x_k</math> <math>(k=0,1,2,3, \ldots)</math> の[[確率分布]]を観測値 <math>y_k</math> <math>(k=0,1,2,3, \ldots)</math> から推定することである。ベイズ推定値 <math>x_k</math> は一つ過去の確率分布 <math>p(x_k|y_0,y_1,\ldots,y_k)</math> から得られるのに対し、マルコフ連鎖法では過去の全てを含む確率分布 <math>p(x_0,x_1,\ldots,x_k|y_0,y_1,\ldots,y_k)</math> より求められる.

==状態空間モデル==
粒子フィルタでは観測不可能な状態 <math>x_k</math> と観測値 <math>y_k</math> は以下のように表せるとする。

観測不可能な状態列 <math>x_0, x_1, \ldots</math> は以下を満たす1階[[マルコフ過程]]。つまり <math>x_k</math> は <math>x_{k-1}</math> で決まる。ただし初期値 <math>x_0</math> の確率分布 <math>p(x_0)</math> を持つものとする。なお、これは2つ前 <math>x_{k-2}</math> の状態を使えないという意味ではなく、必要なら2つ前の状態も <math>x_{k-1}</math> に含めて使うという意味である。
:<math>x_k|x_{k-1} \sim p_{x_k|x_{k-1}}(x|x_{k-1})</math>

観測データ列 <math>y_0, y_1, \ldots</math> は <math>x_0, x_1, \ldots</math> が既知であるという条件の下で独立である。換言すると <math>y_k</math> は <math>x_k</math> のみに従属する。
:<math>y_k|x_k \sim p_{y|x}(y|x_k)</math>

その上で、下記に従う[[状態空間法|状態方程式]]を状態空間モデルと呼ぶ。<ref>{{Cite book |和書 |author=北川源四郎|authorlink=北川源四郎 |year=2005 |title=時系列解析入門 |publisher=岩波書店 |page=209 |isbn=4000054554 }}</ref><ref>{{Cite book |和書 |author=樋口知之 |year=2011 |title=予測にいかす統計モデリングの基礎―ベイズ統計入門から応用まで |publisher=講談社 |page=29 |isbn=4061557955 }}</ref>

:<math>x_k = f_k(x_{k-1}, v_k) \,</math>
:<math>y_k = h_k(x_k, w_k) \,</math>

ただし <math>v_k</math> も <math>w_k</math> も異なる <math>k</math> について互いに独立であり、ある確率分布に従うノイズで、<math>v_k</math> をシステムノイズ、<math>w_k</math> を観測ノイズと呼ぶ。また、<math>f_k(\cdot)</math> と <math>h_k(\cdot)</math> は既知の関数である。

[[カルマンフィルター]]の状態方程式は状態空間モデルの一種であり、<math>x_k</math> と <math>y_k</math> が実数の列ベクトル、関数 <math>f_k(\cdot)</math> と <math>h_k(\cdot)</math> が線形、<math>v_k</math> と <math>w_k</math> が[[多変量正規分布]]に従うとすると、下記形式で表現でき、
:<math>
\begin{align}
x_k &= F_k x_{k-1} + G_k v_k \\
y_k &= H_k x_k + w_k
\end{align}
</math>
<math>x_k</math> の確率分布は多変量正規分布になり、カルマンフィルタによってベイズ推定と厳密に一致する解が得られる。線形でない場合などは、カルマンフィルタは1次近似に過ぎない。粒子フィルタも[[モンテカルロ法]]による近似には変わりないが、十分な数の粒子があれば高い精度が得られる。

==モンテカルロ近似==
粒子法は他のサンプルリング法と同様に、フィルタリング確率分布<math>p(x_k|y_0,\dots,y_k)</math>を近似する点列を生成する。したがって <math>P</math> 個のサンプル点があれば、フィルタリング確率分布による期待値は

:<math>\int f(x_k)p(x_k|y_0,\dots,y_k)dx_k\approx\frac1P\sum_{L=1}^Pf(x_k^{(L)})</math>

によって近似される。そして通常のモンテカルロ法と同様に <math>f(\cdot)</math>を適切に調整することで, 近似の程度に応じた次数までの [[モーメント (確率論)|モーメント]] を得ることができる。  

== Sampling Importance Resampling (SIR) ==
''Sampling importance resampling (SIR)'' アルゴリズムは、モンテカルロフィルタ（[[北川源四郎]] 1993<ref name="Kitagawa1993"/>）やブートストラップフィルタ（Gordon et al. 1993<ref name="Gordon1993"/>）による本来の粒子フィルタであるが、よく使われる粒子フィルタアルゴリズムである。ここではフィルタリング確率分布 <math>p(x_k|y_0,\ldots,y_k)</math> を<math>P</math>個の重みつき粒子によって近似する。

: <math>\{(w^{(L)}_k,x^{(L)}_k)~:~L=1,\ldots,P\}</math>.

''重み'' <math>w^{(L)}_k</math> は以後の <math>\sum_{L=1}^P w^{(L)}_k = 1</math> である粒子の相対事後分布(密度)の近似となっており、
<math>\sum_{L=1}^P w^{(L)}_k = 1</math>を満たす。

SIRは[[重点サンプリング]] の逐次版 (つまり、再帰的) と言える。
重点サンプリングにあるように、関数 <math>f(\cdot)</math> の推定値は重みつき平均
:<math>    \int f(x_k) p(x_k|y_0,\dots,y_k) dx_k \approx \sum_{L=1}^P w^{(L)} f(x_k^{(L)}) </math>
で近似できる。

粒子の個数が有限である場合、このアルゴリズムの性能は提案分布<math>\pi(x_k|x_{0:k-1},y_{0:k})</math>の選択に依存する。
最適な提案分布は目的分布
:<math>    \pi(x_k|x_{0:k-1},y_{0:k}) = p(x_k|x_{k-1},y_{k}) </math>
である。

しかしながら事前遷移確率
:<math>    \pi(x_k|x_{0:k-1},y_{0:k}) = p(x_k|x_{k-1}) </math>
を提案分布として用いることが多い。
なぜならそこからは容易に粒子をサンプルすることができるし、その後に重みを計算することもできるからである。

提案分布として事前遷移確率を用いるSIRフィルタは一般にブートストラップフィルタ・コンデンセーションアルゴリズムとして知られている。

唯一つを除いた全ての重みがゼロに近付くこと ― アルゴリズムの縮退 ― を防ぐために再サンプルが行なわれる。
このアルゴリズムの性能は再サンプリング方式の選びかたにもかかっている。
北川源四郎 (1993<ref name="Kitagawa1993"/>) によって提案された[[層化抽出法]]は分散の意味で最適である。

SIR法の1ステップは以下の様になる。

:1) <math>L=1,\ldots,P</math> について, 提案分布から粒子をサンプルする

:: <math>
x^{(L)}_k \sim \pi(x_k|x^{(L)}_{0:k-1},y_{0:k})
</math>

:2) <math>L=1,\ldots,P</math> について、重みを更新し、正規化定数を得る。
::<math>
        \hat{w}^{(L)}_k = w^{(L)}_{k-1} \frac{p(y_k|x^{(L)}_k) p(x^{(L)}_k|x^{(L)}_{k-1})} {\pi(x_k^{(L)}|x^{(L)}_{0:k-1},y_{0:k})}
</math>
このとき<math>\pi(x_k^{(L)}|x^{(L)}_{0:k-1},y_{0:k}) = p(x^{(L)}_k|x^{(L)}_{k-1})</math>
ならば更新式は
::<math>
            \hat{w}^{(L)}_k = w^{(L)}_{k-1} p(y_k|x^{(L)}_k), 
</math>
のように簡単になる。
 
:3) <math>L=1,\ldots,P</math>について、正規化された重みを計算する。

:: <math>
w^{(L)}_k = \frac{\hat{w}^{(L)}_k}{\sum_{J=1}^P \hat{w}^{(J)}_k}
</math>

:4) 有効粒子数の推定量を計算する。

:: <math>
\hat{N}_\mathit{eff} = \frac{1}{\sum_{L=1}^P\left(w^{(L)}_k\right)^2}
</math>


:5) もし有効粒子数が与えられた閾値より少なかったら <math>\hat{N}_\mathit{eff} < N_{thr}</math> 再サンプルを実行する。

::a) 現在の粒子の集合から、重みに比例した確率で <math>P</math> 個の粒子を描く。現在の粒子の集合を新しい粒子の集合で置き換える。

::b) <math>L=1,\ldots,P</math> について、 <math>w^{(L)}_k = 1/P</math> とする。


''Sequential Importance Resampling'' 等の用語は SIR フィルターの意味で用いられることがある。

== 逐次的 Importance Sampling (SIS) ==
再サンプルが無い点を除いて SIR と同様である。

==直接法==
直接法は他の粒子フィルタに比べて簡単で、合成と[[棄却サンプリング|棄却]]を利用する。
<math>k</math> 番目の一つのサンプル <math>x</math> を <math>p_{x_k|y_{1:k}}(x|y_{1:k})</math> から生成するのに以下の手順を踏む。

:1) <math>n = 1</math>とする。

:2) <math>\{1,..., P\}</math>上の[[離散一様分布|一様分布]]から<math>L</math>を生成する

:3) テスト粒子 <math>\hat{x}</math> を確率分布 <math>p_{x_k|x_{k-1}}(x|x_{k-1|k-1}^{(L)})</math> から生成する

:4) <math>\hat{y}</math> の確率を <math>\hat{x}</math> を用いて <math>p_{y|x}(y_k|\hat{x})</math> から生成する。ただし、<math>y_k</math> は観測値である

:5) 別の数<math>u</math>を<math>[0, m_k]</math>上の一様分布からを生成する。ただし、<math>m_k = \sup_x p_{y|x}(y_k|x) </math>

:6) <math>u</math> と <math>\hat{y}</math> を比較

::6a) <math>u</math> が大きければ step 2 以降を繰り返す

::6b) <math>u</math> が小さかったら <math>\hat{x}</math> を <math>x{k|k}^{(p)}</math> として保存して <math>n</math> を1増やす

:7) <math>n > P</math> ならば終了

目的は<math>k</math>における P 個の粒子を、<math>k-1</math> だけから生成することである。 
これにはマルコフ方程式が <math>x_{k-1}</math> だけから <math>x_k</math> を生成できるように記述できなければならない。
このアルゴリズムは <math>k</math> における粒子を生成するために<math>k-1</math>における<math>P</math>個の粒子の合成を利用しており、<math>k</math>において<math>P</math>個の粒子が生成されるまでステップ 2--6 以降を繰り返す。

これは <math>x</math> を2次元配列とみなすと、より視覚的に理解できる。
一つの次元は <math>k</math> であり、もう一つの次元は粒子番号<math>L</math>である。
例えば <math>x(k,L)</math> は <math>k</math> における<math>L</math>番目の粒子であり、<math>x_k^{(L)}</math> と書ける
(上記のアルゴリズムの記述に用いた通り)。
ステップ 3 は <math>k-1</math> においてランダムに選ばれた粒子<math>x_{k-1}^{(L)}</math> から''潜在的な'' <math>x_k</math> を生成し、ステップ 6で棄却 / 受領の判定が行われる。言い替えれば、<math>x_k</math> の値はそれまでに生成された <math>x_{k-1}</math> を用いて生成される。

== その他の粒子フィルタ ==
* [[:en:Auxiliary particle filter|Auxiliary particle filter]]
* [[:en:Gaussian particle filter|Gaussian particle filter]]
* [[:en:Unscented particle filter|Unscented particle filter]]
* [[:en:Monte Carlo particle filter|Monte Carlo particle filter]]
* [[:en:Gauss-Hermite particle filter|Gauss-Hermite particle filter]]
* [[:en:Cost Reference particle filter|Cost Reference particle filter]]

== 脚注 ==
{{脚注ヘルプ}}
{{Reflist}}

==参考文献==
* ''モンテカルロフィルタおよび平滑化について'', by 北川源四郎. 統計数理, Vol.44, No.1, pp. 31--48, 1996
* ''Sequential Monte Carlo Methods in Practice'', by A Doucet, N de Freitas and N Gordon.  Springer, 2001. ISBN 978-0387951461
* ''On Sequential Monte Carlo Sampling Methods for Bayesian Filtering'', by A Doucet, C Andrieu and S. Godsill, Statistics and Computing, vol. 10, no. 3, pp. 197-208, 2000 [http://citeseer.ist.psu.edu/doucet00sequential.html  CiteSeer link]
*  ''Tutorial on Particle Filters for On-line Nonlinear/Non-Gaussian Bayesian Tracking (2001)''; S. Arulampalam, S. Maskell, N. Gordon and T. Clapp; [http://citeseer.ist.psu.edu/maskell01tutorial.html CiteSeer link]
* ''Particle Methods for Change Detection, System Identification, and Control,'' by Andrieu C., Doucet A., Singh S.S., and Tadic V.B., Proceedings of the IEEE, Vol 92, No 3, March 2004. [http://www-sigproc.eng.cam.ac.uk/%7Esss40/papers/andrieu04_pfForChangeDetectionSysIdControl.pdf SMC Homepage link]
* ''A Tutorial on Bayesian Estimation and Tracking Techniques Applicable to Nonlinear and Non-Gaussian Processes'', A.J. Haug, The MITRE Corporation; [http://www.mitre.org/work/tech_papers/tech_papers_05/05_0211/index.html Mitre link]
* ''Distributed Implementations of Particle Filters'', Anwer Bashi, Vesselin Jilkov, Xiao Rong Li, Huimin Chen, [http://ece.engr.uno.edu/isl/Reprints/C121.pdf Available from the University of New Orleans]
* ''A survey of convergence results on particle filtering methods for practitioners'', Crisan, D. and Doucet, A., IEEE Transactions on Signal Processing 50(2002)736-746   [https://doi.org/10.1109/78.984773 doi:10.1109/78.984773]
* ''Beyond the Kalman Filter: Particle Filters for Tracking Applications'', by B Ristic, S Arulampalam, N Gordon. Artech House, 2004. ISBN 1-58053-631-X.

== 関連項目 ==
* [[アンサンブルカルマンフィルタ]]
* [[カルマンフィルタ]]、ガウス分布の解析的な推定器
* [[反復ベイズ推定]]
* [[北川源四郎]]

==外部リンク==
* [http://www-sigproc.eng.cam.ac.uk/smc/ Sequential Monte Carlo Methods (Particle Filtering)] homepage on University of Cambridge
* [http://www.cs.washington.edu/ai/Mobile_Robotics/mcl/ Dieter Fox's MCL Animations]
* [http://www.cs.ubc.ca/~nando/software.html Nando de Freitas' free software]
* [http://web.engr.oregonstate.edu/~hess/ Rob Hess' free software]

{{確率分布の一覧}}
{{DEFAULTSORT:りゆうしふいるた}}

[[Category:モンテカルロ法]]
[[Category:アルゴリズム]]
[[Category:数値解析]]
[[Category:計算科学]]
[[Category:ロボット制御]]
[[Category:統計学]]
[[Category:数学に関する記事]]