統計数理（2018）
第 66 巻 第 2 号 339–351©2018 統計数理研究所［統計ソフトウェア］

P3：Python による並列計算機用粒子フィルタライブラリ
中野 慎也1,2 ・有吉 雄哉1,3 ・樋口 知之1,2（受付 2018 年 2 月 20 日；改訂 8 月 13 日；採択 10 月 2 日）

要

旨は，多数の粒子を用いたモンテカルロ計算に基づく状態推定手法であり，粒子フィルタ（PF）
非線型，非ガウスの問題に適用できることから広く様々な目的で用いられるようになってきている．一方で，PF には，推定に必要となる粒子の数が状態変数の自由度に対して指数関数的に増大するため，計算量も指数関数的に増大してしまうという欠点がある．並列計算機の利用
は，PF の計算量に対処する手段の一つとして有効であると考えられる．しかし，並列計算機を使うには並列プログラミングの知識が必要であり，また，PF には並列化の困難な処理が含まれているため，並列プログラミングの知識があるユーザにとっても，PF で高い並列化効率を実現するのは容易ではない．そこで，並列化効率の高い PF アルゴリズムを手軽に利用でき
3
Python Parallelized Particle Filter Library）という Python ライブラリ
るようにするために P（を開発した．本稿では，P3 で利用できる PF の並列アルゴリズムについて述べ，構成の概要や特徴を紹介する．
キーワード：粒子フィルタ，並列計算，Python．
1.

はじめに

（Gordon et al., 1993; Kitagawa, 1993, 1996; Doucet
粒子フィルタ（particle ﬁlter; 以下 PF）
は，非線型・非ガウスの状態空間モデルに適用可能な状態推定手法で，非線型時系
et al., 2001）
列解析や動画上のターゲット追跡，さらにはデータ同化など，様々な目的で用いられる．PF では，状態変数の確率分布を多数の粒子で表現し，モンテカルロ法の考え方に基づいて逐次ベイズ推定の計算を行う．しかし，PF では，状態変数の自由度（独立と見なせる状態変数の数）が大きくなると，いわゆる次元の呪いの問題が顕在化するという問題がある（Daum and Huang,
2003; Bengtsson et al., 2008; Snyder et al., 2008）．また，詳細な事後分布の情報を得たい場合には，さらに多数の粒子を用いる必要があり，12 変数の推定に 108 個の粒子を使用した事例もある（Nakamura et al., 2009）．PF の計算量は，少なくとも粒子数 N のオーダーになるため，
PF では計算時間が大きな問題となってくる．
多数の粒子を用いて高速に推定を行う手段としては，並列計算機の利用が考えられる．しかし，マルチコア，マルチノードの並列計算機が安価に購入できる状況にある一方，それを使い
1

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

2

総合研究大学院大学 複合科学研究科：〒 240–0193 神奈川県三浦郡葉山町湘南国際村

3

現 日本文理大学 工学部：〒 870–0397 大分県大分市一木 1727

統計数理

340

第 66 巻 第 2 号 2018

こなすには並列プログラミングの知識が必須であり，誰もが直ちに並列計算機上で PF を実装できるものではない．また，PF で通常用いられるフィルタリング手続きには，並列化の困難な処理が含まれているため，並列プログラミングの知識があるユーザにとっても，高い並列化
効率を実現するのは容易ではない．フィルタリング手続きを改良し，高い並列化効率を実現す
る方法はいくつか提案されている（Bolić et al., 2005; Nakano, 2010; Nakano and Higuchi, 2010）
ものの，その手続きは煩雑であり，プログラムの実装に掛かる手間の問題は依然として残る．
3
P（
Python Parallelized Particle Filter Library）
は，並列化効率の高い PF 手法を利用しやすくするために，Python ライブラリとして整備したものである．Python は，インタプリタ型のプログラミング言語ではあるが，numpy, scipy をはじめとする高速な数値計算用ライブラリや，
統計計算，機械学習のライブラリが豊富に用意されている他，mpi4py のような並列計算の手段も用意されており，科学技術計算の分野にもかなり普及した感がある．一般的な状態空間モデルに適用可能な逐次ベイズ推定手法を実装した Python ライブラリとしては，すでに FilterPy（Labbe, 2015）などがあるが，P3 では並列計算機を活用して比較的大きな規模の非線型・非ガウス状態空間モデルを扱うことを想定している．以下では，まず PF の基本的な形のアルゴリ
ズムについて述べ，P3 で実装されている並列計算用の PF アルゴリズムについて説明した後，P3 の構成や特徴を述べる．
2.

