第 11 回横幹連合コンファレンス
2020.10.8-9 統計数理研究所

粒子フィルタ〜線形ガウスの枠を超えた汎用な状態推定法
生駒 哲一 （日本工業大学）
Particle Filter – Universal State Estimation Method Beyond Linear-Gaussian Restrictions

Norikazu Ikoma (Nippon Institute of Technology)
Abstract– State space representation of dynamical system is widely used in prediction and control, etc., andstate estimation is a fundamental problem to be solved in this formulation. While analytical state estimations
via Kalman filter, etc., are available for systems consisting of linear equations and Gaussian distributions,
there is no analytical solution to the state estimation problem for nonlinear and/or non-Gaussian state spacemodels in general case. As a breakthrough to this issue, in early 1990s, a particle filter called Monte
Carlo Filter has been proposed in Japan as the first universal approximation method of state estimation
for nonlinear and/or non-Gaussian models by utilizing many realizations in the state space to represent
probability distribution of posterior state. Due to the universal property and allowed flexibility in modeling,
now, particle filters have become standard methods in many fields, such as natural science, social science,engineering, and so on.
Index terms– particle filter, nonlinear/non-Gaussian state space model, state estimation in dynamicalsystem1

推薦対象および推薦理由の概要「粒子フィルタ」をコトつくり至宝に相応しい対象として推薦し，その理由等について論説する．その概略は，次の通りである．動的システムの状態推定という汎用的な問題設定において，既存のカルマンフィルタ
では線形性とガウス性の二つの制約があったが，これら制約を超える方法論として「粒子フィルタ」が提案された．汎用的な問題設定である事から，さまざまな分
野において定番の手法となり，各分野において，データ分析などの発展に大きく寄与した．このことから，「粒子フィルタ」は，コトつくりの至宝として十分な価値を持つと考えられる．本稿ではこれを，問題設定とその提案の経緯，発展と応用を示しつつ，論説する．
予測や制御などで広く用いられる動的システムの状
態空間表現にて，モデルが線形式とガウス分布から成る場合はカルマンフィルタ等での解析的な状態推定ができる一方，非線形非ガウスのモデルに対する汎用な近似解法として，状態空間中の多数の実現値で状態推定の確率分布を表現する「粒子フィルタ」（Particle Filter）
が 1990 年代前半に日本で提案された．その高い汎用性から，自然科学，社会科学，工学など多くの分野での標準的な手法になっている．
本稿の構成は，次の通りである．まず問題設定として，一般状態空間モデルと状態推定の問題の定式化と，その形式的な解を示した後，粒子フィルタが提案されるに至った経緯を説明する．そこでは，まず，カルマン
フィルタにより解析的に解くことのできる線形ガウス
状態空間モデルと，カルマンフィルタの時間更新の式を示した後，非線形・非ガウスの状態空間モデルと，その近似解を得るいくつかのアプローチを紹介する．粒子群による近似のアプローチをとる「粒子フィルタ」として最初に提案された「モンテカルロフィルタ」のア
ルゴリズムを，確率分布表記および方程式表記の二つについて示す．
モンテカルロフィルタに端を発する，様々な粒子フィ
ルタ達への発展として，動画像追跡に特化した粒子フィルタである「CONDENSATION」による「検出前追跡」の具現化，初期の工夫である補助変数粒子フィルタ（APF:Auxiliary Particle Filter），理論的・方法論
的な成熟形である「逐次モンテカルロフィルタ」，複数対象の同時追跡における PHD フィルタとその各種実装，可能性と必然性を扱う粒子フィルタ，平滑化の推定と固定パラメータ推定について，それぞれ論説する．
多彩な応用課題での活躍として，レーダ観測等におけるターゲット追跡，経済学や社会科学などの時系列
解析，動画像中の対象物の追跡，自然科学においてシ
ミュレーションにデータを取り込む「データ同化」，ロボティクス分野における自律移動ロボットの自己位置推定や地図学習，センサ融合について，概観する．
本稿の末尾にて，「まとめ」として，粒子フィルタの果
たしてきた役割を概観した後，「先導力」「意味力」「規範力」「解決力」の観点から，「粒子フィルタ」がコトつ
くり至宝に相応しい対象である理由を述べる．

2

一般状態空間モデルと状態推定

動的システムの状態空間表現である「状態空間モデ
ル」を，ここでは一般的な式として，確率的な表記で示す．離散時刻を整数 k で表す．直接には観測できな
い，「隠れ状態」xk に対し，その時間変化を司る条件付
き確率分布
f (xk |xk−1 ; Θf )(1)
を「システムモデル」とし，同一時刻内で隠れ状態か
ら「観測」yk が派生する条件付き確率分布

h(yk |xk ; Θh )(2)

を「観測モデル」とする．初期時刻 k = 0 においては，隠れ状態の「初期確率分布」として次を想定する

p(x0 ; Θ0 ).

(3)

上記のそれぞれの確率分布を規定するパラメータをまとめて，「固定パラメータ」とし，次で記述する

Θ = {Θ0 , Θf , Θh } .

(4)

これらの固定パラメータは，以下にて，特に必要のない場合には，記述を省略することにする．

「状態推定」とは，観測系列

y1:k ≡ {y1 , y2 , · · · , yk }(5)が与えられた下での，隠れ状態の条件付き確率（事後確率）分布を求める課題である．これを，時刻 k = 1, 2, · · ·の順に効率よく求める手順として，式 (3) の「初期確率分布」から開始し，「１期先予測」Zp(xk |y1:k−1 ) = f (xk |xk−1 )p(xk−1 |y1:k−1 )dxk−1(6)と「ろ波（フィルタ）」p(xk |y1:k ) =h(yk |xk )p(xk |y1:k−1 )p(yk |y1:k−1 )ΘNYp(yk |y1:k−1 ; Θ)

(8)

k=1

と，ベイズ推定にて隠れ状態との同時事後確率分布p(x0:N , Θ|y1:N ) ∝ p(y1:N |x0:N ; Θ)p(x0:N , Θ)

(9)を求めるアプローチの，二つがある．後者では，Θ について周辺化した事後確率分布について，それを最大化する MAP 推定値

