統計数理（2008）
第 56 巻 第 2 号 225–234
c 2008 統計数理研究所


特集「データマイニングと統計数理」［研究ノート］

Merging Particle Filter とその特性
中野 慎也

1,2

・上野 玄太1,2

・中村 和幸

1,2

・樋口 知之1,2（受付 2008 年 1 月 4 日；改訂 2008 年 3 月 31 日）

要旨

逐次データ同化への応用が検討されつつある手法の一つに particle filter があるが，この手法はアンサンブルの縮退という問題が起こるために高次元の問題に対してはあまり有効に機能しない．そこでこの問題を回避するために merging particle filter（MPF）という手法が提案された．MPF において，フィルタ分布を表現するアンサンブルの構成粒子は，予測分布を表現するアンサンブルから抽出した複数のサンプルの重みつき和によって生成する．このとき，フィルタ分布の平均と共分散が保存されるよう適切に重みを与える必要があるが，その与え方には任意性がある．本研究では，重みの与え方によって MPF の性質がどう変化するかを調べるた
めに，2 種類の重みの与え方のもとでデータ同化実験を行った．その結果，低次元のモデルに対しては，重みのうちの一つを 1 に近い値に，その他を 0 に近い値にした方が基本的に精確な推定ができるが，比較的次元の高いモデルに対しては，そのような重みの設定では，アンサン
ブルを構成する粒子の数を多くしないとよい推定が得られないことがわかった．
キーワード： データ同化，particle filter，merging particle filter．
1.

はじめに

（e.g., Kalnay, 2003; 中村 他, 2005; Evensen, 2006）とは，地球
データ同化（data assimilation）
物理学的現象のモデリングの際に用いられる考え方で，物理法則を記述する数値シミュレー
ションモデルに実際の観測データの情報を組み込むことで，実際の現象をよく再現する信頼性の高いモデルを構築しようとすることを言う．データ同化は，すでに気象予報の精度向上などの目的で応用されているほか，さらに様々な目的での応用が検討されており，地球物理学における重要な研究トピックとなっている．一口にデータ同化と言っても様々なアプローチがあるが，中でもシステムの時間発展をシミュレーションモデルで計算しながら順次観測の情報を取り込んでいくアプローチは逐次データ同化と呼ばれ，地球物理学的現象のダイナミクスをモデリングするための強力な方法として盛んに研究がなされている．
逐次データ同化を行うためのアルゴリズムには様々なものがあり，状況に応じて使い分けられる．まず，システムのダイナミクスが線型で記述でき，かつシステムの状態の確率分布がGauss 分布で記述できる場合には，Kalman filter を用いて逐次データ同化を行うことができる．
しかし，線型の方程式で記述できるような地球物理現象はさほど多くはなく，Kalman filter を適用できる状況は限られてくる．システムのダイナミクスが非線型であっても，システム方程
式を線型近似して状態の共分散行列を計算する extended Kalman filter が使えるが，不安定性
12

統計数理研究所：〒106–8569 東京都港区南麻布 4–6–7
科学技術振興機構，戦略的創造研究推進事業

226

統計数理

第 56 巻