粒子フィルタの基本的なアルゴリズム

2.1 非線型状態空間モデル
時刻 tk の状態を xk ，時刻 tk に得られる観測を y k と表す．時刻 tk−1 から時刻 tk の間の状
態遷移は，関数 f k を用いて
(2.1)

xk = f k (xk−1 ) + v kという形のシステムモデルで記述されるものとする．但し，確率変数 v k は確率的な変動を表し，システムノイズと呼ばれる．一方，観測 y k は，状態 xk との間に以下の関係があるものと
する:
(2.2)

y k = hk (xk ) + wk .

関数 h は xk を観測に対応づける関数，wk は観測ノイズである．式（2.1），（2.2）をまとめて，非線型状態空間モデルと呼ぶ．
この非線型状態空間モデルは，以下に示す確率分布の形に書き直すこともできる:(2.3a)

xk ∼ p(xk |xk−1 ),(2.3b)y k ∼ p(y k |xk ).

例えば，式（2.1）
は，v k ，wk がガウス分布(2.4a)v k ∼ N (v k ; 0, Q),(2.4b)

wk ∼ N (wk ; 0, R)に従うとき，(2.5a)p(xk |xk−1 ) = N (xk ; f k (xk−1 ), Q),(2.5b)p(y k |xk ) = N (y k ; hk (xk ), R)と書くことができる．

P3 ：Python による並列計算機用粒子フィルタライブラリ341

式（2.3）の形を導入すると，時刻 tk までの観測データ y 1:k = {y 1 , . . . , y k } が与えられた時のxk は，t = 0 における状態 x0 の分布 p(x0 ) を与えておけば，次の 2 つの式を用いて推定することができる:

(2.6a)(2.6b)p(xk |y 1:k−1 ) =p(xk |y 1:k ) = p(xk |xk−1 )p(xk−1 |y 1:k−1 ) dxk−1 ,p(y k |xk )p(xk |y 1:k−1 )p(y k |xk )p(xk |y 1:k−1 ) dxk.

ここで，p(xk |y 1:k−1 ) は予測分布，p(xk |y 1:k ) はフィルタ分布と呼ばれる．式（2.6a）では予測分布を得るのにフィルタ分布が用いられ，式（2.6b）では逆にフィルタ分布を得るのに予測分布が（2.6b）を時刻 t0 から交互に適用することで，各，用いられていることに注意すると，式（2.6a）（2.6a），（2.6b）を逐次的に適用して，時刻のフィルタ分布 p(xk |y 1:k ) が得られる．このように式x1 , . . . , xk を推定する方法を逐次ベイズ推定と呼ぶ．PF は逐次ベイズ推定の枠組みに従って
xk を推定するアルゴリズムの一つである．
2.2 粒子フィルタ（PF）の概要
PF では，確率分布 p(xk |y 1:k−1 ) や p(xk |y 1:k ) を N 個の粒子で表し，式（2.6）を近似計算する．PF には様々な変形版が存在する（e.g., Doucet et al., 2001; van Leeuwen, 2009）が，最も基本的なアルゴリズムは次に述べるとおりである．
の計算をモンテカルロ法で行う．時刻 tk−1 において，p(xk−1 |y 1:k−1 ) に従う
まず，式（2.6a）(i)(i)
N 個のサンプル（粒子）{xk−1|k−1 }N
i=1 が得られていたとする．p(xk |xk−1|k−1 ) に従う乱数は，(i)p(v k ) から生成した乱数 v k を用いて
(2.7)

(i)

(i)

(i)

xk|k−1 = f k (xk−1|k−1 ) + v k