Θ̂MAP = arg max p(Θ|y1:N )

(10)

Θ

などが用いられる．これらの式は，形式的な求解式であり，それに基づき，具体的な求解アルゴリズムを導出する必要がある．

3粒子フィルタへの道のり

xk = Fk xk−1 + Gk vk , vk ∼ N (0, Qk )状態推定課題の具体的な求解アルゴリズムとして，線形ガウス状態空間モデルの場合には，解析的な解である「カルマンフィルタ」がある．一方，非線形，あるいは非ガウスの状態空間モデルの場合には，一般的な求解アルゴリズムは存在しない為，解である事後確率分布に対して何等かの近似を行った上で，一般的な解を求める『近似アルゴリズム』が複数提案されている．
粒子フィルタは，それら近似アルゴリズムの中でも，汎用性の高い方法である．
以下では，まず，線形ガウス状態空間モデルとカル
マンフィルタについて説明した後，非線形非ガウス状
態空間モデル（方程式による表記）と，粒子フィルタの最初の提案として知られている「モンテカルロフィ
ルタ」について説明する．

(11)

で記述される．vk を「システムノイズ」と呼ぶ．また，
式 (2) の観測モデルは，線形ガウスの「観測方程式」

yk = Hk xk + wk , wk ∼ N (0, Rk )(7)

を，時刻を進めながら交互に適用する方法がある．
これら二式 (6)(7) は，状態推定課題の「形式的な解」
であり，解を明確に求める具体的なアルゴリズムではない．具体的な求解アルゴリズムは，これら二式や，あるいは方程式に基づく議論などから，別途，導出する必要がある．
固定パラメータは，モデルを規定する値として所与となる部分もあるが，観測データからの推定を必要とする部分もある．固定パラメータの推定は，固定長の観測系列 y1:N に対して，最尤法によるアプローチ

Θ̂ML = arg max

3.1 カルマンフィルタ
隠れ状態がおよび観測が，それぞれ，実数ベクトル
xk ∈ ℜn および yk ∈ ℜm の場合，式 (1) のシステムモ
デルと式 (2) の観測モデルは，代数的な方程式で記述できる．方程式が線形で，確率項がガウス分布（正規分布）に従う場合，式 (1) のシステムモデルは，線形ガウスの「システム方程式」(12)で記述される．wk を「観測ノイズ」と呼ぶ．更に，式(3) の初期確率分布に，ガウス分布p(x0 ) : N (x̂0|0 , V0|0 )(13)を想定する．式 (11) (12) (13) にて，「線形ガウス状態
空間モデル」が構成される．
線形ガウス状態空間モデルにおいて，システムノイ
ズ，観測ノイズ，初期確率分布に対する適切な独立性の前提の下で，状態推定で求めるべき事後確率分布は，正規分布に従うことが導かれる．すなわち，１期先予
測分布についてp(xk |y1:k−1 ) : N (x̂k|k−1 , Vk|k−1 )(14)と表され，ろ波（フィルタ）分布について，

p(xk |y1:k ) : N (x̂k|k , Vk|k )(15)と表される．式 (11) 〜(15) にて，N (x̂, V) は，平均ベ
クトル x̂，分散共分散行列 V の正規分布を示す．
正規分布（ガウス分布）は，平均と分散の二つの要素
で完全に規定される為，状態推定のアルゴリズムとして，平均と分散を時間更新する手順が導出できる．これが「カルマンフィルタ」1)2) である．そこでは，１期先予測として，½
x̂k|k−1 = Fk x̂k−1|k−1(16)
Vk|k−1 = Fk Vk−1|k−1 F′ k + Gk Qk G′ kと更新し，ろ波（フィルタ）として，¡¢½
x̂k|k = x̂k|k−1 + Kk yk − Hk x̂k|k−1
Vk|k = [I − Kk Hk ] Vk|k−1(17)と更新する．行列 Kk は「カルマンゲイン」と呼ばれ£¤−1
Kk = Vk|k−1 H′ k Hk Vk|k−1 H′ k + Rk
(18)である．式 (16) (18) にて，H′ は行列 H の転置を表す．

3.2 モンテカルロフィルタ
システム方程式や観測方程式に非線形要素が存在す
る場合，あるいは，システムノイズや観測ノイズが非ガウス分布に従う場合，下記の「非線形非ガウス状態
空間モデル」となる
xk = fk (xk−1 , vk ) , vk ∼ qk (v),(19)

yk = hk (xk , wk ) , wk ∼ rk (w).

(20)

粒子フィルタの導出に際し，wk について陽に解いて

wk = h̃k (yk , xk )(21)と記述が可能な，関数 h̃k の存在を前提する．
非線形非ガウス状態空間モデルでは，事後確率分布のベイズ推測に対し，共役性（共役事前分布などとも呼ばれる）が一般には成り立たない．そのため，パラメトリックな事前確率分布から推定を開始しても，パラメトリックな事後確率分布が得られるとは限らない．
勿論，事前分布と同じパラメトリックな確率分布族に，事後分布が含まれるという保証も，一般には得られな
い．一方，カルマンフィルタの場合には，確率項が全てガウス分布であり，線形計算によるガウス分布の再生性から，共役性が成り立つため，解析的なアルゴリ
ズムが導出可能であった．
解析的な状態推定のアルゴリズム導出は，非線形非ガウス状態空間モデルでは難しい為，近似的なアプロー
チを考えることになる．そこでは，推定すべき事後確率分布を，どう近似するかが，大きなポイントである．
混合ガウス分布による近似とカルマンフィルタに基づくアルゴリズム導出 3) や，関数形状の数値的な近似と数値積分によるアプローチ 4) などが提案される中，確率分布に従う多数の実現値（粒子，パーティクル）による近似が考案され，そこから粒子フィルタのアルゴ
リズム 5) 7) が導出された．これが，粒子フィルタとして最初に提案された「モンテカルロフィルタ」である．
なお，同時期に「ブートストラップフィルタ」6) という名称で，同一のアルゴリズムが提案されており，世界的にはこの名称が広く通用しているようである．
モンテカルロフィルタでは，初期状態の確率分布を，

oM
n
(i)
p(x0 ) ≈ x̂0|0

