統計数理（2019）
第 67 巻 第 2 号 241–253©2019 統計数理研究所

創立 75 周年記念号［研究詳解］

粒子フィルタとデータ同化
上野 玄太1,2（受付 2019 年 1 月 4 日；改訂 6 月 17 日；採択 6 月 21 日）

要

旨粒子フィルタによる状態推定には，粒子と呼ばれる状態ベクトルの実現値が多数必要である．一方で，数値シミュレーションモデルを基盤とするデータ同化では，時間積分の計算コストが大きいことから，データ同化のアルゴリズムの計算コストを抑えることが要求される．多数の粒子を生成することと計算コストを抑えることは両立が難しい要求であるため，粒子フィルタによるデータ同化の方法は，現在のところ現業段階には至っておらず，研究段階にあるといってよい．粒子フィルタを適用する際の課題の一つは，各粒子に割り当てられる重みが 1 粒子に集中する，いわゆる退化の問題を限られた数の粒子でいかに克服するかである．
本稿では，重点サンプリング法により対処する方法として，陰粒子フィルタ（implicit particle
および等重粒子フィルタ（equivalent-weights particle ﬁlter）の考え方とアルゴリズムを解
ﬁlter）
説した．
キーワード：粒子フィルタ，データ同化，重点サンプリング，陰粒子フィルタ，等重粒子フィルタ．
1.

はじめに

数値シミュレーションモデルに基づくデータ同化システムの最大の特徴は，変数が多いことである．典型的には，状態ベクトル xt の次元は k = 105 ∼ 108 ，観測ベクトル y t の次元は
 = 103 ∼ 107 である．このため，計算アルゴリズムに気を配る必要があり，典型的な例を挙げ
ると，k × k 次元の密行列の数値計算を行うことはご法度とされる．これが意味することは，状態空間モデルの推定で最も有名なアルゴリズムであるカルマンフィルタ，また，その延長にある拡張カルマンフィルタは適用できないということである．そのような状況で，データ同化の方法の研究の方向性は大きく二つある．第 1 の方向性は，推定誤差を求めないことである．その代表的なアルゴリズムは 4 次元変分法（アジョイント法）である．4 次元変分法は，事後分布
最大解（MAP 解）を求めることに注力し，解の誤差評価は行わないという計算コストの選択と集中を実施することで，大規模なシミュレーションモデルに対してもデータ同化を可能にしている．第 2 の方向性は，数学的な厳密さに多少目をつぶることである．その代表はアンサンブ
ルカルマンフィルタである．アンサンブルカルマンフィルタは，推定誤差は求めることができる一方で，厳密性は欠くが変数の分布をガウス分布とみなして推定を行うアルゴリズムである．
一方で，急激な変化を伴う現象を幅をつけて予報するなど，状態変数の非ガウス分布を忠実（Kitagawa, 1996）の適用が検討されに扱いたい場合には，第 3 の方向性として，粒子フィルタ
1

統計数理研究所：〒190–8562 東京都立川市緑町 10–3

2

総合研究大学院大学 複合科学研究科統計科学専攻：〒190–8562 東京都立川市緑町 10–3

統計数理

242

第 67 巻 第 2 号 2019ている．ところが，データ同化と粒子フィルタは実際のところ相性が悪い．粒子フィルタにより適切に状態推定を行うためには，多くの粒子，具体的には粒子数を N としたときに，N  k
を満たす粒子数が必要であることが予想される．この条件は，単に k 次元ガウス分布を推定する場合であっても，正定性を満たすサンプル分散共分散行列を得るためには k + 1 個のサンプルが必要になることを考えれば，より複雑な非ガウス分布を対象にしたい以上は妥当なものであると考えられる．しかし，各粒子を生成するためにシミュレーションモデルの時間積分の実施が必要であること，その時間積分には大きな計算コストが要求されるのが一般的であることから，N  k はもちろん，N > k を要求するアルゴリズムの実施も不可能である．そのため，
たとえば N = 2k + 1 であるアンセンテッドカルマンフィルタ（Julier and Uhlmann, 2004）も適用不可である．
すべての状態変数間の共分散の正確な推定は難しいかもしれないが，みなしガウス分布のような便宜的な操作を施さず，シミュレーションモデルによる状態変数の非線形発展をそのまま生かす粒子フィルタの特徴は魅力である．実際，アルゴリズム的には粒子数 N は小さくても粒子フィルタの実行は可能であるから，粒子のばらつきから幅を持った予測もできるはずである．ところが，粒子フィルタを実際に適用してみると，観測ベクトルの次元が高いことに起因して，リサンプリングの際の重みが 1 粒子に集中してしまう，いわゆる退化の問題が深刻であることが明らかになる．1 粒子では確率分布の推定はもちろん，幅を表現することすら不可能
である．
そのため，データ同化の分野においては，第 3 の方向性として，粒子フィルタの適用というよりも，粒子フィルタの改良の研究が進められているという段階である．改良の方向は二つに大別でき，一つは粒子フィルタによる退化を緩和させ，粒子の多様性を失わせない処理を付け加えるものである（例えば融合粒子フィルタは Nakano et al., 2007，レビューは van Leeuwen,2009）．近年は，変数間の相関距離を短く打ち切る，いわゆる局所化の処理を加えた粒子フィルタの研究が盛んである（Farchi and Bocquet, 2018; van Leeuwen et al., 2019）．融合にせよ局
所化にせよ，付け加える処理は，状態空間モデルが表現する変数の数学的な構造に，多少目をつぶった結果の操作である．
もう一つの改良の方向は，重点サンプリングの方法をデータ同化の問題向けに展開するものである．ここに分類されるアルゴリズムは，変数の構造の数学的な厳密さを失わずに構築されていることが特徴である．ただし，粒子数 N をそれほど大きくはとれない以上，数学的厳密性の有無が，推定・予測性能に直結するかどうかは自明ではないことは注意しておきたい．本稿では，後者に分類されるアルゴリズムとして，陰粒子フィルタおよび等重粒子フィルタを解説する．現時点で筆者の知る限り，大規模なデータ同化システムに適した，重点サンプリングを基礎としたアルゴリズムはこの二つのフィルタに限られる．
2.