から得られる．式（2.7）にしたがい，各 i に対して x(i)k|k−1 を生成すれば，予測分布 p(xk |y 1:k−1 )Nに従う N 個の粒子 {x(i)k|k−1 }i=1 が得られる．

N
次に式（2.6b）にしたがって，フィルタ分布を求める．式（2.6b）は，{x(i)k|k−1 }i=1 の各粒子に尤

度 p(y k |x(i)k|k−1 ) で重みづけすれば，フィルタ分布 p(xk |y 1:k ) の形状を表現できることを示して(i)いる．ここで，各粒子 x(i)k|k−1 が p(y k |xk|k−1 ) に比例する確率
(i)

(2.8)

p(y k |xk|k−1 )(i)βk = N(i)p(y k |xk|k−1 )i=1N
で抽出されるように N 回復元抽出を行い，新たに N 個の粒子の集合 {x(i)k|k }i=1 を作る．そう(i)(i)
すると，{x(i)}N は，元の粒子 xk|k−1 の複製を βk にほぼ比例する個数だけ含んでおり，重k|k i=1
みづけをしなくてもフィルタ分布 p(xk |y 1:k ) を表現していることになる．このように N 回復
元抽出によって新たに N 個の粒子の集合を得る操作をリサンプリングと呼ぶ．
アルゴリズムをまとめると，以下のようになる：(i)t = 0 における状態 x0 の分布 p(x0 ) にしたがう乱数 x0|0 ∼ p(x0 ) を N 個生成する．
（1）（2）毎時間ステップ k (k = 1, . . . , K) において以下を実行する．
（a）予測• 各 i (i = 1, . . . , N ) について p(v k ) にしたがう乱数 v (i)k ∼ p(v k ) を N 個生成する．
• 各粒子 i (i = 1, . . . , N ) について

統計数理

342

第 66 巻 第 2 号 2018

(i)

(i)

(i)

xk|k−1 = f k (xk−1|k−1 ) + v k(i)
により，x(i)k−1|k−1 から xk|k−1 を生成する．
（b）フィルタリング〈尤度・重みの計算〉 各粒子 i について，尤度 p(y k |x(i)k|k−1 ) を計算し，重み
(i)

(i)

N


βk = p(y k |xk|k−1 )(i)p(y k |xk|k−1 )

i=1

を求める．
(N )
(i)
(i)
〈リサンプリング〉 粒子の集合 {x(1)k|k−1 , . . . , xk|k−1 } から，各粒子 xk|k−1 が βk の確率で(N )抽出されるように復元抽出を N 回繰り返す「リサンプリング」を行い，{x(1)k|k , . . . , xk|k }を生成する．

このように，尤度に比例する重みでリサンプリングを行い，フィルタ分布を計算するのが PFの基本的な形式（Gordon et al., 1993; Kitagawa, 1996）である．このような形式のアルゴリズ
ムを特にブートストラップフィルタ，あるいはモンテカルロフィルタと呼ぶ場合もある．
3.