(22)

i=1

と，粒子群で近似する．同様な近似を，１期先予測

oM
n
(i)
p(xk |y1:k−1 ) ≈ x̂k|k−1

i=1

(23)

と，ろ波（フィルタ）

oM
n
(i)
p(xk |y1:k ) ≈ x̂k|k

i=1

(24)

についても行う．
確率分布に対する解析的な更新ではなく，代わりに，粒子群に対する更新を行うことで，状態推定を近似的に進めることができる．これを行うのが，「モンテカルロフィルタ」のアルゴリズムである．そこでは，まず，初期分布に従う粒子群を生成する(i)x̂0|0 ∼ p(x0 ) , i = 1, 2, · · · , M.

(25)

次に，時刻 k を適切に進めながら，時間更新の手順を適用する．手順は３つから成り，一つ目は「１期先予測の粒子群生成」(i)(i)

x̂k|k−1 ∼ f (xk |x̂k−1|k−1 ) , i = 1, 2, · · · , M(26)

Algorithm 1 モンテカルロフィルタ：確率分布表記
Input: y1:N
oM ¾N
oM n
oM ½n
n
(i)
(i)
(i)
, x̂k|kx̂k|k−1,
Output: x̂0|0
i=1

i=1

i=1

Initialization :1: for i = 1 to M do(i)2:x̂0|0 ∼ p(x0 )3: end for
Loop of Time Update :4: for k = 1 to N do
5:
for i = 1 to M do
(i)
(i)
6:
x̂k|k−1 ∼ f (xk |x̂k−1|k−1 )
7:

8:
9:
10:

(i)

k=1

(i)

ωk ∝ h(yk |x̂k|k−1 )end forfor i = 1 n
to M do
oM
(j)
(i)
(j)
x̂k|k ∼ x̂k|k−1 with prob. ωk

j=1

end for
12: end for
oM ½n
oM n
oM ¾N
n
(i)
(i)
(i)
13: returnx̂0|0,
x̂k|k−1, x̂k|k
11:

i=1

i=1

i=1

k=1

である．これを方程式で表記すると，下２式となる．
(i)

(i)

(i)

(27)

vk ∼ qk (v) , i = 1, 2, · · · , M(28)x̂k|k−1 = fk (x̂k−1|k−1 , vk ) , i = 1, 2, · · · , M(i)

手順の二つ目は，観測モデルに基づき，生成された１期先予測粒子について「尤度評価」を
(i)

(i)

ωk ∝ h(yk |x̂k|k−1 ) , i = 1, 2, · · · , M(29)と行う．これを方程式で表記すると，下２式となる．
(i)

(i)

ek = h̃k (yk , x̂k|k−1 ) , i = 1, 2, · · · , M
(i)

(i)

ωk ∝ rk (ek ) , i = 1, 2, · · · , M

(30)
(31)

式 (29) (31) にて，比例記号 ∝ を用いたのは，右辺の定数値の計算が省略可能である事と，ω を後で確率値として使うために正規化を行い，総和を１にする
PM (i)（ i=1 ωk = 1）事が，その理由である．
手順の三つ目は，「再サンプリング」と呼ばれる，復元抽出の手順を，１期先予測の粒子群に対して実施する．

oM
n
(j)
(j)
(i)
, i = 1, 2, · · · , Mx̂k|k ∼ x̂k|k−1 with prob. ωk
j=1

(32)
すなわち，手順の二つ目で評価した尤度の値に対し，そ(j)
れに比例する確率 ωk で，復元抽出が行われる．
上で説明した時間更新の手順を，アルゴリズムとして明記すると，確率分布による表記では Algorithm 1，
方程式による表記では Algorithm 2 となる．アルゴ
リズム中では，状態空間モデル（システムモデルおよび観測モデル）が素直に適用され，カルマンフィルタにあるような逆行列の計算は行わない．多数の粒子を用いるものの，その計算量は，再サンプリングの箇所を除き，粒子数 M のオーダーである．再サンプリン
グの箇所は，工夫により，M log M やそれ以下のオー

Algorithm 2 モンテカルロフィルタ：方程式表記
Input: y1:N
oM ¾N
oM n
oM ½n
n
(i)
(i)
(i)
, x̂k|kx̂k|k−1,
Output: x̂0|0
i=1

i=1

i=1

Initialization :1: for i = 1 to M do(i)2:x̂0|0 ∼ p(x0 )3: end for
Loop of Time Update :4: for k = 1 to N do
5:
for i = 1 to M do(i)
6:vk ∼ qk (v)
(i)
(i)
(i)
7:
x̂k|k−1 = fk (x̂k−1|k−1 , vk )

8:
9:
10:
11:
12:

k=1

(i)

(i)

ek = h̃k (yk , x̂k|k−1 )
(i)

(i)

ωk ∝ rk (ek )end forfor i = 1 n
to M do
oM
(j)
(j)
(i)
x̂k|k ∼ x̂k|k−1 with prob. ωk

j=1

end for
14: end for
oM ¾N
oM n
oM ½n
n
(i)
(i)
(i)
, x̂k|k,x̂k|k−1
15: return
x̂0|0
13:

i=1

i=1

i=1

k=1

ダーに抑えることが可能である．また，粒子ごとに独立に処理が可能な部分が多く，並列計算による実装で効率を上げる事ができる．再サンプリングでの復元抽出のための確率分布（式 (32) の右辺）を準備する箇所
を除き，粒子ごとに並列化できる．
粒子フィルタでは，状態 xk および観測 yk は，必ずしも実数ベクトル（xk ∈ ℜn や yk ∈ ℜm ）でなくてもよい．ガウス性を仮定しない為，ベクトルには，離散値を取る要素を含むこともできる．なお，状態ベクト
ルの要素が全て離散値を取る場合は，「隠れマルコフモ
デル」に相当する．更には，数値でない要素も可能であり，リスト構造など構造を持ったデータや，テキス
ト（文章）などを扱うこともできる．

4