第 2 号 2008の強い非線型システムを扱う際には有効に機能しない（e.g., Evensen, 1992）
うえ，行列の計算を多く含むために高次元のモデルでは著しく計算時間が掛かるという問題もある．計算時間の問
題については， singular evolutive extended Kalman（SEEK）filter（Pham et al., 1998b）によって解決できるが，これも線型近似を用いており，不安定性の強いシステムには適用できない．
不安定性の強い非線型システムのモデリングには，確率分布を多数の実現値で構成されるアン
サンブルで近似するのが有効である．Ensemble Kalman filter（EnKF）（Evensen, 1994; Burgers
et al., 1998）は，確率分布をアンサンブルで近似して一期先予測に利用するアルゴリズムの一つであり，非線型モデルに対するデータ同化によく用いられている．ただ，EnKF もシステムの状態と観測との間に線型の関係が成り立つことを仮定しており，その関係が成り立たない場
合には妥当な結果が得られない場合もある．また，Kalman filter や extended Kalman filter と同様，行列の計算を多く含むために高次元のモデルでは著しく計算時間が掛かるという問題もある．アンサンブル近似に基づく別の手法として，SEEK filter の変形である singular evolutiveという方法も提案されている．これはinterpolated Kalman（SEIK）filter（Pham et al., 1998a）
EnKF よりも計算効率がよいと言われている（Nerger et al., 2005）が，システムの状態と観測との間に線型の関係が成り立たない場合は，EnKF 同様，あまり有効ではない．
そこで，システムの状態と観測との間の関係が線型であることを仮定しない方法として挙げられるのが，particle filter（PF）（Gordon et al., 1993; Kitagawa, 1993, 1996; Kitagawa andGersch, 1996; Higuchi and Kitagawa, 2000; van Leeuwen, 2003）である．PF も確率分布のアン
サンブル近似に基づくアルゴリズムの一つであるが，この方法は，システム，観測の線型性の仮定や Gauss 分布の仮定を必要としないため，広く様々なデータ同化問題に適用可能である．
PF では，アンサンブルを構成するサンプル（粒子と呼ばれる）がどれだけデータに適合しているかを尤度によって評価し，尤度の比に応じて粒子の複製をつくることで，フィルタ分布を近似するアンサンブルを得る．このとき，尤度がきわめて小さい粒子は破棄されることになる．
しかし，尤度が小さい粒子が破棄されて他の粒子の複製に置き換えられることで，フィルタリ
ングを何度も繰り返すうちにアンサンブルを構成する粒子の多くが同一もしくは互いにきわめて近い値を取るようになり，状態の確率分布がうまく表現できなくなってしまう場合があるという問題がある．これは EnKF などのアルゴリズムにはみられない問題で，アンサンブルの縮退と呼ばれる．アンサンブルの縮退は，原理的にはアンサンブルに含まれる粒子の数を十分多くすることで回避可能なのだが，一期先予測に多大な計算コストがかかるような問題の場合には粒子数を増やすこと自体が困難であり，実際に PF の応用できる範囲が極めて限定されてしまう大きな原因となっている．観測の線型性を仮定する必要がなく，かつアンサンブルの縮退が起こらないアルゴリズムとして，事後分布を一旦 Gauss 分布で近似し，その Gauss 分布から粒子をサンプリングし直す方法（Kotecha and Djurić, 2003），各粒子の示す実現値を中心とする
Gauss 分布の重ね合わせで事後分布を表現し，サンプリングし直す kernel filter などと呼ばれる
方法（Hürzeler and Künsch, 1998; Anderson and Anderson, 1999; Musso et al., 2001）など， PFから派生した手法も幾つか提案されているが，高次元のモデルでは適切な共分散を持つ Gauss
分布からのサンプリング自体に多大な計算コストを要する場合があり，これらの方法もあまり有効とは言えない．
上記のような背景から，PF の利点を活かしつつも，アンサンブルの縮退が起こりにくく，またモデルが高次元になっても計算コストがあまり増大しない方法として提案されたのが merging）Nakano et al., 2007）である． MPF では，フィルタ分布を表現するアンサ
particle filter（MPF（
ンブルの各構成粒子を，予測分布を表現するアンサンブルから複数のサンプルを抽出した上で，その重みつき和を取ることによって生成する．この重みつき和を取る際に，重みを適切に調整することで，フィルタ分布の平均と共分散の情報がアンサンブルで保持されるようにする．以Merging Particle Filter とその特性

227

下では，まず第 2 節において MPF のアルゴリズムについて説明した上で，第 3 節で幾つかのデータ同化実験を行ってその結果を示し，第 4 節でその結果を踏まえてこのアルゴリズムの特性についての考察を加える．
2.

Merging particle filter

まず，(2.1a)

xk = Fk (xk−1 , v k )(2.1b)

yk = Hk (xk , wk )のような状態空間モデルを考える．ここで，xk ，y k はそれぞれ時刻 T = tk (k = 1, . . . ) における状態ベクトル，観測ベクトルであり，v k ，w k はそれぞれシステムノイズ，観測ノイズを示
す．演算子 Fk はシステムの時間発展を記述し，Hk はある状態ベクトル，観測ノイズに対して得られるべき観測を与える．
PF を適用する場合と同様，時刻 T = tk−1 でのフィルタ分布 p(xk−1 |y 1:k−1 ) を，アンサンブ
ル {x(i)}N を用いてk−1|k−1 i=1