並列アルゴリズムPF では，必要な粒子数 N が xk の次元に対して指数関数的に増大する．そこで，多数の粒子を高速で処理する手段として，並列計算機の使用が考えられる．粒子フィルタのアルゴリズの予測の手続きは，各粒子の処理が独立しているため，容易に並列化できる．
ムにおいて，（a）
しかし，（b）のフィルタリングを行う際に必要となるリサンプリングの手続きが，並列化の際に問題となる．
図 1 は，通常の粒子フィルタのリサンプリングを並列計算機上で行った場合の概念図である．ここでは，互いにメモリを共有しない複数のプロセスによって並列に処理が行われるマル
チプロセス型の並列計算を考える．リサンプリングの手続きでは，重み βk(i) の小さい（観測と合わない）
粒子は破棄され，重み βk(i) に応じて粒子の複製が生成されるため，個々のプロセスが保持する粒子数にばらつきが生じることになる．次の予測の手続きに進んだときに高い並列化
効率を得るには，粒子数の増えたプロセスから粒子数の減ったプロセスに粒子を移動し，各プロセスの保持する粒子の数を均等に再配分する操作が必要になる．粒子数を決める重み βk(i) の値にはランダム性があることから，各プロセスの粒子の増減数はランダムに決まる．したがって，粒子の再配分に伴うプロセス間通信には規則性がなく，並列処理が難しい．
そこで，高い並列化効率を実現する方法として，粒子をグループに分割してリサンプリングを行い，グループごとに異なる重みを持たせる方法が提案されている（Bolić et al., 2005）．単にグループに分けただけでは，各グループに割り当てられた個数の粒子のみで推定を行うのと同等になってしまうため，全体で高い精度の推定結果を得るには，グループ間の情報交換が必要になる．情報交換は，グループ間で粒子を部分的に交換する方法（Balasingam et al., 2011;
や，グループ間の相互作用を考える方法（Hlinka et al., 2013; Savic et al., 2014）
Bai et al., 2016）なども研究が進められているが，P3 では，並列化効率を重視し，次に述べる 2 種類の手法を実装している．
3.1 階層的リサンプリング
P3 で用意しているリサンプリング法の 1 つ目は，階層的リサンプリング（Nakano, 2010）
である．この方法は，図 2 に示すように，まず各プロセス内でリサンプリングを行い，次に各プ

P3 ：Python による並列計算機用粒子フィルタライブラリ343

図 1．通常のリサンプリング手法を並列計算機上で行った場合．

図 2．階層的リサンプリングの概念図（Nakano, 2010）．

ロセスに属する粒子の集合を超粒子と見なして，超粒子をリサンプリングするというものである．
階層的リサンプリングを行うには，まず各プロセス μ に割り当てられた粒子の集合 Gμ の重
みは
(3.1)

(μ)
Ωk =

λ


([μ−1]λ+ν)

βk

ν=1

となる．但し，λ は各プロセスに属する粒子の数で，ここでは全プロセスに均等に λ 個の粒子が割り当てられているものとする．したがって，プロセス数を Λ として，

344

統計数理

(3.2)

第 66 巻 第 2 号 2018

λ = N/Λ

である．各粒子の重み βk([μ−1]λ+ν) を集合 Gμ の重み Ω(μ)
で割ると，集合 Gμ 内での各粒子の重
k
みは
(3.3)

(ν)

ωk

([μ−1]λ+ν)

= βk

(μ)

/Ωk

となる．次に，各集合 Gμ の中で，重み {ωk(ν) } を用いて粒子のリサンプリング (ローカル・リ
サンプリング) を行う．ローカル・リサンプリングの手続きは，各プロセスで閉じており，異なるプロセスのローカル・リサンプリングは並列に実行できる．最後に，各集合 Gμ を超粒子
（メタ・リサンプリング）
する．これにより，観と見なし，重み {Ω(μ)
k } を用いてリサンプリング
測と合わない粒子しか保持していない集合は破棄され，観測と合う粒子の集合に置き換えられる．P3 では，ローカル・リサンプリングの手続きを local_resampling() という関数で，メタ・リサンプリングの手続きを meta_resampling() という関数で提供している．
なお，各 μ に対する Ω(μ)
の値のばらつきが大きくない場合は，メタ・リサンプリングを行k(μ)
わず，Ωk で重み付けしたままの方がフィルタ分布 p(xk |y 1:k ) をよりよく表現できると考えらのばらつきを評価するためにエントロピー
れる．そこで P3 では，Ω(μ)k(3.4)

SΩ,k = −

Λ


(μ)

(μ)

Ωk log Ωk

μ=1

から
(3.5)

Λeﬀ,k = eSΩ,kという量を計算し，Λeﬀ,k がある閾値より小さくなった場合のみに，メタ・リサンプリングの手続きが実行される．
3.2 Alternately lattice-pattern switching（ALPS）法3
P で提供するリサンプリング法の 2 つ目は，alternately lattice-pattern switching（ALPS）法である．この方法では，複数のプロセスに割り当てられた粒（Nakano and Higuchi, 2010, 2012）
子の集合をまとめてグループにし，各グループの中でローカル・リサンプリングと同じ操作を適用する．グループ分けを行う際には，まず，2m × 2n 個のプロセスを図 3 に示すような 2 次（プロセス）
元トーラス状のグラフの各ノードに割り当てる．そして，図 3 左のように，ノード