様々な粒子フィルタ達への進化「モンテカルロフィルタ」により創始された粒子フィ
ルタは，その後，様々な発展的あるいは派生的な粒子フィルタ達 13) へと進化し，より先進的な研究が現在でも展開されている 29) ．粒子フィルタは，理論的あるいは数学的な基盤が堅固であるため，このような発展性を持つことができたと考えられる．以下では，そのような発展的あるいは派生的な粒子フィルタ達について，重要と思われるものをいくつか紹介する．

4.1

検出前追跡：TBD(Track-Before-Detect)

粒子フィルタを用いて，動画像中の特定対象物の追跡を行う方法として，CONDENSATION(CONditional
DENSity propagATION) 8) がある．そこで初めて提案された，画像中の対象物らしさの評価方法は，旧来の状態推定によるターゲット追跡とは異なるアプローチ
である「検出 “前” 追跡」（TBD:Track-Before-Detect）
を，粒子フィルタにより具現化したものである．
従来のターゲット追跡は，粒子フィルタでは，まず(i)
式 (30) で観測の予測誤差 ek を算出し，次に式 (31)で予測誤差の尤度評価を行う．カルマンフィルタでは，
式 (17) 中の yk − Hk x̂k|k−1 が予測誤差の算出であり，その一部分の Hk x̂k|k−1 が，予測の平均値計算である．
一般には「観測のデータ値 yk が所与」は自明なのだ
が，動画像追跡では「画像」は多次元情報であり，これを観測のデータ値 yk とするのは得策ではない．予測値や予測誤差の算出，尤度評価において，計算コストが高くなり，現実的ではないからである．
現実的な方法として，画像において検出処理を行い，追跡対象の画像上での位置座標を算出して，これを観測のデータ値 yk として用いるアプローチがある．まず検出処理を行い，その後，状態推定で追跡を行うので，「検出 “後” 追跡」と呼ぶことができる．このアプロー
チでは，動画像追跡は，旧来のターゲット追跡と同等の課題に還元される．しかし，検出処理における誤検
出や未検出などを孕み，より複雑な状態空間モデリン
グが必要となる（これについては，4.4 で再び触れる）．
画像に対して検出処理をせず，状態推定による追跡を画像に対して直接行い，結果的に検出が遂行されるアプローチが「検出 “前” 追跡」である．CONDENSATION(CONditional DENSity propagATION) 8) が具現化したのは，追跡対象の典型形状を事前に登録しておき，それをアフィン変換して画像にあてはめ，適合度を算出する方法である．画像にあてはめた形状のエッヂ情報に基づき，尤度の評価を行う．粒子はアフィン変換パラメータであり，追跡対象の位置，スケール，向きの情報を持つ．追跡対象の部分的な隠れも考慮し，頑健な尤度評価の工夫もなされている．
より一般には，粒子フィルタでの尤度評価は，式 (29)に従い，観測の予測値や予測誤差を経由せずに，直接に行ってよい．例えば，数式での陽な表現をせずとも，尤
度値の評価を行うアルゴリズム等が明確であれば，それを直接に適用してよい．すなわち，粒子フィルタは，誤差等の確率分布に非ガウス分布を想定でき，非線形な演算を含む複雑な式を用いることができるに留まらず，更に自由なモデリングを許すものとなっている．そこでは，観測における予測とその誤差評価が不要な為，式 (21) で前提とした関数 h̃k の存在は問わない．

4.2

補助変数粒子フィルタ：APF

粒子フィルタでは，１期先予測を行った後に尤度の評価を行うため，予測が大幅に外れた場合には，状態推定の近似アルゴリズムが破綻する．これを避けるた
めの，初期の頃の改良の一つとして，補助変数粒子フィルタ（APF:Auxiliary Particle Filter）10) がある．
上述したアルゴリズムの破綻を避けるアプローチに，「重点サンプリング」がある．そこでは，現在の観測値
yk を考慮した「プロポーザル分布」から有望な粒子を生成して用い，１期先予測とは異なる分布を用いた事で生じる偏りを「インポータンス重み」により是正し，所望の確率分布に対する近似を担保する．
これを粒子フィルタに素朴に適用すると，インポー
タンス重みの算出において，１期先予測の確率分布の評価が必要となり，それは全粒子に対する計算を伴う．
その評価を，全ての粒子に対し行う為，計算のコスト
が粒子数 M の二乗となってしまい，実用的ではない．
この問題点を克服するために，補助変数を導入して，インポータンス重みの算出を容易にした効果的なアルゴリズムが，補助変数粒子フィルタ（APF:Auxiliary
Particle Filter）である．そこでの補助変数は，１時刻(i)
前のろ波粒子 x̂k−1|k−1 の粒子番号 i である．推定すべき『ろ波』の確率分布 p(xk |y1:k ) を，状態 xk と粒子番
号 i との同時確率分布 p(xk , i |y1:k ) へと拡張し，性質の良いプロポーザル分布 q(xk , i |y1:k ) を採用して，重点サンプリングを定式化している．粒子フィルタのア
ルゴリズムに従って生成された粒子は，状態 xk と粒子
番号 i の対となるが，粒子番号は捨て去り，状態のみを扱うだけで，所望の結果が得られる．

4.3 逐次モンテカルロフィルタ
粒子フィルタの改良として，より包括的なアプロー
チが 1990 年代の後半に提案され 11) 12) ，成熟形に到達した．これを「逐次モンテカルロフィルタ」と呼ぼう．
そこでは，状態系列の事後確率分布の更新式

h(yk |xk )f (xk |xk−1 )p(yk |y1:k−1 )(33)に基づき，１期先予測を経由せずに，時間更新を行う．
粒子の時間系列に重みを付与した，重み付き粒子系列群
n³
´oM
(i)
(i)
x0:k , ωk
により，式 (33) 左辺の同時事後確
p(x0:k |y1:k ) = p(x0:k−1 |y1:k−1 )i=1率分布を近似する．重点サンプリングの考え方に従い，
´oM
n³
(i)
(i)
１時刻前の重み付き粒子系列群
x0:k−1 , ωk−1
i=1

に対し，新時刻の粒子をプロポーザル分布から生成
(i)

(i)

xk ∼ q(xk |x0:k−1 , y1:k )(34)

して粒子の時間系列 x0:k−1 に追加し，追加に伴うインポータンス重みの更新を
(i)

(i)

(i)