モデルと粒子フィルタ

非線形のシステムモデル，線形の観測モデルからなる状態空間モデルを考える．システムノ
イズ，観測ノイズはガウス分布に従うとする．
(2.1)

xt = f t (xt−1 ) + v t(2.2)

v t ∼ N (0, Qt )(2.3)

y t = Ht x t + w t(2.4)

wt ∼ N (0, Rt )ここで，xt , xt−1 , v t は k 次元ベクトル，y t , wt は  次元ベクトル，f t は k 次元非線形関数，

粒子フィルタとデータ同化

243

Ht は  × k の観測行列であるとする．Qt , Rt はそれぞれ，k × k,  ×  の正定値対称行列で，既知であるとする．
状態空間モデルで表現される状態の統計的推測に関しては，カルマンフィルタや粒子フィル
タに代表される逐次的な状態推定アルゴリズムが知られている．このアルゴリズムでは，一期先予測分布 p(xt |y 1:t−1 )，およびフィルタ分布 p(xt |y 1:t ) が逐次的に得られる．すなわち，初期
状態 x0 の確率分布 p(x0 ) を与えたとき，以降の時刻での状態 xt の一期先予測分布，フィルタ分布はそれぞれ



(2.5)
(2.6)

p(xt |y 1:t−1 ) =p(xt |xt−1 )p(xt−1 |y 1:t−1 )dxt−1

p(xt |y 1:t ) = p(y t |xt )p(xt |y 1:t−1 )p(y t |xt )p(xt |y 1:t−1 )dxtと与えられる．
（N 個とする）の実現値 {x(n) }N
粒子フィルタでは，x の確率分布を多数
n=1 を用いて. 1 p(x) =
δ(x − x(n) )
N
N

n=1

.
とモンテカルロ近似により表現する．ここで，= はモンテカルロ近似等号である．個々の実現
値 x(n) を粒子という．以降，一期先予測分布 p(xt |y 1:t−1 )，フィルタ分布 p(xt |y 1:t ) を近似す
(n)る粒子をそれぞれ x(n)t|t−1 , xt|t と書き，
1 (n)δ(xt − xt|t−1 )
N
N

(2.7)

p(xt |y 1:t−1 )=
˙

n=1

˙
p(xt |y 1:t )=

(2.8)

N
1 

N

(n)

δ(xt − xt|t )n=1と近似表現する．
フィルタ分布（2.6）にアンサンブル近似（2.7）
を用いると，.  (n)(n)p(xt |y 1:t ) =βt δ(xt − xt|t−1 )
N

(2.9)

n=1
(n)

(n) defβt =(2.10)p(y t |xt|t−1 )

N

n=1

(n)

p(y t |xt|t−1 )となる．この式を（2.8）
式のように，各粒子で等しい重み (1/N ) を持つ粒子 x(n)t|t で表現するにNは，重み βt(n) (n = 1, . . . , N ) で {x(n)（復元抽出）
を行い，t|t−1 }n=1 から計 N 回のリサンプリング

N得られた N 個のサンプルを {x(n)t|t }n=1 とすることで実現できる．

3.

重点サンプリングの導入

フィルタ分布（2.6）に一期先予測分布（2.5）
を代入し，xt−1 のフィルタ分布の粒子近似（（2.8）
を行うと，次のようになる．
式で t − 1 としたもの）

N

(3.1)

(n)

n=1

p(xt |y 1:t )=˙  Nn=1

p(y t |xt )p(xt |xt−1|t−1 )(n)p(y t |xt )p(xt |xt−1|t−1 )dxt

統計数理

244