図 3．Alternate lattice-pattern switching 法の概念図（Nakano and Higuchi, 2010）．

P3 ：Python による並列計算機用粒子フィルタライブラリ345

図 4．P3 の構成の概略図．

を 2 × 2 個のグループに分割する．リサンプリングは，各グループで並列に実行される．これだけでは，単なるローカル・リサンプリングと大差ないが，ALPS では，次のステップで，図3 右に示すようにグループ分けのパターンを変え，その上でグループごとのリサンプリングを行う．ステップごとにグループ分けパターンを変えることで，観測に合う粒子の情報が全プロ
セスに行き渡るようにし，推定性能を向上させる．
4.

P3 の構成カルマンフィルタを適用できる線型状態空間モデルの場合，状態遷移行列，観測行列などの行列を与えれば記述できる．しかし，PF の扱う一般の非線型状態空間モデルの定義の仕方は，個々の問題によって様々であり，非線型微分方程式が使われる場合もあるし，非ガウスの確率分布が含まれることも考えられる．このような様々な場合に対応できるようにするため，P3 で（2.2）の状態空間モデルの内は，既存の Python ライブラリとの互換性には配慮せず，式（2.1），（2.2）の計算を P3 のモ
容をユーザが自由に定義できるように設計されている．但し，式（2.1），
ジュールから実行できるようにするため，ユーザは，状態ベクトル，状態空間モデルの定義を
system.py という名前のファイルに記述しておく必要がある．状態ベクトル，状態空間モデルの定義は，それぞれ class State，class System という名前で所定の形式に従って記述する．
PF を実行するために必要な関数は，ファイル pcubed.py 中で定義されている class Filterのメンバ関数の形で用意されている．並列計算やリサンプリングの処理は，class Filter の中に記述されており，ユーザは class State，class System を記述さえすれば，面倒なプログ
ラミング作業を行わずとも class Filter 中の関数を使って状態推定の計算が実行できる．図
4 に，P3 の構成を示す．
なお，従来の逐次ベイズ推定のための計算ライブラリでは，状態ベクトルを 1 つの配列にまとめた上で計算を行う実装が多かった．しかし，高次元の状態空間モデルにおいては，状態変数のそれぞれが別々の意味を持つことが少なくない．このような場合，異なる意味を持つ状態
変数は，プログラム上でも異なる変数名で扱った方が，可読性が向上すると考えられる．特に，シミュレーションモデルからシステムモデルを構成するデータ同化においては，元となるシミュレーションプログラムの変数名が流用できるため，システムモデルのプログラムを書く際にも都合がよい．そのため P3 では，状態ベクトルを配列ではなく，class State というクラス（構造体）の形で保持する設計になっている．
以下では，実際に P3 を使うために必要な情報として，class State と class System の定義

346

統計数理

第 66 巻 第 2 号 2018

表 1．class System で定義すべき関数．