(i)

ωk = ωk−1

(i)

(i)

h(yk |xk )f (xk |xk−1 )
(i)

(i)

q(xk |x0:k−1 , y1:k )(35)と行う．なお，この更新式は，式 (33) から導かれる．
重み付き粒子系列群において，重みがほぼ均等である場合は，上述した２つの手続き（式 (34) と式 (35)）
を，時間を進めながら適用してゆくことで，同時事後確率分布の時間更新が近似的に為される．だが，これを進めてゆくと，重みの値が粒子系列間で偏り，ごく少数の粒子系列だけが有意な重みの値を持ち，他の粒
子系列は重みの値が小さくなってしまう「重みの縮退」が生じる．これを回避するには，再サンプリングの手続きを適用すればよい．重みが偏っているかどうかを判定し，偏っている場合にのみ再サンプリングの手続き実施することで，無用なモンテカルロ誤差を避ける．
再サンプリング実施後は，重みを均等化する．
マ ル コ フ 連 鎖 モ ン テ カ ル ロ 法（MCMC:Markov
Chain Monte Carlo）は，粒子フィルタ発明の前から知られている，確率分布に従う粒子を生成する常套的な方法である．これを，粒子フィルタのアルゴリズム中に組み入れ，粒子の配置が所望の確率分布により近くなるよう改良する方法も提案されている．ただし，マル
コフ連鎖モンテカルロ法では，粒子の生成を反復的に実施する為，実時間性を要する課題には適していない．
「ラオ・ブラックウェル化」と呼ばれる，粒子フィル
タと，カルマンフィルタ等の解析的フィルタとの，ハ

イブリッドなアルゴリズムも提案されている．そこでは，事後確率分布の推定に，状態のある一部分に対して粒子フィルタの手続きを適用し，残りの部分については粒子所与の下で解析的フィルタの手続きを適用す
る．状態空間モデルで記述した際に，解析的フィルタ
が部分的に適用できる場合には，ラオ・ブラックウェ
ル化を行うことで，状態の要素全てを粒子フィルタで推定した場合と比べて精度のよい推定が可能となる．

4.4 PHD フィルタ：複数対象の同時追跡レーダ等で観測した情報から，飛翔体や船舶などの位置や速度等を精度よく推定する課題が「ターゲット
追跡」である．ターゲット追跡は，カルマンフィルタが発明される前からの古典的な課題であるが，そこに内包される諸問題は奥が深く，現在でも様々な研究が活発に進められている．古典的な文脈においても，レー
ダによる極座標系での観測や，電波の性質を考慮した観測過程のモデル化など，非線形あるいは非ガウスの状態空間モデルが必要となる場面は多い．
レーダで取得した原信号に対して，追跡対象の検出処理を行うと，「誤検出」や「未検出」などが発生し得るため，それらを含めた複雑な状態空間モデリングが必要になる．これは，動画像追跡において検出処理を伴う場合にも生じる事項である．現実的なシーンでは，追跡対象は複数存在し，観測領域への侵入等に伴う「出現」や「消滅」があり得るので，対象の個数は一定ではない．従って，これらの可変個数の複数対象を，区別しながら，同時に追跡することが必要となる．
複数対象の状態や観測を，確率的な要素を持つ「ランダム有限集合」で表し，状態空間モデルとして定式化することができる．状態推定は非常に複雑となり現実的ではないが，その１次モーメント（すなわち単一対象の空間）である PHD(Probability Hypothesis Density：
確率仮説密度）により推定を行う「PHD フィルタ」が，解析的な更新式として提案されていた 14) ．これに基づく具体的なアルゴリズムとして，PHD を粒子群で近似する逐次モンテカルロ実装が提案された 15) 16) ．
その後，これが更に発展してゆき，線形ガウスモデ
ルに対する解析的なフィルタアルゴリズムとして，ガ
ウス混合分布による PHD フィルタである「GM-PHD
フィルタ」19) ，対象個数（Cardinality）を明示的に扱う「C-PHD フィルタ」20) ，ガウス混合分布で対象個数を明示的に扱う PHD フィルタである「GM-CPHD
フィルタ」21) などが提案された．より最近では，対象ごとに固有のラベルを付与し，ポアソン分布ではなくベルヌーイ分布に基づくより広い定式化の「ラベル付き複数ベルヌーイフィルタ」28) などが提案され，今も活発に研究が進められている．

4.5 可能性必然性フィルタ
一般状態空間モデルは確率モデルの一つであり，確率論に基づき論理が展開されている．確率の３つの公理のうち，「加法性」（完全加法性）を，「単調性」に緩めることで，可能性と必然性を表す不確実性測度論が展開可能である．このような，数学的制約を緩めた測度は，例えば人の選好行動のような，必ずしも合理的とは言えない状況のモデル化に貢献できる．そこでは，確率論のように全てを見通すモデル化ではなく，無知である事を知識として記述することができる．

可能性と必然性を表す不確実性測度に基づく状態空間モデリングと，一般化されたベイズ推測に基づき，粒子群の集合を単位（すなわち粒子）とした粒子フィル
タのアルゴリズムが提案されている 24) ．推定対象は，離散的な変数のみに限られるが，例えば SLAM 問題におけるグリッドマップの占有判定などで，ラオ・ブラッ
クウェル化にて解析的フィルタを併用し，有益な結果が得られることが報告されている．

4.6 平滑化の推定
固定区間のデータ y1:N に対し，時刻 0 ≤ k ≤ N の状態 xk について周辺事後確率分布 p(xk |y1:N ) を求めることを「固定区間平滑化」といい，その形式的解はp(xk |y1:NZ ) = p(xk |y1:k )×p(xk+1 |y1:N )f (xk+1 |xk )
dxk+1
p(xk+1 |y1:k )(36)