第 67 巻 第 2 号 2019（3.1）式を重点サンプリングの方法に基づき近似する．x(n)t−1|t−1 をパラメータとして持つ，xt
(n,1)
(n,2)
(n,M )のある提案分布 q(xt ; x(n), xt
, . . . , xtとす
t−1|t−1 ) からの M 個の独立なサンプルを xt
ると，（3.1）式の分子および分母に現れる項は(n)(n)
p(y t |xt )p(xt |xt−1|t−1 ) =p(y t |xt )p(xt |xt−1|t−1 )(n)q(xt ; xt−1|t−1 )(n)q(xt ; xt−1|t−1 ). p(y t |xt )p(xt |xt−1|t−1 ) 1 (n,m)δ(xt − xt)=(n)Mq(xt ; xt−1|t−1 )
m=1
M

(n)

)p(xt
|xt−1|t−1 )1  p(y t |xt(n,m))δ(xt − xt
=
(n,m)
(n)
M
q(xt
;x
)
M

(n,m)

(n,m)

(n)

t−1|t−1m=1と表される．これを（3.1）
式に代入すると，.   (n,m)(n,m)
p(xt |y 1:t ) =βtδ(xt − xt
)
N

(3.2)

M

n=1 m=1となる．ここで，(n,m)

λ(n,m) def

(3.3)

t M

= N

βt

n=1

(n,m)

p(y t |xt(n,m) def

(3.4)

=

λt

(n,m)

m=1

λt

(n,m)

)p(xt

(n,m)

q(xt

(n)

|xt−1|t−1 )(n); xt−1|t−1 )と置いた．（3.2）式は，n と m に関する二重の和で，N M 個の粒子によるモンテカルロ近似になっていることに注意されたい．
(n)粒子フィルタ（2 節）は，提案分布を q(xt ; x(n)t−1|t−1 ) = p(xt |xt−1|t−1 )，各提案分布からの独立なサンプル数を M = 1 として，合計 N 個の粒子でフィルタ分布を表現したものである．一
得ている状況において方，4 次元変分法に基づき最適解を一つだけ（x(1)
t−1|t−1 ，N = 1 に相当）
も，提案分布 q(xt ; x(1)t−1|t−1 ) から M 個の粒子を生成し，フィルタ分布の表現が可能である．具の方法による（Atkins et al., 2013）．
体的には陰粒子フィルタ
（5 節）
4.

最適提案分布

提案分布として
(4.1)

(n)

(n)

q(xt ; xt−1|t−1 ) = p(xt |y t , xt−1|t−1 )を選ぶ．ベイズの定理
(n)

(4.2)

(n)

p(xt |y t , xt−1|t−1 ) =p(y t |xt )p(xt |xt−1|t−1 )(n)p(y t |xt−1|t−1 )を用いると，（3.3），（3.4）
式より，
(n)

(4.3)

(n,m)

βt

p(y t |xt−1|t−1 )

=
M

N

n=1

(n)

p(y t |xt−1|t−1 )と な る ．こ の と き の 重 み βt(n,m) は m に 依 存 し な い ．こ れ は ，n を 固 定 し た も と で は ，

粒子フィルタとデータ同化

245

m = 1, . . . , M に関する重みの分散は最小値であるゼロとなる，と言い換えることができ
る．等しい重みであれば，リサンプリングによる退化は完全に克服できる．この分散が最小を最適提案分布という（Doucet et al., 2000, Proposition 2）．たになるという性質から，（4.1）
だし，n を n = 1, . . . , N と変化させ，M = 1 とする粒子フィルタのアルゴリズムにおいてはn = 1, . . . , N に関する重みの分散は最小というわけではない．
（2.4）–のもとでは，（4.1），（4.3）式に現れる確率分布が解析的に表現
状態空間モデル（2.1）
でき，
(n)

(n)

p(xt |y t , xt−1|t−1 ) = φ(xt ; ζ t , Ut )(4.4)

p(y t |xt−1|t−1 ) = φ(y t ; Ht f t (xt−1|t−1 ), Ht Qt Ht + Rt )
(n)

(4.5)

(n)

となる．ここで，φ(x; μ, Σ) は平均ベクトル μ，分散共分散行列 Σ の正規分布に従う確率変数
x の密度関数とし，以下のようにおいている．
(n) def

ζt

(4.6)

(n)

(n)

= f t (xt−1|t−1 ) + Kt (y t − Ht f t (xt−1|t−1 ))

def

(4.7)

Ut = (I − Kt Ht )Qt(4.8)

Kt = Qt Ht (Ht Qt Ht + Rt )−1def

導出方法は，樋口 他（2011, A.3 節）を参照のこと．提案分布（4.4）からのサンプルは，Ut のコ
レスキー分解によって得られる LL = Ut を満たす k × k 下三角行列 L を用いれば，(4.9)(n,m) def

xt

(n)

= ζt

(n,m)

+ Lξ tn = 1, . . . , N ; m = 1, . . . , M

として生成できる（北川, 2005, 16.3 節）．ただし，ξ(n,m)
は k 次元標準正規分布に従う確率変数tξ t ∼ N (0, I)(4.10)の実現値である．
5.

陰粒子フィルタ

ところが，（4.9）
式に基づいて x(n,m)
を生成するには非常に計算コストがかかる．最も計算t式の右辺で用いる L を得るための，コレスキー分解 LL = Ut でコストがかかる部分は，（4.9）
ある．k × k 次元の密行列のコレスキー分解はデータ同化分野ではご法度とされている演算である．そこで，計算コストを抑えることを意識して，（4.9）式と同値な関係式を導くことを考える．
5.1 サンプリングの方法
は，ξ(n,m)を与えたもとで，次の関係式を満たすものと考える（4.9）式を変形すると，x(n,m)
ttことができる．
1 (n,m)1 (n,m)  (n,m)
(n)
(n,m)
(n)
(5.1)
− ζ t ) Ut (xt− ζ t ) = (ξ t
) ξt
(x
2 t
2

ここで，左辺はk
1 (n,m)
k
1
1
(n)
(n,m)
(n)
− ζ t ) Ut (xt− ζ t ) − log 2π − log |Ut |log 2π + log |Ut | + (xt
2
2
2
2
2
(n,m)
(n)
(n)
(n)
; ζ t , Ut ) − [− log φ(ζ t ; ζ t , Ut )]= − log φ(xt(n,m)

= − log φ(xt

(n)

(n)

; ζ t , Ut ) − min[− log φ(xt ; ζ t , Ut )]
xt

統計数理

246

第 67 巻 第 2 号 2019と書けることに着目する．
（4.4）
式およびベイズの定理（4.2）を用いると，この式はさらに，(n,m)− log p(xt

(n)

xt

(n,m)

= F (n) (xt

(5.2)

(n)

|y t , xt−1|t−1 ) − min[− log p(xt |y t , xt−1|t−1 )](n)) − Fminとなる．ここでdef(n)F (n) (xt ) = − log p(y t |xt )p(xt |xt−1|t−1 )(5.3)(n) def

Fmin = min F (n) (xt )

(5.4)

xt

と定義した．したがって，（5.1），（5.2）
式より，k 次元標準正規分布に従う確率変数の実現値(n,m)ξt(n = 1, . . . , N ; m = 1, . . . , M ) を与えたもとで，(n,m)

F (n) (xt

(5.5)

(n)

) − Fmin =1 (n,m)  (n,m)
) ξt
(ξ
2 t

を解くことで，x(n,m)
を生成することができる．この方法を陰粒子フィルタ（implicit particletという（Chorin et al., 2010, 2013; Atkins et al., 2013）．（5.5）式を計算コストを抑えて解く
ﬁlter）
方法が，Morzfeld and Chorin（2012）にある．
（2.4）–のもとでは，（5.3）式は状態空間モデル
（2.1）
(5.6)

F (n) (xt ) =1
1
log 2π + log |Rt | + (y t − Ht xt ) Rt−1 (y t − Ht xt )
2
2
2
1
k
1
(n)
(n)
+ log 2π + log |Qt | + (xt − f t (xt−1|t−1 )) Q−1
t (xt − f t (xt−1|t−1 ))222となる．第 3 項，第 6 項以外は xt に依存しないので，変分法で通常用いる目的関数
def 1
(5.7)
J (n) (xt ) = (y t − Ht xt ) Rt−1 (y t − Ht xt )
2
1
(n)
(n)
+ (xt − f t (xt−1|t−1 )) Q−1
t (xt − f t (xt−1|t−1 ))2(5.8)(n) def

Jmin = min J (n) (xt )
xtを用いると，（5.5）
式の代わりに
(5.9)

(n,m)

J (n) (xt(n)) − Jmin =1 (n,m)  (n,m)
) ξt
(ξ
2 t

を解いて x(n,m)
を求めればよい．
t
5.2 重みの計算（3.3），（3.4）式に基づき，重みを求める．
) を考える．
（4.9）式より，xt は標準正規分布（3.4）式の分母に現れる提案分布 q(xt ; x(n)t−1|t−1（4.10）に従う確率変数 ξt から，次の関係を経て生成される確率変数であると解釈できる．
(5.10)(n)xt = ζ t+ Lξ t（5.10）式により生成される xt の確率分布，すなわちサンプリングを行う提案分布 q(xt ; x(n)t−1|t−1 )は(n)φ(ξ t ; 0, I) = q(xt (ξ t ); xt−1|t−1 ) ∂xt (ξt )  ∂ξt 粒子フィルタとデータ同化

247

(n)

= q(xt (ξ t ); xt−1|t−1 )L

(5.11)

を満たす．（5.11）
式，および（5.5）
式を用いると，（3.4），（3.3）式より，重みは次のように表さ
れる．
(n,m)

(n,m)

λt

=

exp{−[− log p(y t |xt

(n,m)

q(xt

(n,m)

=

exp[−F (n) (xt(n); xt−1|t−1 )

)]

L



1 (n,m)  (n,m)
) ξt
(ξ
2 t

 L
1 k/2
1 (n,m)  (n,m)
exp − (ξ t
) ξt
2π
2
(n)

exp −Fmin −= (2π)k/2 exp(−Fmin )L(n)(5.12)(5.13)(n)|xt−1|t−1 )]}

(n,m)
φ(ξ t
; 0, I)



=

(n,m)

)p(xt

(n,m)

βt

(n)

exp(−Fmin )

=
M

N

n=1

(n)

exp(−Fmin )状態空間モデル
（2.1）
（2.4）
–
のもとでは，重み（5.13）は
(n,m)

βt

(5.14)(n)exp(−Jmin )

=
M

N

n=1

(n)

exp(−Jmin )となる．
特に，4 次元変分法のように N = 1 の場合で，m = 1, . . . , M によって多数の粒子を生成すは m によらない定数になるため，生成した粒子 M 個をリサる場合には，重み（5.13），（5.14）
ンプリングなしでそのままフィルタ分布を表現する粒子として扱うことができる．
6.

等重粒子フィルタ

，（5.13）式で与えられるように，最適提案分布ならびに陰粒子フィルタでは，重み βt(n,m)（4.3）
は m に依存しないが，n を変えれば重みは x(n)t−1|t−1 に依存して変化する．そのため，M = 1
として N 個の粒子を扱う粒子フィルタでは，最適提案分布から受ける恩恵は限定的である．
実際，観測ベクトル y t の次元  が大きい場合には，最適提案分布を用いた重みであっても退化が起こる．等重粒子フィルタ（equivalent-weights particle ﬁlter, van Leeuwen, 2010; Ades andは，重み βt(n,m) を n にも依存しない定数になるように粒子を生成す
van Leeuwen, 2013, 2015）
る方法である．以降の節で見るように，得られる重みは厳密に等しくなるわけでなく，ほとん
ど等しい重みとなる．そのため，「ほとんど等重粒子フィルタ」と呼ぶのが適切で，提案当初
は almost が添えられていた（van Leeuwen, 2010）が，最近は単に等重粒子フィルタと呼ばれている．
等重粒子フィルタも粒子フィルタと同様に，M = 1 として N 個の粒子を扱う．そのため，本節の以下の説明では m = 1 の場合のみを考えているとご理解頂きたい．
6.1 サンプリングの考察を表現する（3.4）式は，提案分布にもとづき生成される粒子
重みの構成要素である λ(n,m)
t
(n,m)
xt
の関数の分数として表される．そこで，提案分布をうまく設計して，分子と分母それ

統計数理

248

第 67 巻 第 2 号 2019

ぞれを粒子に依存しない値にすることを考える．λ(n,m)の値が粒子に依存しなければ，重み
t
(n,m)
βt
も n, m によらない一定値となる．
以下で見るように，まず，分子が一定値となるために，各粒子が満たすべき条件を求める．
提案分布はその条件を満たす粒子を生成しやすいもの，すなわち，条件を満たす粒子の確率密度の値が大きいものを設定する．次に，分母も一定値とするために，これらの条件を満たす粒子の確率密度は共通の値をとるとする．
（4.2）
および
（4.4），（4.5）式を用いて変形すると，（3.4）式の分子を，ベイズの定理(n,m)

p(y t |xt(n,m)

= p(xt

(n,m)

)p(xt

(6.1)

=

(n)

(n)

|y t , xt−1|t−1 )p(y t |xt−1|t−1 )(n,m)= φ(xt(n)|xt−1|t−1 ); ζ t , Ut )φ(y t ; Ht f t (xt−1|t−1 ), Ht Qt Ht + Rt )

1 k/2
2π

(n)

1
|Ut |

(n)

1
2π

/2

1
|Ht Qt Ht + Rt |(n,m)

exp[−J (n) (xt)]となる．ここで，(6.2)

J (n) (xt ) =

(6.3)

ηt

1
1 (n)
(n)
(n)
(n)
(xt − ζ t ) Ut−1 (xt − ζ t ) + (η t ) (Ht Qt Ht + Rt )−1 η t
22(n) def(n)= y t − Ht f t (xt−1|t−1 )と置いた．（6.2）
式で（5.7）式と同じ J (n) (xt ) を用いているのは，これらの 2 式が恒等式である

−1 (n)
式より，J (n) (xt ) は xt = ζ (n)のときに最小値 12 (η (n)ηtことによる．（6.2）
tt ) (Ht Qt Ht + Rt )をとることが分かる．また，xt をうまく選ぶことで，J (n) (xt ) はこの最小値以上の任意の値をとることが可能であることに注意しておく．
(n,m)
式に含まれる J (n) (x(n,m)) がある定数 C となるように xtを生成するこ
ここでは，
（6.1）
t
(n,m)
(n)
の表現方法には任意性があるが，ある実数 αを導入してとを考える．xt
(6.4)

(n,m)

xt

(n)

(n)

= f t (xt−1|t−1 ) + α(n) Kt η tと表すことにする．
（6.2），（6.4）
式を用いると，J (n) (x(n,m)) = C となるために α(n) が満たすtべき条件は，
1
1 (n)
(n)
(n)
(n)
(6.5)
C = (α(n) − 1)2 (η t ) Kt Ut−1 Kt η t + (η t ) (Ht Qt Ht + Rt )−1 η t
2
2と な る ．さ ら に ，逆 行 列 補 題（ 例 え ば 樋 口 他, 2011, A.1 節 ）に よ り ，Kt = (Q−1
+
t
−1
 −1
−1
 −1
+
H
R
H
)
と表されることを用いると，
K
U
=
R
H
Ht Rt−1 Ht )−1 Ht Rt−1 , Ut = (Q−1
t
t
t t
t t
t
t
となることから，
（6.5）
式は，
1
1 (n)
(n)
(n)
(n)
(n)
(n)
(6.6)
C = (α(n) )2 (η t ) Rt−1 Ht Kt η t − α(n) (η t ) Rt−1 Ht Kt η t + (η t ) Rt−1 η t
2
2となる．ここで，右辺第 3 項の導出においては，次の変形を行った．
Rt−1 Ht Kt + (Ht Qt Ht + Rt )−1 = Rt−1 Ht QHt (Ht Qt Ht + Rt )−1 + (Ht Qt Ht + Rt )−1= (Rt−1 Ht QHt + I)(Ht Qt Ht + Rt )−1= Rt−1 (Ht QHt + Rt )(Ht Qt Ht + Rt )−1= Rt−1（6.6）式は α(n) についての 2 次方程式なので，厳密解が得られ，

粒子フィルタとデータ同化

α(n) = 1 ±

(6.7)

1−

(η t ) Rt−1 η t
(n)

249

(n)

− 2C

(n)
(n)
(η t ) Rt−1 Ht Kt η tとなる．（6.7）式には複号で示される 2 根があるが，正符号の根を採用する．α(n) が確実に正の値となることで各粒子が観測 y t に近づくからである（Ades and van Leeuwen, 2013, Section．
3.3）
式を用いて粒子（6.4）を生成すれば，（6.2）式は J (n) (x(n,m)) = C とな
このようにして，（6.7）
tり，（6.1）式は
(6.8)

1 k/2
2π

1

1
2π

|Ut |

/2

1
|Ht Qt Ht + Rt |exp(−C)となる．すなわち，（3.4）
式の分子は x(n,m)に依存しない値をとる．
t
定数 C は，N 個の粒子のうちいくつをフィルタ分布のアンサンブル近似に残したいかというユーザーの意向に応じて決めることができる．例えば，N = 100 で，80 個の粒子を残したいな(n) −1 (n) 100
1
η t }n=1 から，らば，J (n) (xt ) 最小値の集合 {minxt J (n) (xt )}100
n=1 = { 2 (η t ) (Ht Qt Ht + Rt )(n)
値が 80 番目に小さい要素を選び，その値を C80 とする．minxt J (xt ) の値が C80 以下の 80
個の n に対しては，（6.7）
式により α(n) を選び，
（6.4）式を計算することで，J (n) (x(n,m)) = C80t(n,m)(n)となる xtが生成される．minxt J (xt ) の値が C80 より大きい残り 20 個の n に対して
(n,m)
（重みは C80 から決まる共は，その 20 個の予測粒子 f t (x(n)
t−1|t−1 )，生成された 80 個の xt通の値）を合わせた計 100 個の粒子から 20 個をリサンプリングして補充する．
6.2 提案分布の設計とサンプリング
式の分母を考える．分母には提案分布が入るが，その分布の性質として期待
つづいて，
（3.4）
で与えられる x(n,m)
が生成されやすいこと，そのときの確率密度の値がされることは，（6.4）
tn によらず一定になることである．素朴には，デルタ関数を用いた
(6.9)

(n)

(n)

(n)

q(xt ; xt−1|t−1 ) = δ(xt − f t (xt−1|t−1 ) − α(n) Kt η t )がこの条件にふさわしそうである．提案分布（6.9）を用いれば，必ず xt は（6.4）で与えられる
(n,m)
xt
となり，
（6.8）
式から決まる各粒子の重みは等しくなる．
しかし，（6.9）式は，重点サンプリングの基本的な条件として，提案分布の台が分子の台を含むという条件を満たしていないという点で問題がある．そこで，重みが厳密に等しくなるという条件を緩和した次善策として，(6.10)

(n,m)

xt

(n)

(n)

= f t (xt−1|t−1 ) + α(n) Kt η t(n,m)+ Q1/2 ξ tと x(n,m)を生成することにする．ここで，ξ(n,m)
は次の密度関数を持つ確率変数 ξ t の実現値ttであるとする．
k(6.11)defp(ξ t ; , γU , γN ) = (1 − )
i=1



1
2γU


2
1l(|(ξ t )i | ≤ γU ) + φ(ξ t ; 0, γN
I)

(n,m)
ここで，(ξ t )i は ξ t の i 番目の要素を表す．
（6.10）式の右辺第 3 項として，ノイズ Q1/2
t ξtが導入されているところが，台の問題を回避するポイントである．
（6.11）式で与えられる ξ t の密度関数は，ξ t = 0 周りに限定した [−γU , γU ]k に台を持つ一様
分布と，平均が ξ t = 0 である正規分布の混合分布である．混合比を表すパラメータ 0 < < 1を小さく設定することで，確率密度の値は一様分布で与えられる n に依存しない一定の値

統計数理

250

第 67 巻 第 2 号 2019

(1 − )/(2γU )k でほとんど決まることになる．Ades and van Leeuwen（2013）では， = 10−3 /Nを採用している．γU > 0, γN > 0 も小さい数を設定することで，確率密度を ξ t = 0 の近傍に集中させることができる．具体的な値は 6.3 節で紹介する．また，懸案であった台の大きさは，正規分布の台 (−∞, ∞)k が分子の台を確実に含むことで解決している．（6.10）での Q1/2の乗t
算は，各成分が独立な ξt に相関を持たせ，状態変数間の物理的なバランスを考慮したい意図がある．
6.3 重みの計算（3.4）式に基づき，重みを求めよう．
（6.1）
式で与えられ，（6.1）式に含まれる J (n) (x(n,m))は（6.2）式で与えまず，（3.4）式の分子はtを（6.2）
式に代入し，（6.2）と（5.7）が恒等式であることを用いると，られる．粒子（6.10）(n,m)

J (n) (xt)1(n,m)  −1(n,m)
) Rt (y t − Ht xt)= (y t − Ht xt
2
1
1/2
(n)
+ [α(n) Qt Ht (Ht Qt Ht + Rt )−1 η t + ξ (n,m) ]
2
1/2
(n)
× [α(n) Qt Ht (Ht Qt Ht + Rt )−1 η t + ξ (n,m) ](6.12)となる．（6.10）
式で与えられる粒子 x(n,m)
は，（6.12）式が C に近いほぼ定数値になるように設t
計されている．その結果，（6.12）
式を用いる（6.1）式もほぼ定数，すなわち（3.4）式の分子もほぼ定数になる．
式を得たときと同様に考えると，粒子の生成（6.10）は，
次に，（3.4）式の分母を考える．（5.10）
確率変数 xt と ξt の関係
(n)

(n)

xt = f t (xt−1|t−1 ) + α(n) Kt η t(6.13)+ Q1/2 ξ tに対応し，∂xt /∂ξt = Q1/2 となる．よって，ξ t の密度関数（6.11）を用いると，（3.4）式の分母は次のようになる．
(n,m)

q(xt

(6.14)(n); xt−1|t−1 ) =1−=(2γU )k



1

(n,m)

; , γ U , γN )
p(ξ t


 ∂xt (n,m)   (ξt
)
 ∂ξt

(n,m)

|Q|

1l(ξ t∈ [−γU , γU ]k ) k/2 
+

2
π

1−

γU
γN

k



1
(n,m)  (n,m)
exp − 2 (ξt
) ξt
2γN



よって，重みは（3.3），（3.4）
式より，λ(n,m)の定数部分を除いて定義する λ̃(n,m)
を用いて，
t
t
次のように表される．
(n,m)(6.15)

βt

= N
n=1

(6.16)

(n,m)

λ̃

t M

m=1

(n,m)

λ̃t

(n,m)

exp[−J (n) (xtdefλ̃(n,m) = k/2 

1l(ξ

(n,m)

∈ [−γU , γU]k ) +

1−

2
π

γU
γN

)]

k



1
exp − 2 (ξ (n,m) ) ξ (n,m)
2γN

粒子フィルタとデータ同化

251

ただし，
（6.16）
式の分子において，J (n) (x(n,m))は（6.12）式で，x(n,m)
は（6.10）式でそれぞれ与tt
えられるものとする．
γU , γN の設定値について述べておく．
（6.16）式に見るように，γU は一様分布の台を，また，
比 γU /γN が正規分布の寄与を制御している．Ades and van Leeuwen（2013）では，γU = 10−5と設定している．正規分布の寄与は重みの値の違いをもたらすので，退化を起こさないようにしたいという観点からは，できるだけ小さいほうが望ましい．すでに は小さい値を設定しているので（ = 10−3 /N など，6.2 節），(2/π)k/2 (γU /γN )k = 1 であれば十分である．これより，
γN = (2/π)γU と設定することができる．
謝

辞査読者からの有益なコメントに感謝いたします．本研究は統計数理研究所共同研究プログラの助成を受けたものです．
ム（27 - 共研 - 2007）

参

考

文

献

Ades, M. and van Leeuwen, P. J. (2013). An exploration of the equivalent weights particle ﬁlter, Quarterly Journal of the Royal Meteorological Society, 139(672), 820–840, DOI:10.1002/qj.1995.
Ades, M. and van Leeuwen, P. J. (2015). The equivalent-weights particle ﬁlter in a highdimensional system, Quarterly Journal of the Royal Meteorological Society, 141(687), 484–503,
DOI:10.1002/qj2370.
Atkins, E., Morzfeld, M. and Chorin, A. J. (2013). Implicit particle methods and their connection withvariational data assimilation, Monthly Weather Review, 141, 1786–1803, DOI:10.1175/MWR-D12-00145.1.
Chorin, A., Morzfeld, M. and Tu, X. (2010). Implicit particle ﬁlters for data assimilation, Communications in Applied Mathematics and Computational Science, 5(2), 221–240.
Chorin, A. J., Morzfeld, M. and Tu, X. (2013). A survey of implicit particle ﬁlters for data assimilation,
State-Space Models: Applications in Economics and Finance (eds. Y. Zeng and S. Wu), Statistics
and Econometrics for Finance 1, 63–88, Springer, New York, DOI:10.1007/978-1-4614-7789-1_3.
Doucet, A., Godsill, S. and Andrieu, C. (2000). On sequential Monte Carlo sampling methods forBayesian ﬁltering, Statistics and Computing, 10(3), 197–208.
Farchi, A. and Bocquet, M. (2018). Review article: Comparison of local particle ﬁlters and new implementations, Nonlinear Processes in Geophysics, 25(4), 765–807, https://www.nonlin-processesgeophys.net/25/765/2018/, DOI:10.5194/npg-25-765-2018.
樋口知之, 上野玄太, 中野慎也, 中村和幸, 吉田 亮 (2011).『データ同化入門』, 予測と発見の科学, 6, 朝倉書店, 東京.
Julier, S. J. and Uhlmann, J. K. (2004). Unscented ﬁltering and nonlinear estimation, Proceedings ofthe IEEE, 92, 401–422, DOI:10.1109/JPROC.2003.823141.
Kitagawa, G. (1996). Monte Carlo ﬁlter and smoother for non-Gaussian nonlinear state space models,
Journal of Computational and Graphical Statistics, 5(1), 1–25.
北川源四郎 (2005).『時系列解析入門』, 岩波書店, 東京.
Morzfeld, M. and Chorin, A. J. (2012). Implicit particle ﬁltering for models with partial noise, and an
application to geomagnetic data assimilation, Nonlinear Processes in Geophysics, 19(3), 365–382,
https://www.nonlin-processes-geophys.net/19/365/2012/, DOI:10.5194/npg-19-365-2012.
Nakano, S., Ueno, G. and Higuchi, T. (2007). Merging particle ﬁlter for sequential data assimilation,
Nonlinear Processes in Geophysics, 14, 395–408.
van Leeuwen, P. J. (2009). Particle ﬁltering in geophysical systems, Monthly Weather Review, 137(12),4089–4114, DOI:10.1175/2009MWR2835.1.

252

統計数理

第 67 巻 第 2 号 2019

van Leeuwen, P. J. (2010). Nonlinear data assimilation in geosciences: An extremely eﬃcientparticle ﬁlter, Quarterly Journal of the Royal Meteorological Society, 136(653), 1991–1999,
DOI:10.1002/qj.699.
van Leeuwen, P. J., Künsch, H. R., Nerger, L., Potthast, R. and Reich, S. (2019). Particle ﬁlters forhigh-dimensional geoscience applications: A review, Quarterly Journal of the Royal Meteorological Society, 145, 2335–2365, DOI:10.1002/qj.3551.

Proceedings of the Institute of Statistical Mathematics Vol. 67, No. 2, 241–253 (2019)253

Particle Filter and Data Assimilation
Genta Ueno1,2
1 The Institute of the Statistical Mathematics
2 Department of Statistical Science, School of Multidisciplinary Sciences, The Graduate University for

Advanced Studies, SOKENDAIState estimation by the particle filter requires many realizations of the state vector,
called particles, for representing non-Gaussian distributions. On the other hand, data assimilation based on a numerical simulation model cannot adopt an assimilation algorithm
that is computationally expensive because computational resources need to be assigned
to time integration of the state vector by the simulation model. Therefore, the methodology for data assimilation by the particle filter is considered to still be in the research
stage. One important issue is how to avoid so-called the filter degeneracy, in which weights
for resampling concentrate to a single particle. This article reviews particle filter algorithms on the basis of importance sampling, specifically the implicit particle filter and theequivalent-weights particle filter.

Key words: Particle ﬁlter, data assimilation, importance sampling, implicit particle ﬁlter, equivalent-weights
particle ﬁlter.