の仕方を説明し，class Filter の内容について説明した後，メインプログラムの書き方の概要を説明する．
4.1 class State の定義
P3 の class Filter を使う際には，予め system.py の中で class State と class System を（2.1）に出てくる状態ベクトル xk の全要
定義する必要がある．このうち，class State は，式素（つまり全状態変数）をまとめたもので，推定するべき状態変数は，すべて class State のメンバ変数として列挙する．
各メンバ変数は，numpy.ndarray とし，0 で初期化するようにする．例えば，a0,k , a1,k , . . . ,a4,k と b0,k , b1,k , b2,k という 8 つの変数をまとめて状態ベクトル xk として扱いたいときは，
class State:def __init__(self):
self.a = np.zeros((5))
self.b = np.zeros((3))のように定義する．
4.2 class System の定義
class State は，状態ベクトルの定義を与えるだけで，式（2.1）のようなシステムモデルの定義は，class System で行う．class System の中では，pcubed.py から呼ばれる関数を所定の名前で定義する必要がある．表 1 が，class System で定義すべき関数である．以下に具体的な定義の仕方を述べる．
system_config(cls)この関数は，クラスメソッドとして定義され，事前に実行しておくべき前処理や，パラメー
タ値の設定をここに記述する．前処理等が必要ない場合も，system_config() をダミー関数にするなど，何らかの形で定義する必要がある．
init_state(self, ydata, x)各粒子の初期値を設定する．引数 x は class State 型の変数であり，x のメンバ変数に状態
変数の初期値が代入されるように関数 init_state() を定義する必要がある．PF において，各粒子の初期値は，初期状態の確率分布 p(x0 ) にしたがう乱数で与えるので，基本的に p(x0 ) にしたがう乱数が x のメンバ変数に代入されるように定義すればよい．
ydata は System.nobs の長さを持つ numpy.ndarray クラスのオブジェクトで，初期化に用いる観測データをここに与えることができる．これは，観測可能な状態変数をデータから与えるというようなことを想定している．つまり，p(x0 ) ではなく，p(x0 |y 0 ) から各粒子の初期値を

P3 ：Python による並列計算機用粒子フィルタライブラリ347

生成できるようになっている．
step_system_model(self, x)
1 ステップ分の粒子の時間発展を計算する．データ同化を行う場合，シミュレーションモデ
ルの 1 ステップ分をここで定義することになる．
x は入力と出力を兼ねた class State の変数で，粒子 i の時刻 tk−1 における状態変数の値(i)(i)
xk−1|k−1 を x として与えると，f k (xk−1|k−1 ) の値が State クラスのオブジェクトの形で x に代入されるようにする．
add_noise(self, x)粒子にシステムノイズ v (i)k を付加する．粒子の持つ状態変数の値を State クラスの形で x に入力すると，それにノイズが付加された値が x に代入されるように定義する．add_noise() を
step_system() と組み合わせることで
(i)

(i)

(i)

xk|k−1 = f k (xk−1|k−1 ) + v kの計算が実行できる．
log_likelihood(self, x, ydata)(i)
観測データ y k が与えられた下での粒子 x(i)k|k−1 の対数尤度 p(y k |xk|k−1 ) を計算する．入力として，粒子の状態変数の値を State クラスの形で x に，観測値を長さ System.nobs のnumpy.ndarray クラスの形で ydata に与えるようにする．計算された対数尤度は，実数値（ス
カラー）の返り値として返すようにする．
4.3 class Filter の内容
class Filter は，表 2 に示す関数で構成されており，メインプログラムから，表 2 の関数を呼ぶことで粒子フィルタの計算が実現できる．以下で，各関数の役割を述べる．
init_ensemble(nptot, yinit)
を入力し，yinit には初期化のN 個の粒子を初期化する．引数 nptot には粒子数 N（整数値）
ための観測データを入力する．yinit の型は numpy.ndarray で System.nobs の長さを持つ．具体的な初期化の方法は，system.py に記述する必要がある．
finalize()
終了処理をする．class Filter を使い終わった後に呼ぶ．
step_ensemble()
各粒子を 1 ステップ進める．これは式（2.7）の f k (x(i)
k−1|k−1 ) の部分の計算に相当する．f k の
表 2．class Filter で用意されている関数．

統計数理

348

第 66 巻 第 2 号 2018

定義は system.py に記述する．
add_noise_to_ensemble()
各粒子にシステムノイズを加える．これは式（2.7）で v (i)k を足す部分に相当する．
ensemble_mean()
粒子が保持する状態 x(i)k|k の平均1  (i)xk|k
NNx̄k|k =