である．この式に基づき，平滑化のアルゴリズムを導くことができる．そこでは，予め，１期先予測とろ波の推定を行ってその結果を記憶しておき，それらを参照しながら，時間を遡って平滑化の推定を順に求める．
線形ガウス状態空間モデルの場合には，平滑化の確率分布はガウス分布となる為，カルマンフィルタと同様に，平滑化分布の平均と分散を時間更新する解析的なアルゴリズムが導出されている．一方，非線形非ガ
ウス状態空間モデルの場合には，解析的な一般アルゴ
リズムは導出できず，粒子群による平滑化分布の近似が有用となり，様々なアルゴリズムが提案されている．
4.3 の逐次モンテカルロフィルタでは，粒子系列を保持し，過去の状態も同時に推定している．しかし，再サンプリングの手続きを実施する度に，系列中の粒子の多様性が減少し，それは過去に遡るほど顕著になる．
過去の時刻の周辺分布として，現在時刻から一定時間
L > 0 だけ遡った p(xk−L |y1:k ) を「固定ラグ平滑化」というが，このラグ L が小さい範囲においては，逐次
モンテカルロフィルタは有効である．

4.7 固定パラメータの推定
状態空間モデルの固定パラメータ（式 (4)）は，状態遷移のダイナミクスを規定するもの（線形ガウスの場
合は状態遷移行列 Fk ），観測過程を規定するもの（線形ガウスの場合は観測行列 Hk ），初期分布やシステム
ノイズ，観測ノイズを規定するもの（線形ガウスの場
合は初期分布の平均 x̂0|0 と分散 V0|0 ，システムノイズの係数行列 Gk ，システムノイズの分散行列 Qk ，観測ノイズの分散行列 Rk ）など，各種の役割を担う，固定的な要素である．典型的なモデリングにおいては，固定パラメータのうちの一部を未知とし，データからの推測でそれら未知要素の値を求める．
最尤法に基づくアプローチ（式 (8)）は，カルマン
フィルタでは厳密な尤度関数値が算出されるので，これを目的関数として数値的最適化の手法が適用可能である．一方，粒子フィルタは，乱数を用いたアルゴリ
ズムである為，そこで得られる尤度関数の近似値はモンテカルロ誤差を含み，数値的最適化の手法を直接適用すると困難が生じる．その解決策として，尤度関数の勾配ベクトル（スコア統計量）やヘッセ行列（フィッシャー情報行列）を粒子群から適切に求める工夫 23) などが提案されている．

固定パラメータと状態との同時確率分布をベイズ推定するアプローチ（式 (9)）では，初期のものとして，状態を拡張して固定パラメータを組み入れ，固定パラメータに対しても微小な変化を許す恣意的なダイナミクスを導入した，自己組織型状態空間モデリング 9) と，その各種の改良などがある．そこでは，状態推定を通して，固定パラメータの値が求められる．
ベイズ推定のアプローチのうち，より発展的で実用性の高いものとして，マルコフ連鎖モンテカルロ法の枠組みに従い，粒子フィルタに基づき工夫して構築された状態系列のプロポーザル分布から固定パラメータ
所与下での状態系列の実現値を生成（更新）し，固定
パラメータのプロポーザル分布と共にギブス・サンプ
ラを構成する方法 22) などが提案されている．

5

多彩な応用課題での活躍

一般状態空間モデルは広いクラスの動的システムを記述でき，粒子フィルタが基づく状態推定により，直接には観測できない隠れ状態の確率分布を求めることができる．粒子フィルタは，高い汎用性を持つことから，様々な分野において，多彩な応用課題に適用されている．以下では，その典型例をいくつか概観する．

5.1 ターゲット追跡
古典的な課題であるが，有意なノイズを伴う原信号からの追跡対象の位置や速度等の推定は，ウイナーフィ
ルタに端を発し，カルマンフィルタにて発展し，粒子フィルタにより非線形性や非ガウス性を積極的に扱えるようになった．
実シーンでは一般に複数の追跡対象が存在し，またレーダ等の原信号における対象の検出処理に伴う誤検出や未検出があり，これらを的確に表す複数ターゲッ
ト追跡の定式化が，複数のアプローチで行われている．
そこでの状態推定法として，粒子フィルタをきっかけとする実用的な状態推定アルゴリズムの研究は，今なお発展している分野である 26) ．
5.2 時系列解析
経済学や社会科学，医学や生物学などの幅広い分野
において，時系列データや時空間データに対する粒子フィルタの応用も盛んである．例えば，時系列データの，トレンドや周期パターン，不規則な周期変動などの成分への分解は，線形ガウス状態空間モデルとカル
マンフィルタでも行われていたが，粒子フィルタの登場により，非線形非ガウスの状態空間モデルが積極的に使えるようになった．
より具体的な例として，計量経済学 25) における確率的ボラティリティモデルでは，観測ノイズの分散が状態変数に依存する定式化となる．これは，例えば株価などが，その額面の大きさに応じて大きな変動をする現象を的確に捉える為であるが，強い非線形性を持つことになる．このように非線形性の強いモデリングは，カルマンフィルタ等では適切な状態推定が困難であったが，粒子フィルタはそれを可能にした．
他にも，生体の非線形現象や，量的尺度と質的尺度の混在した医学データなどを，適切に扱えるようになった．
5.3 動画像における追跡
動画像中の対象物を認識してその動きを追跡する課
題において，CONDENSATION 8) が，粒子フィルタ

として最初に顕著な成果をあげた．これは，画像における対象の検出処理を伴わない「検出前追跡」の一つの実現形態である．
近年，この分野では，深層学習により高性能な識別が可能となり，検出精度が格段に向上した為，誤検出や未検出とその対処があまり意識されなくなってきている．各時刻ごとの検出結果を，そのまま時間順に並べて利用するという，素朴な措置が，散見される．
しかし実際には，誤検出や未検出が相応の確率で生じ，検出される値（位置や大きさ等）の精度も完璧ではない為，状態推定によるフィルタリング処理が必要である．粒子フィルタは柔軟である為，深層学習ベースの高性能識別器を，新たな高性能センサとみなして，アルゴリズム中に自然に取り込むことが可能である．
画像フレームごと個別の認識結果に対して，時間的な対応付けを行う必要があるが，これは粒子フィルタをはじめとする状態推定が得意とするところである．これを，リカレント型のニューラルネットワークで行うアプローチも採られているが，学習データの質に依存したブラックボックス的な方法がボトルネックとなる．
状態推定での明示的な数式モデリングと，ニューラル
ネットワーク等でのデータからの学習モデルとを融合したアプローチが，今後更に発展すると期待したい．