1  
(i)
δ xk−1 − xk−1|k−1
N i=1
N

(2.2)

p(xk−1 |y 1:k−1 ) ≈のように近似する．ここで，δ は Dirac のデルタ関数であり，N はアンサンブルを構成するサンプルの数である．ここから，次のタイムステップ T = tk での予測分布のアンサンブル近似が

1  
(i)
δ xk − xk|k−1
N i=1
N

(2.3)

p(xk |y 1:k−1 ) ≈というように得られる．ただし，予測分布のアンサンブルを構成する各粒子 x(i)k|k−1 は，アンサ(i)(i)}N を構成する各粒子にシステム方程式を Fk (xk−1|k−1 , v k ) のように適用す
ンブル {x(i)k−1|k−1 i=1ることで得られる．
(i)(i)
MPF では，フィルタ分布のアンサンブル近似 {xk|k }N
i=1 を構成する各粒子 xk|k を生成するたNめに，アンサンブル {x(i)k|k−1 }i=1 から復元抽出によって得た複数個の粒子を用いる．したがっを生成するのに n 個の粒子を用いるとすると，N 個の粒子のアンサンブルて，もし粒子 x(i)
k|k
(i)

(1,1)

(n,1)

(1,N)

(n,N)

{xk|k }Ni=1 を生成するために，n × N 個のサンプル {x̂k|k , . . . , x̂k|k , . . . , x̂k|k , . . . , x̂k|k

} を予

(i)
}N から復元抽出することになる．各粒子 xk|k は，この n × N 個
測分布アンサンブル {x(i)k|k−1 i=1
(i,j)
(1,i)
(n,i)
のサンプル {x̂k|k }i,j のうちの n 個の粒子からなるサブセット {x̂k|k , . . . , x̂k|k} の重みつき和

(2.4)

(i)

xk|k =n
j=1

によって生成する．このとき，重み αj を
(2.5a)n

αj = 1j=1(2.5b)

n

j=1

(j,i)

αj x̂k|kα2j = 1

統計数理

228

第 56 巻

第 2 号 2008

Nを満たすように与えれば，新たに得られるアンサンブル {x(i)k|k }i=1 がフィルタ分布 p(xk |y 1:k )と近似的に等しい平均・共分散を持つことになる（Nakano et al., 2007 参照）．そこで MPF では，このアンサンブルを 2 次のモーメントまでを保持したフィルタ分布の近似と見なし，次の一期先予測のステップに進む．
N
なお，MPF によって得られるアンサンブル {x(i)k|k }i=1 には，フィルタ分布の 2 次のモーメントまでの情報が近似的に保存されるが，3 次以上のモーメントの情報は保存されない．したがって，MPF ではフィルタ分布の形状についての情報が一般には失われてしまうことになる．
2 次までのモーメントを保存し，かつ 3 次以上のモーメントも保存できるような実数の重みの組み合わせ {αj }nj=1 が存在しないことは，以下のようにして確認できる．簡単のため確率変数
xk を 1 次元として，m 次のモーメントを考えると，
(xk − µk|k )m p(xk |y1:k ) dxk
(2.6)




1  
(i)
δ xk − xk|k dxk
N i=1
⎛
⎞
m
N
N 
n


1  (i)
1(j ,i)⎝=(x − µk|k )m =αj x̂ l − µk|k ⎠N i=1 k|k
N i=1 l=1 j =1 l k|k≈N(xk − µk|k )m

l


n
N

m 

1  
(j,i)
(j,i)
m
x̂
αm
−
µ
=
α
δ xk − x̂k|k dxk(xk − µk|k )m
≈
k|k
j
j
k|k
N i=1 j=1
N i=1
j=1

n
 m
≈
αj (xk − µk|k )m p(xk |y1:k ) dxk
n
N
1 

j=1

となるから，m 次のモーメントが保存されるためには，
n


(2.7)

αm
j =1

j=1

が成り立つ必要がある．ただし，式（2.6）の式変形の過程で
1    (jl ,i)x̂k|k − µk|k ≈ 0
N i=1
N

m

(unless j1 = j2 = · · · = jm )l=1と近似した．一方，0 でない実数の組 {αj } が式
（2.5b）
を満足するためには，0 < |αj | < 1 がすべての j に対して成り立たなくてはならないから，2 より大きい整数 m に対して，|αj |2 > |αj |m(2.8)

が言える．αj が実数なら α2j > 0 だから，αj > 0，αj < 0 に関わらず
α2j > αm
j .

(2.9)

これがすべての j について言えるので，(2.10)

n


α2j >

j=1

n


αm
j .

j=1

したがって，0 でない実数の組 {αj } が式
（2.5b）を満足するとき，m > 2 であれば式（2.10）から(2.11)

n

j=1

αm
j <1

Merging Particle Filter とその特性

229

となり，式
（2.7）の m 次のモーメントが保存されるための条件が成り立たなくなることがわかる．つまり，2 次のモーメントを保存するような重みで 3 次以上のモーメントが保存されることはない．ここでは確率変数 xk を 1 次元としたが， xk が多次元の場合でも m 次のモーメン
トが保存されないことが同様にして確認できる．
重みつき和を取る粒子の数 n は 3 以上の数にする必要がある．もし n = 1 ならば通常の PFそのものとなり，また n = 2 の場合も，式（2.5a）（，2.5b）
を満足するには片方の重みを 1，もう片方を 0 に取るしかなく，やはり通常の PF と同じものになってしまうからである．一方，n ≥ 3であれば，式（2.5a），（2.5b）を満たすような重み {α1 , . . . , αn } の与え方は無数に存在する．基本
的には，アンサンブルの多様性が維持されるように，すべての重みを互いに異なる値に設定するのが望ましいと考えられるが，重みをどのように与えるべきかについての具体的な指針は今のところ特に提案されていない．ただし，重みの与え方によって MPF の性質は少し変化する．
以下ではそのことを示すために 2 種類の重みの与え方のもとでデータ同化実験を行い，その結果を比較していく．
3.

データ同化実験

3.1 Lorenz 63 モデルによる実験
まず，Lorenz（1963）
によって発表された 3 変数の簡単なカオスのモデル（以下 Lorenz 63 モデ
ルと呼ぶ）
を使ってデータ同化の実験を行う． Lorenz 63 モデルは，以下のような方程式によって記述される．
(3.1a)(3.1b)(3.1c)dx= −s(x − y)dtdy
= rx − y − xzdtdz
= xy − bz.
dt

ここでは，パラメータを元の論文にしたがって s = 10，r = 28，b = 8/3 と与える．時間の 1 ス
テップは 0.01 とする．同化するデータとしては，Lorenz 63 モデルにある初期値を与えて得られた結果にノイズを加えることで人工的に生成したデータを用いた．データは 20 ステップごとに 1 回，x，y ，z の 3 成分すべてについて得られるものとし，各成分に加わるノイズは標準偏差 2 の Gaussian ノイズとした．また，MPF でも通常の PF でも，フィルタリングの際，観測 y k が得られたときの予測ベクトル xk|k−1 = (xk|k−1 , yk|k−1 , zk|k−1 ) の尤度 p(yk |xk|k−1 ) を計算する必要があるが，ここでは，観測ノイズ wk が平均 0，共分散行列が diag(σ 2 , σ 2 , σ 2 ) とな
るような Gauss 分布に従うものとして，(3.2)||y k − xk|k−1 ||21
p(yk |xk|k−1 ) = √
exp −
3
2σ 2
( 2πσ)のように尤度を計算する．ただし，σ = 3 とする．すなわち，尤度の計算の際の観測ノイズの標準偏差は，人工データに加えられた観測ノイズの標準偏差とは異なる値としている．シス
テムノイズに関しては，平均 0，共分散行列が diag(0.01, 0.01, 0.01) の Gaussian ノイズを仮定する．時間 T = t1 における x1|1 の分布は，そのときの観測データを平均とする共分散行列がdiag(16, 16, 16) の Gauss 分布になるものとした．
MPF を適用する際の重みについては，和を取る粒子の数 n を 3 としたうえで重みをそれぞれ√√
3
13 + 113 − 1(3.3)
α2 =,α3 = −
α1 = ,
4
8
8

統計数理

230
表 1.

第 2 号 2008

Root-mean-square deviations from the true state over 50,000 time steps for an
experiment using the Lorenz 63 model.

と与えた場合，及び，それぞれ(3.4)

第 56 巻

α1 =

19
,
20

√
α2 =

77 + 1,40√
α3 = −

77 − 140と与えた場合の 2 通りを考える．式（3.4）の重みは，式（3.3）の重みに比べて α1 が 1 に近くなっており，通常の PF に近い重みづけになっている．これら 2 種類の MPF と通常の PF を加えた計 3 つの方法について，粒子数 N を様々な値に変化させながらデータ同化実験を行い，それぞれどのくらいの粒子数で縮退の影響を受けるのかを調べた．
表 1 は，データ同化の結果得られた状態の推定値を人工データを生成したシミュレーションの上での実際の値と比較したものであり，推定値と実際の状態の値との差を，50,000 ステップ分について root-mean-square をとったもので示している．MPF 1 は式
（3.3）で重みを定めたもの，MPF 2 は式
（3.4）で重みを定めたものである． N が大きい場合には PF の方が精度がよい
が，N が小さい場合には概して PF よりも MPF の方が精度のよい推定ができるということが分かる．また，MPF のうちでも，より PF に近い式
（3.4）の重みづけの方が，式（3.3）の重みづ
けよりもよい精度で推定が得られることが分かる．
3.2 Lorenz 96 モデルによる数値実験
次に，Lorenz and Emanuel（1998）のより次元の大きいモデル（以下 Lorenz 96 モデルと呼ぶ）
を使った実験を行う． Lorenz 96 モデルの各成分 xj (j = 1, . . . , J) の時間発展は，以下のような方程式によって記述される
(3.5)

dxj
= (xj+1 − xj−2 )xj−1 − xj + f.
dt

但し，x−1 = xJ −1 ，x0 = xJ ，xJ +1 = x1 とし，f = 8 とする．また，モデルの次元 J は 40，時間の 1 ステップは 0.005 とした．同化するデータは，この Lorenz 96 モデル自体を走らせて得られた結果に観測ノイズを加えることで人工的に生成した．データの取得は，まずモデルを初期値(3.6a)

xj = 8.0(for j = 20)(3.6b)

xj = 8.008(for j = 20)から擾乱を十分発達させるために 2,000 ステップ走らせた上で，そこから 10 タイムステップ
ずつ行った．但し，このモデルの変数 xj のうちの j が偶数のもののみが観測されるものとし，かつ観測で得られるのは xj の絶対値のみであるものとした．データに加えた観測ノイズには，
平均 0，標準偏差 1.5 の Gauss 分布にしたがうノイズを用いた．

Merging Particle Filter とその特性
表 2.

231

Root-mean-square deviations from the true state from time step 3,000 to time
step 20,000 for an experiment using the Lorenz 96 model.

このようにして生成したデータを同化する際，システムノイズは平均 0，共分散 diag(0.25, . . . ,0.25) の Gauss 分布にしたがうものと仮定し，また尤度は2||y k − H(xk|k−1 )||1p(y k |xk|k−1 ) = √
exp −
20
2σ 2
( 2πσ)(3.7)とする．ただし，観測演算子 H は，H(xk|k−1 ) = (|x2,k|k−1 | |x4,k|k−1 | . . . |x40,k|k−1 |)T となるよ
うに定義されたものであり，また σ = 3 とする．先の Lorenz 63 モデルでの実験と同様，MPFを適用するにあたっては重みが式（3.3）で与えられる場合と式（3.4）で与えられる場合との 2 通りを考え，これに通常の PF を加えた 3 種類の方法について実験を行い，比較した．
表 2 は，データ同化の結果得られた状態の推定値を人工データを生成したシミュレーションの上での実際の値と比較したものであり，推定値と実際の状態の値との差を，最初の観測から5,000 ステップ目から 20,000 ステップ目までについて root-mean-square をとったもので示している．表 1 と同様，MPF 1 は式
（3.3）
で重みを定めたもの，MPF 2 は式
（3.4）
で重みを定めたものである．PF を用いた場合，N を 32768 まで増やしても，MPF ほど精度のよい推定はできない．一方，MPF については，2 種類を比較すると，N が大きい場合には，より PF に近い式
（3.4）の重みづけの方が，式（3.3）の重みづけよりもよい精度で推定が得られるが，N = 1024 もしくは N = 2048 では，式（3.3）
の重みを使った場合の方が式（3.4）の重みの場合よりもよい推定結果が得られている．
4.

考察

Nakano et al.（2007）
で議論されたように， MPF は PF と比べて少ない粒子数で精度よい状態の推定ができるが，前節で示されたように MPF の特性はその重みの与え方によっても変化する．MPF で重みつき和を取る際，仮に重みのうちの 1 つ α1 を 1 とし，残りの重みを 0 とす
ると通常の PF と同値になるが，そこから重みの値を少し動かして，α1 を 1 に近い値にとり，}N
残りは 0 に近い値となるようにすると，通常の PF で得られるアンサンブル {x̂(1,i)
i=1 の各粒k|k

N
子の場所にアンサンブル ( nj=2 αj x̂(j,i)
k|k )i=1 で表現される幅の小さい分布を配置した形になる．
つまりこれは kernel filter（Hürzeler and Künsch, 1998; Anderson and Anderson, 1999; Musso etal., 2001）と類似のものを MPF の文脈で実現したものと見なすことができる．MPF では，第2 節でも述べたように一般にはフィルタ分布の形状が保存されないが，このように α1 を 1 に近い値にとれば，kernel filter 同様，フィルタ分布の形状が大きく歪められることはないものと考えられる．
Lorenz 63 モデルの実験では， MPF の重みのうち， α1 を 1 に近い値にとった方がよい推定

統計数理

232

第 56 巻

第 2 号 2008

結果が得られるという傾向が現れた．これは，Lorenz 63 モデルでは，強い非線型性のために状態の確率分布が Gauss 分布とは大きく異なる形になってしまう場合が多く，フィルタ分布の形状を歪めずに扱った方が精確な状態の推定ができることを示しているものと考えられる．実際に地球物理モデルにデータ同化を適用する際にも，しばしば状態の確率分布が Gauss 分布から大きく外れてくる場合があるが，そのような場合には，こうした α1 を 1 に近い値にとるよ
うな重みづけが有効となるだろう．
一方，Lorenz 96 モデルの実験では，粒子数 N を少なくした場合，α1 を 1 に近づけない方がよい推定結果を得た．このことは，高次元のモデルにおいては，粒子数 N を少なくすると，フィルタ分布の形状を保存するよりもアンサンブルの縮退の方が問題になってくるためと考えられる．但し，粒子数 N を多く取れば，α1 を 1 に近づけた場合の方がより精確な推定値を得ることができた．これは，十分な数の粒子を使うことができるならば，フィルタ分布の形状をうまく表現できる手法を用いた方がより精確な推定ができるということを示している．勿論，α1 をあまり 1 に近づけ過ぎると，通常の PF とほとんど同じになってしまって縮退の影響が出やすくなると考えられるので，粒子数 N を多く取れる場合でも，α1 に縮退の影響が出にくい
程度の適切な値を選ぶことが大切であろう．
謝

辞

本研究は科学技術振興機構・戦略的基礎研究推進事業（JST/CREST）ならびに科学研究費補助金（19740306）の助成を受けたものである．

参 考 文 献Anderson, J. L. and Anderson, S. L.（1999）. A Monte Carlo implementation of the nonlinear filtering
problem to produce ensemble assimilations and forecasts, Monthly Weather Review, 127, 2741–
2758.
Burgers, G., van Leeuwen, P. J. and Evensen, G.（1998）. Analysis scheme in the ensemble Kalmanfilter, Monthly Weather Review, 126, 1719–1724.
. Using the extended Kalman filter with a multilayer quasi-geostrophic model,
Evensen, G.（1992）
Journal of Geophysical Research, 97（C11）, 17905–17924.
. Sequential data assimilation with a nonlinear quasi-geostrophic model using
Evensen, G.（1994）
Monte Carlo methods to forecast error statistics, Journal of Geophysical Research, 99（C5）,10143–10162.
. Data Assimilation: The Ensemble Kalman Filter, Springer-Verlag, Berlin andEvensen, G.（2006）
Heidelberg.
Gordon, N. J., Salmond, D. J. and Smith, A. F. M.（1993）. Novel approach to nonlinear/non-Gaussian
Bayesian state estimation, IEE Proceedings F, 140, 107–113.
. Knowledge discovery and self-organizing state space model,
Higuchi, T. and Kitagawa, G.（2000）
IEICE Transactions on Information and Systems, E83-D, 36–43.
. Monte Carlo approximations for general state space models,Hürzeler, M. and Künsch, H. R.（1998）
Journal of Computational Graphical Statistics, 7, 175–193.
. Atmospheric Modeling, Data Assimilation and Predictability, Cambridge UniverKalnay, E.（2003）
sity Press, Cambridge.
Kitagawa, G.（1993）. A Monte Carlo filtering and smoothing method for non-Gaussian nonlinear
state space models, Research Memorandum, No. 462, The Institute of Statistical Mathematics,

Merging Particle Filter とその特性

233

Tokyo.
. Monte Carlo filter and smoother for non-Gaussian nonlinear state space models,
Kitagawa, G.（1996）
Journal of Computational Graphical Statistics, 5, 1–25.
.
Kitagawa, G. and Gersch, W.（1996）

Smoothness Priors Analysis of Time Series, Chapter 6,

Springer-Verlag, New York.
. Gaussian particle filtering, IEEE Transactions on Signal
Kotecha, J. H. and Djurić, P. M.（2003）
Processing, 51, 2592–2601.
Lorenz, E. N.（1963）. Deterministic nonperiodic flow, Journal of the Atmospheric Sciences, 20, 130–141.
. Optimal sites for supplementary weather observations:
Lorenz, E. N. and Emanuel, K. A.（1998）
Simulations with a small model, Journal of the Atmospheric Sciences, 55, 399–414.
. Improving regularized particle filters, Sequential
Musso, C., Oudjane, N. and Le Gland, F.（2001）
Monte Carlo Methods in Practice（eds. A. Doucet, N. de Freitas and N. Gordon）, Chapter 12,
247–271, Springer-Verlag, New York.

中村和幸，上野玄太，樋口知之（2005）. データ同化：その概念とアルゴリズム，統計数理，53, 211–229.
Nakano, S., Ueno, G. and Higuchi, T.（2007）. Merging particle filter for sequential data assimilation,
Nonlinear Processes in Geophysics, 14, 395–408.
. A comparison of error subspace Kalman filters, Tellus,
Nerger, L., Hiller, W. and Schröter, J.（2005）
57A, 715–735.
Pham, D. T., Verron, J. and Gourdeau, L.（1998a）．Filtres de Kalman singuliers évolutifs pourl’assimilation de données en océanographie, Comptes rendus de l’Académie des sciences, Série
2 : Sciences de la Terre et des Planètes, 326, 255–260.
Pham, D. T., Verron, J. and Roubaud, M. C.（1998b）. A singular evolutive extended Kalman filter
for data assimilation in oceanography, Journal of Marine Systems, 16, 323–340.
. A variance-minimizing filter for large-scale applications, Monthly Weather
van Leeuwen, P. J.（2003）
Review, 131, 2071–2084.

234

Proceedings of the Institute of Statistical Mathematics Vol. 56, No. 2, 225–234(2008)

Merging Particle Filter and Its Characteristics
Shin’ya Nakano, Genta Ueno, Kazuyuki Nakamura and Tomoyuki Higuchi
The Institute of Statistical Mathematics;
Japan Science and Technology Agency

A signiﬁcant problem with the basic particle ﬁlter algorithm is degeneration. The
merging particle ﬁlter algorithm has recently been proposed to overcome this problem
at a reasonable computational cost. In an MPF, each member of a ﬁltered ensemble isgenerated from a weighted sum of multiple samples from the forecast ensemble such thatthe mean and covariance of the ﬁltered distribution are approximately preserved. In thisstudy, we performed data assimilation experiments using an MPF with two diﬀerent setsof merging weights. When one merging weight is set to near 1 and the other weights
are set small, better estimates are obtained for data assimilation into a low dimensionalmodel. For data assimilation into a relatively large dimensional model, such a weight set
requires a large ensemble size.

Key words: Data assimilation, particle filter, merging particle filter.