i=1を求める．これにより PF による xk の推定値が得られる．
dist_reweight(ydata)
観測データ ydata を参照して，式（2.8）の重み βk(i) を計算する．ydata は numpy.ndarray で
System.nobs の長さを持つ．
local_resample()
dist_reweight で求めた重みにしたがい，各プロセス内で粒子のローカル・リサンプリングを行う．
meta_resample()
各プロセスに割り当てられた粒子の集合をまとめて超粒子と見なし，dist_reweight で求めた重みから計算される超粒子の重みにしたがって，メタ・リサンプリングを行う．但し，この関のばらつきを評価するために，式（3.5）の Λeﬀ を計算し，Λeﬀ ≥ 0.5
数では，超粒子の重み Ω(μ)
kの時は何もしない（つまりメタ・リサンプリングを行わない）
ようになっている．
regional_resample()
ALPS 法によるリサンプリングを実行する．ALPS のグループ分けパターンは，この関数を呼ぶ度に自動的に切り替えられる．
4.4 プログラムの書き方のまとめ
class Filter を用いて，PF の予測ステップを行うには，メインプログラム中でFilter.step_ensemble()
Filter.add_noise_to_ensemble()のように記述する．フィルタリングについては，階層的リサンプリングを用いる場合，
Filter.dist_reweight( ydata )Filter.local_resample()
Filter.meta_resample()のように記述すればよく，ALPS を用いる場合，
Filter.dist_reweight( ydata )Filter.regional_resample()のように記述すればよい．
5.

おわりに

PF は，状態変数が数個程度の小規模な問題であれば容易に実装できるため，非線型・非ガウスの状態空間モデルにおいて広く活用されているが，問題の規模に対して，計算量が指数関P3 ：Python による並列計算機用粒子フィルタライブラリ349

数的に増大するという問題がある．P3 は，粒子フィルタの適用対象としては，比較的規模の大きい非線型の問題に対しても，並列計算機を利用して状態推定を実現できる手段を提供する．粒子フィルタを並列計算機で実行しようとすると，個々の粒子を並列処理するために並列
プログラミングの知識が必要になる上，リサンプリング時に粒子を再配分する処理の実装に労力が必要となるが，P3 では並列計算やリサンプリングの処理が抽象化されているため，ユー
ザは煩雑な処理を自分でプログラミングしなくても，並列化効率の高い PF アルゴリズムを利用でき，状態空間モデルの構築に専念することができると考えられる．P3 は無償で配布している．利用を希望される方は，次のウェブサイト（http://daweb.ism.ac.jp/support/software/に記載の要領で申し込みいただきたい．
P-cubed/P-cubed.html）
現在は PF のみを実装しているが，用途によっては別の手法を用いた方がよい場合もある．
例えば，数百次元以上の大規模な非線型状態空間モデルを扱う場合には，PF よりもアンサンが有効であるし，システムノイズがガウス的であブルカルマンフィルタ（Evensen, 1994, 2003）
る場合には，混合ガウスフィルタ（Stordal et al., 2011）などが有効であると考えられる．今後
は，このような他の有用な非線型逐次ベイズ推定手法を追加するなどして，機能の充実を図る
予定である．
謝辞

P3 のプログラム開発にあたっては，科学研究費補助金基盤研究 B（課題番号: 26280010）の助成を受けた．ここに感謝の意を表する．

参

考

文

献

Bai, F., Gu, F., Hu, X. and Guo, S. (2016). Particle routing in distributed particle ﬁlters for large-scalespatial temporal systems, IEEE Transactions on Parallel and Distributed Systems, 27, 481–493.
Balasingam, B., Bolić, M., Djurić, P. M. and Míguez, J. (2011). Eﬃcient distributed resampling forparticle ﬁlters, Proceedings of IEEE International Conference of Acoustics, Speech, Signal Processing, 3772–3775.
Bengtsson, T., Bickel, P. and Li, B. (2008). Curse-of-dimensionality revisited: Collapse of the particle ﬁlter in very large scale systems, Probability and Statistics: Essays in Honors of David A.
Freedman, 2, 316–334, Institute of Mathematical Statistics, Beachwood, Ohio.
Bolić, M., Djurić, P. M. and Hong, S. (2005). Resampling algorithms and architectures for distributed
particle ﬁlters, IEEE Transactions on Signal Processing, 53, 2442–2450.
Daum, F. and Huang, J. (2003). Curse of dimensionality and particle ﬁlters, Proceedings of IEEE
Aerospace Conference, 4, 1979–1993.
Doucet, A., de Freitas, N. and Gordon, N. (eds.) (2001). Sequential Monte Carlo Methods in Practice,Springer-Verlag, 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.
Evensen, G. (2003). The ensemble Kalman ﬁlter: theoretical formulation and practical implementation,
Ocean Dynamics, 53, 343–367.
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.
Hlinka, J., Hlawatsch, F. and Djurić, P. M. (2013). Distributed particle ﬁltering in agent networks,
IEEE Signal Processing Magazine, 30, 61–81.
Kitagawa, G. (1993). Monte Carlo ﬁltering and smoothing method for non-Gaussian nonlinear state
space model, Reserch Memo., No.462, The Institute of Statistical Mathematics, Tokyo.