ついては粒子フィルタで推定を進め，地図情報である
ランドマーク（目印）の位置座標を拡張カルマンフィ
ルタで推定するアプローチで，効率的な形式的解を導出し，それをアルゴリズムとして実装している．

5.6

センサ融合

自律移動ロボットや自動運転，その他の工学応用やビッグデータ分析などにおいて，複数の異なる種類のセンサ信号やデータを扱う必要が生じる．すなわち，複数異種のセンサ情報の融合が必要となる．これを粒子フィルタでは，尤度関数値の計算で簡潔に実現できる．
特に，センサのノイズ特性の確率分布が，センサ間で独立と想定できる場合には，センサごとに算出した尤
度関数値を掛け合わせるだけでよい．更に複雑な問題
設定についても，包括的に研究が進められている 27) ．

6

主たる貢献者

世界初の粒子フィルタ「モンテカルロフィルタ」の提案者『北川源四郎氏』（統計数理研究所 元所長，情報・システム研究機構 元機構長，東京大学 特任教授）
が，主たる貢献者の筆頭である．その研究活動を支えた『統計数理研究所』は，有意な貢献をした研究組織である．学術組織としては，『日本統計学会』が，数理的な方法論を公表し議論する場として大きく貢献した．
5.4 データ同化
更に，粒子フィルタに特化した学術組織『パーティ
自然科学の分野において，計算機シミュレーション
クルフィルタ研究会』があり，約１５年間にわたり地は有用であるが，そこに更に実験データを取り込んで， 道な活動が行われた．学際的な位置づけである「粒子シミュレーションのリアリティを向上させるアプロー
フィルタ」を主力とする学術組織は『パーティクルフィ
チが「データ同化」18) である．
ルタ研究会』以外にはなく，その活動は，学術界に留
自然現象のダイナミクス等には非線形な特性を持つ
まらず，産業界などにおいても，多くの人々に粒子フィものも多いが，カルマンフィルタではそれらを直接にルタを啓蒙し，普及させる事に貢献した．
は扱えないのに対し，粒子フィルタならばそれらの非国 際 的 に は ，情 報 融 合 に 関 す る 国 際 学 会 の線形式をそのまま使うことができる．
ISIF(International Society of Information Fuただし，実際の大規模なモデルでの計算においては， sion) と，そこが毎年主催している国際会議 Fusion
計算量が膨大となることから，カルマンフィルタを複（International Conference on Information Fusion）が数用いる等が行われている．粒子フィルタを活用した， 挙げられる．ISIF は欧米のメンバをコアとする組織でより実用的なアプローチの，今後の発展が期待されるあり，国際会議 Fusion は日本からの参加があまり多分野であると言えよう．
くないが，そこでは粒子フィルタに関する最先端の研5.5 ロボティクス
究が多数発表され，活発な議論がなされている．

自律移動ロボットの分野でも，粒子フィルタの貢献は大きい 17) ．環境中でロボットが自分の位置姿勢を把握する基本タスクにおいて，ロボットに搭載されたセンサ信号を各時刻ごとに処理するだけでは，解に多義性が生じる．すなわち，環境中の複数の箇所が解候補として挙がる為，解の確率分布は多峰性を持つ非ガウス分布となり得る．これに対し，ロボットのダイナミ
クスと，センサが持つノイズの非ガウス特性を考慮して，粒子フィルタにより推定を行う MCF(Monte Carlo
Localization：モンテカルロ自己位置推定) を用いることで，高い性能を得ることができる．ただし，ここでは環境，すなわち地図情報を既知としている．
SLAM 問題 (Simultaneous Localization And Mapping) として知られる，自己位置推定と地図学習の同時実施が，より実際的なシーンでは必要となるが，これを効果的に行うアルゴリズムが複数提案されている．
そのうち，「FastSLAM」と呼ばれる自己位置推定と地図学習の同時実施のアルゴリズムでは，ラオ・ブラック
ウェル化に基づき，ロボットの自己位置および姿勢に

7

まとめ

コトつくり至宝の評価基準である「先導力」，「規範力」，「意味力」，「解決力」の４つの観点について，これまで概説してきた粒子フィルタが，どのようにそれらの力を持っていると言えるのか，以下に列挙する．
(1) 先導力
初期の「モンテカルロフィルタ」で提案された方法論が，その後，4 節で見たように，様々な形へと進化し発展してきたことから，高い先導力があるといえる．
粒子フィルタは，理論的あるいは数学的な基盤が堅固であるため，このような発展性を持つことができたが，その特性が今後更に活用され，グラフィカルモデ
ル，強化学習，深層学習など，様々な理論・方法論との融合のプラットフォームになる事が期待される．
(2) 規範力
非線形非ガウス状態空間モデルにおける状態推定の標準的な手法として定着しており，高い汎用性を持つこともあって，5 節で見た通り，自然科学，社会科学，工学など多くの分野で活用されている．このことから，粒子フィルタは，高い規範力を持つといえる．
(3) 意味力
粒子フィルタは数理的な方法論であり，その妥当性として最適性や高い推定性能が特徴的である．これらは，必ずしも人々の共感を集めたり，社会文化活動に直結するものではないが，以下に述べる観点で，一定の意
味力を持っていると考えることがでる．
まず，最初の粒子フィルタの提案をきっかけにして，その後，多くの研究者を，世界規模で巻き込みながら，その考え方が発展してゆき，実応用にも供されてきた事が挙げられる．
次に，アルゴリズムが比較的シンプルであり，多くの人々の理解が得られやすい事が挙げられる．例えば，
カルマンフィルタでは逆行列の計算が必要になるが，粒子フィルタでは，モデル式を通常の向き（順向き）で計算すればよく，素直な考え方で計算の手続きを行う
為，アルゴリズムの理解も容易である．
更には，その名称が，「粒子フィルタ」，「モンテカルロフィルタ」，「パーティクルフィルタ」など，好感の持てる印象を伴い，その名称から興味を持つ場合もあり得る点が挙げられる．そのひとつの実現形態として，「パーティクルフィルタ」を名前に冠した研究会が組織され，約１５年間にわたる活動をつづけた事は，意味力の一つの顕れと言えよう．
(4) 解決力
粒子フィルタ以前の方法（カルマンフィルタ等）では，
線形性およびガウス性の枠内に留まる必要があったが，粒子フィルタはこの枠を取り払い，非線形非ガウスの新しい課題を解決する手段を我々に提供してくれた．このように，多くの分野において，課題解決できる範囲が拡張された事は，課題を解決に導く大きな影響力を持っていると言える．