350

統計数理

第 66 巻 第 2 号 2018

Kitagawa, G. (1996). Monte Carlo ﬁlter and smoother for non-Gaussian nonlinear state space models,
Journal of Computational and Graphical Statistics, 5, 1–25.
Labbe, R. (2015). Kalman and Bayesian Filters in Python, http://ﬁlterpy.readthedocs.io（2018 年 2 月
17 日閲覧）.
Nakamura, K., Yoshida, R., Nagasaki, M., Miyano, S. and Higuchi, T. (2009). Parameter estimationof in silico biological pathways with particle ﬁltering towards peta-scale computing, Proceedingsof Pacific Symposium on Biocomputing, 14, 227–238.
Nakano, S. (2010). Population-based quasi-Bayesian algorithm for high-dimensional sequential problems and hierarchization of it for distributed computing environments, Proceedings of 2010 IEEE
World Congress on Computational Intelligence.
Nakano, S. and Higuchi, T. (2010). A dynamic grouping strategy for implementation of the particle
ﬁlter on a massively parallel computer, Proceedings of 13th International Conference on Information Fusion.
Nakano, S. and Higuchi, T. (2012). Weight adjustment of the particle ﬁlter on distributed computing
system, Proceedings of 15th International Conference on Information Fusion, 2480–2485.
Savic, V., Wymeersch, H. and Zozo, S. (2014). Belief consensus algorithms for fast distributed target
tracking in wireless sensor networks, Signal Processing, 95, 149–160.
Snyder, C., Bengtsson, T., Bickel, P. and Anderson, J. (2008). Obstacles to high-dimensional particle
ﬁltering, Monthly Weather Review, 136, 4629–4640.
Stordal, A. S., Karlsen, H. A., Nævdal, G., Skaug, H. J. and Vallès, B. (2011). Bridging the ensemble Kalman ﬁlter and particle ﬁlters: The adaptive Gaussian mixture ﬁlter, Computational
Geosciences, 15, 293–305.
van Leeuwen, P. J. (2009). Particle ﬁltering in geophysical systems, Monthly Weather Review, 137,
4089–4114.

Proceedings of the Institute of Statistical Mathematics Vol. 66, No. 2, 339–351 (2018)351

P3 : Python Parallelized Particle Filter Library
Shin’ya Nakano1,2 , Yuya Ariyoshi1,3 and Tomoyuki Higuchi1,2
1 The Institute of Statistical Mathematics

2 School of Multidisciplinary Science, SOKENDAI

3 Now at Faculty of Engineering, Nippon Bunri University

Particle ﬁlter (PF) is a class of state-estimation techniques based on Monte Carlocomputation that use a large number of particles. Because PF is applicable even to nonlinear and/or non-Gaussian problems, it is used for a variety of purposes. One serious
problem of PF is its computational time, which is exponential in the degrees of freedomof the state vector. Parallel computing is an eﬀective way to decrease computationaltime, but this approach requires skills in parallel programming. Even for experienced
users, it is challenging to achieve high computational eﬃciency in PF computation because the PF algorithm contains a procedure diﬃcult to parallelize. We developed a
Python library named P3 (Python Parallelized Particle Filter Library), that enables us toreadily use parallel-ready PF algorithms with high parallel eﬃciency. In this paper, we
describe the parallelized PF algorithms available in P3 , as well as explaining the designand characteristics of the library.

Key words: Particle ﬁlter, parallel computing, Python.