参考文献
1) R.E.Kalman, ”A new approach to linear filtering andprediction problems”, Journal of Basic Engineering,
Vol.82(1), pp.35-45, 1960.
2) R.E.Kalman, R.S.Bucy, ”New Results in Linear Filtering and Prediction Theory”, Journal of Basic Engineering, Vol.83(1): pp.95-108, 1961.
3) D.Alspach, H.Sorenson ”Nonlinear Bayesian estimation using Gaussian sum approximations”, IEEE
Trans. on Automatic Control Vol.17, No.4, pp.439448, 1972.
4) G.Kitagawa, ”Non-Gaussian State-Space Modeling of
Nonstationary Time Series”, Journal of the American
Statistical Association, Vol.82, No.400, pp.1032-1041,
1987.
5) G.Kitagawa, ”A Monte Carlo filtering and smoothing method for non-Gaussian nonlinear state spacemodels”, Technical Report, The Institute of Statistical Mathematics, Research Memorandum, No.462,
1993.
6) N.J.Gordon, D.J.Salmond, A.F.M.Smith, ”Novel approach to nonlinear/non-Gaussian Bayesian state estimation”, IEE Proceedings F - Radar and Signal Processing, Vol.140, No.2, pp.107-113, 1993.
7) G.Kitagawa ”Monte Carlo Filter and Smoother for
Non-Gaussian Nonlinear State Space Models”, Journal of Computational and Graphical Statistics, Vol.5,
No.1, pp.1-25, 1996.

8) M.Isard, A.Blake, CONDENSATION ― Conditional
Density Propagation for Visual Tracking International Journal of Computer Vision, Vol.29, pp.5-28,
1998.
9) G.Kitagawa ”Self-Organizing State Space Model”,
Journal of the American Statistical Association,
Vol.93, pp.1203-1215, 1998.
10) M.K.Pitt, N.Shephard, ”Filtering via Simulation:
Auxiliary Particle Filters”, Journal of the American
Statistical Association, Vol.94, No.446, pp.590-599,
1999.
11) A.Doucet, S.J.Godsill, C.Andrieu, ”On sequential
Monte Carlo sampling methods for Bayesian filtering”, Statistics and Computing, Vol.10, pp.197-208,
2000.
12) J.S.Liu, ”Monte Carlo Strategies in Scientific Computing”, Springer, 2001.
13) A.Doucet, N.de Freitas, N.Gordon (eds), ”Sequential
Monte Carlo Methods in Practice” Springer, 2001.
14) R.Mahler, ”Multitarget Bayes filtering via first-ordermultitarget moments”, IEEE Trans. on Aerospaceand Electronic Systems, Vol.39, No.4, pp.1152-1178,
2003.
15) B.N.Vo, S.Singh, A.Doucet, ”Sequential Monte Carlo
Implementation of the PHD Filter for Multi-target
Tracking”, Proc. 6th Annual Conf. on Information
Fusion (Fusion 2003), pp.792-799, 2003.
16) B.N.Vo, S.Singh, A.Doucet, ”Sequential Monte Carlo
methods for Bayesian Multi-target filtering with Random Finite Sets”, IEEE Trans. Aerospace and Electronic Systems, Vol.41, No.4, pp.1224-1245, 2005.
17) S.Thrun,
W.Burgard,
D.Fox,”Probabilistic
Robotics”, The MIT Press, 2005.
18) 中村，上野，樋口，”データ同化:その概念と計算アルゴリズム”，統計数理，第 53 巻，第 2 号，pp.211-229, 2005.
19) B.N.Vo, W.K.Ma, ”The Gaussian mixture Probability Hypothesis Density filter”, IEEE Trans. Signal
Processing, Vol. 54, No. 11, pp.4091-4104, 2006.
20) R.Mahler, ”PHD filters of higher order in target number”, IEEE Trans. on Aerospace and Electronic Systems, Vol.43, No.4, pp.1523-1543, 2007.
21) B.T.Vo, B.N.Vo, A.Cantoni, ”Analytic implementations of the Cardinalized Probability Hypothesis Density Filter”, IEEE Trans. Signal Processing, Vol.55,
No.7, Part 2, pp.3553-3567, 2007.
22) C.Andrieu, A.Doucet, J.R.Holenstein, ”Particle
Markov Chain Monte Carlo Methods”, Journal ofRoyal Statistical Society, B, 72, Part3, pp.269-342,
2010.
23) G.Poyiadjis, A.Doucet, S.S.Singh, ”Particle approximations of the score and observed information matrix in state space models with application to parameter estimation”, Biometrika, Vol.98, No.1, pp.65-80,
2011.
24) T.Reineking, ”Particle filtering in the DempsterShafer theory”, International Journal of Approximate
Reasoning, Vol.52, No.8, pp.1124-1135, 2011.
25) 矢野 浩一，”粒子フィルタの基礎と応用：フィルタ・平
滑化・パラメータ推定”，日本統計学会誌，第 44 巻，第1 号 pp.189-216, 2014.
26) L.D.Stone,
R.L.Streit,
T.L.Corwin,
K.L.Bell,”Bayesian Multiple Target Tracking”, Second
Edition, Artech House, 2014.
27) R.P.S.Mahler, ”Advances in Statistical MultisourceMultitarget Information Fusion”, Artech House, 2014.
28) S.Reuter, B.T.Vo, B.N.Vo, K.Dietmayer, ”The Labeled Multi-Bernoulli Filter”, IEEE Trans. on Signal
Processing, Vol.62, No.12, pp.3246-3260, 2014.
29) J. Heng, A.Bishop, G.Deligiannidis, A.Doucet, ”Controlled Sequential Monte Carlo”, Annals of Statistics,
to appear, 2020.