(57)【要約】【課題】精度よく、推定対象の状態を観測値から推定することができる。【解決手段】フィルタリング部３０が、入力部１０により受け付けた情報に基づいて、観測値ｙｔのノイズが正規分布に従う場合に、ガウス和リサンプリングフィルタを適用し、観測値ｙｔのノイズが正規分布に従わない場合にモンテカルロ粒子フィルタを適用する。【選択図】図２
2016218674.tif 000002
【特許請求の範囲】
【請求項１】
入力された推定対象の各時刻ｔの観測値に基づいて、前記推定対象の各時刻ｔの状態を推定する状態推定装置であって、
各時刻ｔについて、時刻ｔ−１についてのノイズと、ノイズの分布とに基づいて、ノイズを表す準粒子をＭ個サンプリングし、状態のフィルタ推定量のｉ番目の実現値を表す粒子の各々に対し、サンプリングしたＭ個の準粒子を用いて、時刻ｔ−１のｉ番目の前記粒子から予測される時刻ｔの状態の予測値をＭ個計算し、Ｍ個の時刻ｔの状態の予測値の各々に対する確率密度の分布に従って、時刻ｔのｉ番目の前記粒子をサンプリングし、サンプリングされた時刻ｔのｉ番目の前記粒子の各々と、時刻ｔの前記観測値に基づいて計算されるｉ番目の前記粒子の各々に対する重みとに基づいて、時刻ｔの状態推定量を計算するフィルタリング部
を含む状態推定装置。
【請求項２】
入力された推定対象の各時刻ｔの観測値に基づいて、前記推定対象の各時刻ｔの状態を推定する状態推定装置であって、
各時刻ｔについて、
時刻ｔ−１についてのノイズと、ノイズの分布とに基づいて、ノイズの実現値をＭ個サンプリングし、状態のｉ番目のガウシアンコンポーネントの各々に対し、サンプリングしたＭ個のノイズの実現値を用いて、時刻ｔ−１のｉ番目のガウシアンコンポーネントによるフィルタ推定量から予測される時刻ｔの状態の予測値をＭ個計算し、前記状態のｉ番目のガウシアンコンポーネントの各々に対して計算したＭ個の状態の予測値に対応する、前記状態の予測値、及び時刻ｔ−１のｉ番目のガウシアンコンポーネントの誤差分散から計算される時刻ｔの誤差分散の組み合わせを有するＭ個の準ガウシアンコンポーネントから、時刻ｔの前記観測値に基づいて前記状態のｉ番目のガウシアンコンポーネントの各々に対して計算される前記Ｍ個の準ガウシアンコンポーネントの各々に対する重みに従って、ｉ番目のガウシアンコンポーネントの各々をサンプリングし、
サンプリングされたｉ番目のガウシアンコンポーネントの各々と、サンプリングされたｉ番目のガウシアンコンポーネントの各々に対する前記重みとを用いて表される、時刻ｔの状態の各々に対する確率密度の分布を表すフィルタ分布に従って、時刻ｔの状態推定量を計算するフィルタリング部
を含む状態推定装置。
【請求項３】
前記フィルタリング部は、前記観測値のノイズが正規分布に従う場合に、
各時刻ｔについて、前記フィルタ分布に従って、時刻ｔの状態推定量を計算し、
前記観測値のノイズが正規分布に従わない場合に、各時刻ｔについて、時刻ｔ−１についてのノイズと、ノイズの分布とに基づいて、ノイズを表す準粒子をＭ個サンプリングし、状態のフィルタ推定量のｉ番目の実現値を表す粒子の各々に対し、サンプリングしたＭ個の準粒子を用いて、時刻ｔ−１のｉ番目の前記粒子から予測される時刻ｔの状態の予測値をＭ個計算し、Ｍ個の時刻ｔの状態の予測値の各々に対する確率密度の分布に従って、時刻ｔのｉ番目の前記粒子をサンプリングし、サンプリングされた時刻ｔのｉ番目の前記粒子の各々と、時刻ｔの前記観測値に基づいて計算されるｉ番目の前記粒子の各々に対する重みとに基づいて、時刻ｔの状態推定量を計算する請求項２に記載の状態推定装置。
【請求項４】
入力された推定対象の各時刻ｔの観測値に基づいて、前記推定対象の各時刻ｔの状態を推定する状態推定装置における状態推定方法であって、
フィルタリング部が、各時刻ｔについて、時刻ｔ−１についてのノイズと、ノイズの分布とに基づいて、ノイズを表す準粒子をＭ個サンプリングし、状態のフィルタ推定量のｉ番目の実現値を表す粒子の各々に対し、サンプリングしたＭ個の準粒子を用いて、時刻ｔ−１のｉ番目の前記粒子から予測される時刻ｔの状態の予測値をＭ個計算し、Ｍ個の時刻ｔの状態の予測値の各々に対する確率密度の分布に従って、時刻ｔのｉ番目の前記粒子をサンプリングし、サンプリングされた時刻ｔのｉ番目の前記粒子の各々と、時刻ｔの前記観測値に基づいて計算されるｉ番目の前記粒子の各々に対する重みとに基づいて、時刻ｔの状態推定量を計算するステップ
を含む状態推定方法。
【請求項５】
入力された推定対象の各時刻ｔの観測値に基づいて、前記推定対象の各時刻ｔの状態を推定する状態推定装置における状態推定方法であって、
フィルタリング部が、各時刻ｔについて、
時刻ｔ−１についてのノイズと、ノイズの分布とに基づいて、ノイズの実現値をＭ個サンプリングし、状態のｉ番目のガウシアンコンポーネントの各々に対し、サンプリングしたＭ個のノイズの実現値を用いて、時刻ｔ−１のｉ番目のガウシアンコンポーネントによるフィルタ推定量から予測される時刻ｔの状態の予測値をＭ個計算し、前記状態のｉ番目のガウシアンコンポーネントの各々に対して計算したＭ個の状態の予測値に対応する、前記状態の予測値、及び時刻ｔ−１のｉ番目のガウシアンコンポーネントの誤差分散から計算される時刻ｔの誤差分散の組み合わせを有するＭ個の準ガウシアンコンポーネントから、時刻ｔの前記観測値に基づいて前記状態のｉ番目のガウシアンコンポーネントの各々に対して計算される前記Ｍ個の準ガウシアンコンポーネントの各々に対する重みに従って、ｉ番目のガウシアンコンポーネントの各々をサンプリングし、
サンプリングされたｉ番目のガウシアンコンポーネントの各々と、サンプリングされたｉ番目のガウシアンコンポーネントの各々に対する前記重みとを用いて表される、時刻ｔの状態の各々に対する確率密度の分布を表すフィルタ分布に従って、時刻ｔの状態推定量を計算するステップ
を含む状態推定方法。
【請求項６】
コンピュータを、請求項１に記載の状態推定装置の各部として機能させるためのプログラム。
【請求項７】
コンピュータを、請求項２又は請求項３に記載の状態推定装置の各部として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、状態推定装置、方法、及びプログラムに係り、特に、推定対象の状態を観測値から推定する状態推定装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
非特許文献１では逐次ベイズフィルタを多数の粒子の使用により実現するMonte Carlo filter（以下、ＭＣＦと称する）を提案している。非線形非ガウス型の状態空間モデルに適用可能である。
【０００３】
また、非特許文献２には、ガウス和フィルタ（Gaussian sum filter、 以下、ＧＳＦと称する）が提案されている。ＧＳＦは線形ガウス型の状態空間モデルに対してのみ適用可能である。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】北川源四郎. “モンテカルロ・フィルタおよび平滑化について”, 統計数理, 44巻, 第1号, 31-48, 1996.
【非特許文献２】D. L. Alspach and H. W. Sorenson. “Nonlinear Bayesian Estimation Using Gaussian Sum Approximations”,Vol. AC-17, No. 4, pp. 439-448, IEEE Trans. on Automatic Control, 1972.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、非特許文献１の手法では、状態のフィルタ推定量の分布(フィルタ分布と呼ぶ) を表現する際に、状態の予測推定量の分布(予測分布と呼ぶ) からサンプリングされた粒子を使用するため、フィルタ分布と予測分布に大きな差がある場合、推定精度が下がってしまうという問題がある。
【０００６】
また、非特許文献２の手法では、予測分布のガウス和近似を構成する各ガウシアンコンポーネントの分散値が小さい状況でのみ、その後に続くフィルタ分布のガウス和近似の精度が理論的に保証されるが、その様な状況は一般的ではない。
【０００７】
本発明は、上記問題点を解決するために成されたものであり、精度よく、推定対象の状態を観測値から推定することができる状態推定装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記目的を達成するために、第１の発明に係る状態推定装置は、入力された推定対象の各時刻ｔの観測値に基づいて、前記推定対象の各時刻ｔの状態を推定する状態推定装置であって、各時刻ｔについて、時刻ｔ−１についてのノイズと、ノイズの分布とに基づいて、ノイズを表す準粒子をＭ個サンプリングし、状態のフィルタ推定量のｉ番目の実現値を表す粒子の各々に対し、サンプリングしたＭ個の準粒子を用いて、時刻ｔ−１のｉ番目の前記粒子から予測される時刻ｔの状態の予測値をＭ個計算し、Ｍ個の時刻ｔの状態の予測値の各々に対する確率密度の分布に従って、時刻ｔのｉ番目の前記粒子をサンプリングし、サンプリングされた時刻ｔのｉ番目の前記粒子の各々と、時刻ｔの前記観測値に基づいて計算されるｉ番目の前記粒子の各々に対する重みとに基づいて、時刻ｔの状態推定量を計算するフィルタリング部を含んで構成されている。
【０００９】
また、第２の発明に係る状態推定装置は、入力された推定対象の各時刻ｔの観測値に基づいて、前記推定対象の各時刻ｔの状態を推定する状態推定装置であって、各時刻ｔについて、時刻ｔ−１についてのノイズと、ノイズの分布とに基づいて、ノイズの実現値をＭ個サンプリングし、状態のｉ番目のガウシアンコンポーネントの各々に対し、サンプリングしたＭ個のノイズの実現値を用いて、時刻ｔ−１のｉ番目のガウシアンコンポーネントによるフィルタ推定量から予測される時刻ｔの状態の予測値をＭ個計算し、前記状態のｉ番目のガウシアンコンポーネントの各々に対して計算したＭ個の状態の予測値に対応する、前記状態の予測値、及び時刻ｔ−１のｉ番目のガウシアンコンポーネントの誤差分散から計算される時刻ｔの誤差分散の組み合わせを有するＭ個の準ガウシアンコンポーネントから、時刻ｔの前記観測値に基づいて前記状態のｉ番目のガウシアンコンポーネントの各々に対して計算される前記Ｍ個の準ガウシアンコンポーネントの各々に対する重みに従って、ｉ番目のガウシアンコンポーネントの各々をサンプリングし、サンプリングされたｉ番目のガウシアンコンポーネントの各々と、サンプリングされたｉ番目のガウシアンコンポーネントの各々に対する前記重みとを用いて表される、時刻ｔの状態の各々に対する確率密度の分布を表すフィルタ分布に従って、時刻ｔの状態推定量を計算するフィルタリング部を含んで構成されている。
【００１０】
また、第２の発明に係る状態推定装置において、前記フィルタリング部は、前記観測値のノイズが正規分布に従う場合に、各時刻ｔについて、前記フィルタ分布に従って、時刻ｔの状態推定量を計算し、前記観測値のノイズが正規分布に従わない場合に、各時刻ｔについて、時刻ｔ−１についてのノイズと、ノイズの分布とに基づいて、ノイズを表す準粒子をＭ個サンプリングし、状態のフィルタ推定量のｉ番目の実現値を表す粒子の各々に対し、サンプリングしたＭ個の準粒子を用いて、時刻ｔ−１のｉ番目の前記粒子から予測される時刻ｔの状態の予測値をＭ個計算し、Ｍ個の時刻ｔの状態の予測値の各々に対する確率密度の分布に従って、時刻ｔのｉ番目の前記粒子をサンプリングし、サンプリングされた時刻ｔのｉ番目の前記粒子の各々と、時刻ｔの前記観測値に基づいて計算されるｉ番目の前記粒子の各々に対する重みとに基づいて、時刻ｔの状態推定量を計算するようにしてもよい。
【００１１】
第３の発明に係る状態推定方法は、入力された推定対象の各時刻ｔの観測値に基づいて、前記推定対象の各時刻ｔの状態を推定する状態推定装置における状態推定方法であって、フィルタリング部が、各時刻ｔについて、時刻ｔ−１についてのノイズと、ノイズの分布とに基づいて、ノイズを表す準粒子をＭ個サンプリングし、状態のフィルタ推定量のｉ番目の実現値を表す粒子の各々に対し、サンプリングしたＭ個の準粒子を用いて、時刻ｔ−１のｉ番目の前記粒子から予測される時刻ｔの状態の予測値をＭ個計算し、Ｍ個の時刻ｔの状態の予測値の各々に対する確率密度の分布に従って、時刻ｔのｉ番目の前記粒子をサンプリングし、サンプリングされた時刻ｔのｉ番目の前記粒子の各々と、時刻ｔの前記観測値に基づいて計算されるｉ番目の前記粒子の各々に対する重みとに基づいて、時刻ｔの状態推定量を計算するステップを含んで実行することを特徴とする。
【００１２】
第４の発明に係る状態推定方法は、入力された推定対象の各時刻ｔの観測値に基づいて、前記推定対象の各時刻ｔの状態を推定する状態推定装置における状態推定方法であって、フィルタリング部が、各時刻ｔについて、時刻ｔ−１についてのノイズと、ノイズの分布とに基づいて、ノイズの実現値をＭ個サンプリングし、状態のｉ番目のガウシアンコンポーネントの各々に対し、サンプリングしたＭ個のノイズの実現値を用いて、時刻ｔ−１のｉ番目のガウシアンコンポーネントによるフィルタ推定量から予測される時刻ｔの状態の予測値をＭ個計算し、前記状態のｉ番目のガウシアンコンポーネントの各々に対して計算したＭ個の状態の予測値に対応する、前記状態の予測値、及び時刻ｔ−１のｉ番目のガウシアンコンポーネントの誤差分散から計算される時刻ｔの誤差分散の組み合わせを有するＭ個の準ガウシアンコンポーネントから、時刻ｔの前記観測値に基づいて前記状態のｉ番目のガウシアンコンポーネントの各々に対して計算される前記Ｍ個の準ガウシアンコンポーネントの各々に対する重みに従って、ｉ番目のガウシアンコンポーネントの各々をサンプリングし、サンプリングされたｉ番目のガウシアンコンポーネントの各々と、サンプリングされたｉ番目のガウシアンコンポーネントの各々に対する前記重みとを用いて表される、時刻ｔの状態の各々に対する確率密度の分布を表すフィルタ分布に従って、時刻ｔの状態推定量を計算するステップを含んで実行することを特徴とする。
【００１３】
第５の発明に係るプログラムは、コンピュータを、上記第１の発明に係る状態推定装置の各部として機能させるためのプログラムである。
【００１４】
第６の発明に係るプログラムは、コンピュータを、上記第２の発明に係る状態推定装置の各部として機能させるためのプログラムである。
【発明の効果】
【００１５】
本発明の状態推定装置、方法、及びプログラムによれば、精度よく、推定対象の状態を観測値から推定することができる、という効果が得られる。
【図面の簡単な説明】
【００１６】
【図１】時刻ｔの観測データと状態推定値のデータ構造を示す図である。
【図２】本発明の実施の形態に係る状態推定装置の構成を示すブロック図である。
【図３】本発明の実施の形態に係る状態推定装置における状態推定処理ルーチンを示すフローチャートである。
【発明を実施するための形態】
【００１７】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００１８】
＜本発明の実施の形態に係る原理＞
【００１９】
まず、本発明の実施の形態における原理を説明する。
【００２０】
本発明の実施の形態では、解析対象の状態を観測系列から推定するフィルタリングの問題において、従来フィルタより状態推定精度が高いフィルタを２つ提案する。モンテカルロ粒子フィルタ（Monte Carlo particle filter：ＭＣＰＦ) は非線形非ガウス型の状態空間モデルに対して適用可能であり、ガウス和リサンプリングフィルタ（Gaussian sumresampling filter：ＧＳＲＦ）は観測ノイズに正規性が仮定できる場合に適用可能である。
【００２１】
モンテカルロ粒子フィルタ（ＭＣＰＦ）は、非特許文献３に記載のsequential importance resampling（ＳＩＲ）の提案分布をＭＣＦで実現し、その分布からサンプリングされた粒子を用いて状態のフィルタ推定値を計算する。ＭＣＦによる提案分布の構築にはその非ガウス性が考慮されるため、従来法より優れたフィルタリング精度が期待できる。
【００２２】
［非特許文献３］：A. Doucet, N. De Freitas and N. Gordon. “Sequential Monte Carlo Methods in Practice”, Springer, 2001
【００２３】
一方、ガウス和リサンプリングフィルタ(ＧＳＲＦ) は、予測分布を多数の準ガウシアンコンポーネントの平均和で記述するため、従来のガウス和フィルタで用いられるガウス分布より小さい分散値のガウス分布の集合を使用することができ、フィルタ分布の推定精度が改善する。また、準ガウシアンコンポーネントの生成プロセスでガウス分布の多様性が確保されることも推定精度の向上につながる。
【００２４】
本発明の実施の形態では、計測器によって解析対象の時刻ｔにおける観測データが得られるものとする。この観測データ系列から解析対象の隠れた状態を推定することがフィルタリングの問題である。非線形非ガウス型の状態空間モデルは解析対象の時間発展とその観測プロセスを解析者がモデル化したものであるが、本発明の実施の形態では以下（１）式及び（２）式に従う非線形非ガウス型の状態空間モデルを想定する。
【００２５】
【数１】
2016218674.tif 000003
【００２６】
Ｘｔは時刻ｔにおける状態と呼ばれ、その時刻までの観測データの系列｛Ｙｔ：１｝＝｛Ｙｔ，Ｙｔ−１，・・・，Ｙ１｝が与えられた後の最も確からしい状態推定値を求める問題がフィルタリング問題である。また、一期前までの観測データの系列｛Ｙｔ−１：１｝が与えられた後の最も確からしい状態推定値を求める問題は予測問題と呼ばれる。Ｘｔ、及びＹｔはベクトル値の確率変数であり、次元はそれぞれＬＸ、ＬＹとする。ｆ(・)、ｈ(・)は状態とノイズΦｔ−１、Ψｔに関する非線形関数である。ここでΦｔ−１、Ψｔはそれぞれ任意の確率密度関数に従うＬＸ、ＬＹ次元の白色ノイズであるとする。
【００２７】
以降、簡単のため状態、観測値、及び各ノイズはスカラーの場合を想定する。つまり以下（３）式、及び（４）式に従う加法的なノイズの非線形状態空間モデルを想定し、その観測値ｙｔからフィルタリングによって状態ｘｔの推定値を求める。
【００２８】
【数２】
2016218674.tif 000004
【００２９】
次に、モンテカルロ粒子フィルタ（ＭＣＰＦ）及びガウス和リサンプリングフィルタ（ＧＳＲＦ）のそれぞれの原理の詳細について説明する。
【００３０】
＜モンテカルロ粒子フィルタ（ＭＣＰＦ）＞
モンテカルロ粒子フィルタ（ＭＣＰＦ）では、まず、sequential importance resampling（ＳＩＲ）の提案分布を以下（５）式に従って構築する。
【００３１】
【数３】
2016218674.tif 000005
【００３２】
上記（５）式により構築されるＳＩＲの提案分布は、フィルタ推定量ｘｔの各値に対する確率密度を表す確率分布である。ここで、ｘｔ−１ｉは時刻ｔ‐１の状態のフィルタ推定量のｉ番目の実現値であり、粒子とも呼ばれる。ここでは、全部でＮ個の粒子があるものとする。また、時刻ｔ−１の状態のフィルタ推定量から全部でＮ個の実現値ｘｔ−１ｉがサンプリングされているものとする。ｐ（ｙｔ｜ｘｔ）は状態ｘｔが与えられた後の観測値ｙｔの確率として計算する。ｐ（ｙｔ｜￣ｘｔｉ，ｋ）は状態ｘｔが予測値￣ｘｔｉ，ｋの場合の観測値ｙｔの確率であり、予測値￣ｘｔｉ，ｋは、以下の式で表わされる。
【００３３】
￣ｘｔｉ，ｋ＝ｆ（ｘｔ−１ｉ）＋φｔ−１ｋ、φｔ−１ｋ〜ｐ（φｔ−１）  ・・・（６）
【００３４】
φｔ−１ｋ〜ｐ（φｔ−１）は分布ｐ（φｔ−１）からサンプリングされた粒子であり、全部でＭ個の粒子があるとする。また、φｔ−１ｋ，（ｋ＝１，２，・・・，Ｍ）である。φｔ−１ｋのＭ個の粒子を準粒子と呼ぶ。δ（・）はディラックのデルタ関数であり、上記（５）式においてはδ（０）＝１となる関数である。
【００３５】
上記（５）式の提案分布からｉ番目の粒子を１つサンプリングし、それをｉ＝１〜Ｎまで繰り返す。サンプリングされた粒子を￣ｘｔｉ，ｊとすると、その重みｗｔｉは以下（７）式で計算される。
【００３６】
【数４】
2016218674.tif 000006
【００３７】
もし直前のステップでｉ＝１〜Ｎの粒子のリサンプリング処理が行われている場合はｗｔ−１ｉ＝１／Ｎ，（ｉ＝１，２，・・・，Ｎ）なので、（７）式で表される重みは、以下（８）式で表される。
【００３８】
【数５】
2016218674.tif 000007
【００３９】
この（５）式、又は（７）式若しくは（８）式をフィルタリングと呼ぶ。￣ｘｔｉ，ｊをｘｔｉと書き直すと、時刻ｔにおける状態のフィルタ推定値＾ｘｔは以下（９）式で計算される。
【００４０】
【数６】
2016218674.tif 000008
【００４１】
＜ガウス和リサンプリングフィルタ（ＧＳＲＦ）＞
次に、ガウス和リサンプリングフィルタ（ＧＳＲＦ）について説明する。ガウス和リサンプリングフィルタ（ＧＳＲＦ）では、状態のフィルタ分布を以下（１０）式に従って推定する。
【００４２】
【数７】
2016218674.tif 000009
【００４３】
ここで＾ｍｔｉ、及び＾ｐｔｉはｉ番目のガウシアンコンポーネントｘｉ，ｔの時刻ｔにおけるフィルタ推定値と誤差分散であり、以下の（１１）式及び（１２）式により表される状態空間モデルに基づいて算出される。
【００４４】
【数８】
2016218674.tif 000010
【００４５】
ここでＮ（＾ｍｔ−１ｉ，＾ｐｔ−１ｉ）及びＮ（￣ｍｔｉ，￣ｐｔｉ）はｘｉの時刻ｔ−１におけるフィルタ分布と時刻ｔにおける予測分布である。
【００４６】
上記（１０）式のβｔｉは以下の（１３）式で計算される。
【００４７】
【数９】
2016218674.tif 000011
【００４８】
ここで予測値￣ｍｔｉ，ｊ、及び誤差分散￣ｐｔｉ，ｊは、以下の式に従って計算される。
￣ｍｔｉ，ｊ＝ｆ（＾ｍｔ−１ｉ）＋φｔ−１ｉ，ｊ    ・・・（１４）
￣ｐｔｉ，ｊ＝（ｆ´ｉ，ｔ−１）２＾ｐｔ−１ｉ    ・・・（１５）
【００４９】
【数１０】
2016218674.tif 000012
【００５０】
【数１１】
2016218674.tif 000013
【００５１】
とした。また、αｔ−１ｉ＝１であり、ｒは定数である。ノイズφｔ−１ｉ，ｊ（ｊ＝１，２，・・・，Ｍ）は、分布ｐ（φｉ，ｔ−１）からサンプリングされたＭ個の実現値（サンプリング値）である。Ｍ個の実現値をサンプリングする処理によりΣｉ＝１ＮΣｊ＝１Ｍβｔｉ，ｊ＝１となり、βｔｉ，ｊをｉ番目のガウシアンコンポーネントに対するｊ番目の準ガウシアンコンポーネントの生成確率と見なすことができるようになる。次に、βｔｉ，ｊに従い、Ｎ×Ｍ個の準ガウシアンコンポーネントうち、Ｎ個の準ガウシアンコンポーネントをサンプリングする。また、サンプリングされた準ガウシアンコンポーネントを用いて、カルマンフィルタに従い、時刻ｔの状態のフィルタ推定量を計算する。
【００５２】
上記（１０）式及び（１３）式をフィルタリングと呼ぶ。サンプリングされた準ガウシアンコンポーネントのフィルタ推定量を＾ｍｔｉ、当該準ガウシアンコンポーネントに対応する重みβｔｉ，ｊをβｔｉと書き直すと、上記（１０）式から時刻ｔにおける状態のフィルタ推定値＾ｘｔは以下（１６）式に従って計算される。
【００５３】
【数１２】
2016218674.tif 000014
【００５４】
図１に時刻ｔの観測データと状態推定値のデータ構造を示す。
【００５５】
＜本発明の実施の形態に係る状態推定装置の構成＞
【００５６】
次に、本発明の実施の形態に係る状態推定装置の構成について説明する。図２に示すように、本発明の実施の形態に係る状態推定装置１００は、ＣＰＵと、ＲＡＭと、後述する状態推定処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することが出来る。この状態推定装置１００は、機能的には図１に示すように入力部１０と、演算部２０と、出力部５０とを備えている。
【００５７】
入力部１０は、観測機器（図示省略）によって得られた、推定対象の各時刻ｔの観測値ｙｔからなる観測データの入力を受け付ける。また、入力部１０は、観測値ｙｔのノイズが正規分布に従うか否かを示す情報の入力を受け付ける。
【００５８】
演算部２０は、フィルタリング部３０を含んで構成されている。
【００５９】
フィルタリング部３０は、入力部１０により受け付けた情報に基づいて、観測値ｙｔのノイズが正規分布に従う場合に、ガウス和リサンプリングフィルタを適用し、観測値ｙｔのノイズが正規分布に従わない場合にモンテカルロ粒子フィルタを適用する。
【００６０】
フィルタリング部３０が、モンテカルロ粒子フィルタを適用して時刻ｔの状態推定量＾ｘｔを計算する場合には、フィルタリング部３０は、各時刻ｔについて以下の処理を繰り返し行う。
【００６１】
フィルタリング部３０は、まず、時刻ｔ−１についてのノイズφｔ−１ｋと、ノイズの分布ｐ（φｔ−１）とに基づいて、ノイズφｔ−１ｋを表す準粒子をＭ個サンプリングする。次に、フィルタリング部３０は、状態のフィルタ推定量のｉ番目の実現値ｘｔｉを表す粒子の各々に対し、サンプリングしたＭ個の準粒子を用いて、上記（６）式に従って、時刻ｔ−１の当該ｉ番目の粒子から予測される時刻ｔの状態の予測値￣ｘｔｉ,ｋをＭ個計算し、上記（５）式に示す、Ｍ個の時刻ｔの状態の予測値￣ｘｔｉ，ｋの各々に対する確率密度の分布に従って、時刻ｔのｉ番目の粒子をサンプリングする。そして、フィルタリング部３０は、サンプリングされた時刻ｔのｉ番目の粒子の各々と、入力部１０により受け付けた観測値ｙｔに基づいて上記（７）式に従って計算されるｉ番目の粒子の各々に対する重みｗｔｉとに基づいて、上記（９）式に従って、時刻ｔの状態推定量＾ｘｔを計算する。
【００６２】
フィルタリング部３０が、ガウス和リサンプリングフィルタを適用して時刻ｔの状態推定量＾ｘｔを計算する場合には、フィルタリング部３０は、各時刻ｔについて、以下の処理を繰り返し行う。
【００６３】
フィルタリング部３０は、まず、時刻ｔ−１についてのノイズφｉ，ｔ−１と、ノイズの分布ｐ（φｉ，ｔ−１）とに基づいて、ノイズφｉ，ｔ−１の実現値φｔ−１ｉ，ｊをＭ個サンプリングする。次に、フィルタリング部３０は、上記（１４）式に従って、状態のｉ番目のガウシアンコンポーネントの各々に対し、サンプリングしたＭ個のノイズの実現値φｔ−１ｉ，ｊを用いて、時刻ｔ−１のｉ番目のガウシアンコンポーネントによるフィルタ推定量＾ｍｔ−１ｉから予測される時刻ｔの状態の予測値￣ｍｔｉ，ｊをＭ個計算する。フィルタリング部３０は、状態のｉ番目のガウシアンコンポーネントの各々に対して計算したＭ個の状態の予測値￣ｍｔｉ，ｊに対応する、状態の予測値￣ｍｔｉ，ｊ、及び時刻ｔ−１のｉ番目のガウシアンコンポーネントの誤差分散＾ｐｔ−１ｉから計算される時刻ｔの誤差分散￣ｐｔｉ，ｊの組み合わせを有するＮ×Ｍ個の準ガウシアンコンポーネントから、入力部１０により受け付けた時刻ｔの観測値ｙｔに基づいて状態のｉ番目のガウシアンコンポーネントの各々に対して上記（１３）式に従って計算されるＭ個の準ガウシアンコンポーネントの各々に対する重みβｔｉ，ｊに従って、Ｎ個の準ガウシアンコンポーネントを、ｉ番目のガウシアンコンポーネントの各々としてサンプリングする。そして、フィルタリング部３０は、サンプリングされたｉ番目のガウシアンコンポーネントの各々のフィルタ推定量＾ｍｔｉ、及び誤差分散＾ｐｔｉと、サンプリングされたｉ番目のガウシアンコンポーネントの各々に対する重みβｔｉとを用いて表される上記（１０）式のフィルタ分布に従って、時刻ｔの状態推定量＾ｘｔを計算する。上記（１０）式のフィルタ分布に従って、時刻ｔの状態推定量＾ｘｔを計算することでは、具体的には、上記（１６）式に従って、サンプリングされたｉ番目のガウシアンコンポーネントの各々から推定される時刻ｔの状態のフィルタ推定量＾ｍｔｉと、サンプリングされたｉ番目のガウシアンコンポーネントの各々に対する重みβｔｉ，ｊとに基づいて、時刻ｔの状態推定量＾ｘｔを計算する。
【００６４】
そして、フィルタリング部３０は、モンテカルロ粒子フィルタ又はガウス和リサンプリングフィルタにより計算した各時刻ｔの状態推定量＾ｘｔを出力部５０に出力する。
【００６５】
＜本発明の実施の形態に係る状態推定装置の作用＞
【００６６】
次に、本発明の実施の形態に係る状態推定装置１００の作用について説明する。入力部１０において推定対象の各時刻の観測値ｙｔからなる観測データ、及び観測値ｙｔのノイズが正規分布に従うか否かを示す情報の入力を受け付けると、状態推定装置１００は、図３に示す状態推定処理ルーチンを実行する。
【００６７】
まず、ステップＳ１００では、入力部１０において受け付けた観測データの観測値ｙｔを取得する。
【００６８】
次に、ステップＳ１０２では、入力部１０において受け付けた情報に基づいて、観測値ｙｔのノイズが正規分布に従うか否かを判定し、観測値ｙｔのノイズが正規分布に従わない場合にはステップＳ１０４へ移行し、観測値ｙｔのノイズが正規分布に従う場合にはステップＳ１１６へ移行する。
【００６９】
ステップＳ１０４では、ｉ＝１と初期化する。
【００７０】
ステップＳ１０６では、ｉ番目の実現値ｘｔｉを表す粒子に対する、時刻ｔ−１についてのノイズφｔ−１ｋと、ノイズの分布ｐ（φｔ−１）とに基づいて、ノイズφｔ−１ｋを表す準粒子をＭ個サンプリングする。
【００７１】
ステップＳ１０８では、状態のフィルタ推定量のｉ番目の実現値ｘｔｉを表す粒子に対し、ステップＳ１０６でサンプリングしたＭ個の準粒子を用いて、上記（６）式に従って、時刻ｔ−１の当該ｉ番目の粒子から予測される時刻ｔの状態の予測値￣ｘｔｉ,ｋをＭ個計算する。
【００７２】
ステップＳ１１０では、上記（５）式に示す、ステップＳ１０８で計算されたＭ個の時刻ｔの状態の予測値￣ｘｔｉ,ｋの各々に対する確率密度の分布に従って、時刻ｔのｉ番目の粒子をサンプリングする。
【００７３】
ステップＳ１１２では、ｉ＝Ｎであるか否かを判定し、ｉ＝Ｎでない場合にはステップＳ１１４へ移行し、ｉ＝Ｎである場合にはステップＳ１３２へ移行する。
【００７４】
ステップＳ１１４では、ｉ＝ｉ＋１として、ステップＳ１０６へ移行してステップＳ１０６〜ステップＳ１１２の処理を繰り返す。
【００７５】
ステップＳ１１６では、ステップＳ１００で取得した時刻ｔの観測値ｙｔに基づいて、上記（７）式に従って、ｉ番目の粒子の各々に対する重みｗｔｉを計算する。
【００７６】
ステップＳ１１７では、ステップＳ１１０でサンプリングされた時刻ｔのｉ番目の粒子の各々と、ステップＳ１１６で計算されたｉ番目の粒子の各々に対する重みとに基づいて、上記（９）式に従って時刻ｔの状態推定量＾ｘｔを計算する。
【００７７】
ステップＳ１１８では、ｉ＝１と初期化する。
【００７８】
ステップＳ１２０では、ｉ番目のガウシアンコンポーネントに対する、時刻ｔ−１についてのノイズφｔ−１ｋと、ノイズの分布ｐ（φｔ−１）とに基づいて、ノイズφｔ−１ｋを表す実現値φｔ−１ｉ，ｊをＭ個サンプリングする。
【００７９】
ステップＳ１２２では、上記（１４）式に従って、状態のｉ番目のガウシアンコンポーネントに対し、サンプリングしたＭ個のノイズの実現値φｔ−１ｉ，ｊを用いて、時刻ｔ−１のｉ番目のガウシアンコンポーネントによるフィルタ推定量＾ｍｔ−１ｉから予測される時刻ｔの状態の予測値￣ｍｔｉ，ｊをＭ個計算する。
【００８０】
ステップＳ１２４では、時刻ｔ−１のｉ番目のガウシアンコンポーネントによるフィルタ推定量＾ｍｔ−１ｉと誤差分散＾ｐｔ−１ｉとに基づいて、上記（１５）式に従って、Ｍ個の状態の予測値￣ｍｔｉ，ｊに対応する、時刻ｔの誤差分散￣ｐｔｉ，ｊを計算する。
【００８１】
ステップＳ１２６では、ｉ＝Ｎかを判定し、ｉ＝Ｎでない場合にはステップＳ１２８へ移行し、ｉ＝Ｎである場合にはステップＳ１３０へ移行する。
【００８２】
ステップＳ１２８では、ｉ＝ｉ＋１として、ステップＳ１２０へ移行してステップＳ１２０〜ステップＳ１２６の処理を繰り返す。
【００８３】
ステップＳ１３０では、上記（１３）式に従って、ステップＳ１２２で計算された状態の予測値￣ｍｔｉ，ｊ、及びステップＳ１２４で計算された誤差分散￣ｐｔｉ，ｊの組み合わせを有するＮ×Ｍ個の準ガウシアンコンポーネントの各々に対して、ステップＳ１００で取得した時刻ｔの観測データの観測値ｙｔに基づいて、当該準ガウシアンコンポーネントに対する重みβｔｉ，ｊを計算する。
【００８４】
ステップＳ１３２では、ステップＳ１３０で計算されたＮ×Ｍ個の準ガウシアンコンポーネントに対する重みβｔｉ，ｊに従って、Ｎ個の準ガウシアンコンポーネントを、ｉ番目のガウシアンコンポーネントの各々としてサンプリングする。
【００８５】
ステップＳ１３４では、ステップＳ１３２でサンプリングされたｉ番目のガウシアンコンポーネントの各々を用いて、カルマンフィルタに従い、時刻ｔの状態のフィルタ推定量＾ｍｔｉ、及び誤差分散＾ｐｔｉを計算し、計算されたｉ番目のガウシアンコンポーネントの各々のフィルタ推定量＾ｍｔｉ、及び誤差分散＾ｐｔｉと、ステップＳ１３０で計算された、ｉ番目のガウシアンコンポーネントに対応する準ガウシアンコンポーネントの各々に対する重みβｔｉとを用いて、上記（１６）式に従って、時刻ｔの状態推定量＾ｘｔを計算する。
【００８６】
ステップＳ１３６では、ステップＳ１１７又はステップＳ１３４で計算された状態推定量＾ｘｔを出力部５０に出力し、ステップＳ１００へ戻って次の時刻ｔの観測値ｙｔを読み込んで処理を繰り返す。
【００８７】
＜実験例＞
【００８８】
本発明の実施の形態におけるＭＣＰＦ及びＧＳＲＦの有効性を示す実験例について述べる。本実験では以下の状態空間モデルを扱う。
【００８９】
【数１３】
2016218674.tif 000015
【００９０】
ここでシステムノイズφｔ−１の確率密度関数は以下である。
【００９１】
【数１４】
2016218674.tif 000016
【００９２】
本実験では、状態空間モデルから初期値ｘ０＝０で１０００個の観測値系列を生成し、１００回のモンテカルロシミュレーションを行った。フィルタ推定値＾ｘｔの精度は平均絶対誤差（mean absolute error, ＭＡＥ) で測った（１／１０００Σｔ＝１1000|ｘｔ−＾ｘｔ|）。 比較したフィルタ（ＭＣＦ）と、上記の実施の形態で説明したフィルタ（ＭＣＰＦ, ＧＳＲＦ）は以下である。
【００９３】
１. ＭＣＦ（５０）（非特許文献１）
２. ＭＣＰＦ（５０，１０）
３. ＧＳＲＦ（５０，５）
【００９４】
ＭＣＦ（５０）は５０個の粒子を使用したＭＣＦで、ＭＣＰＦ（５０，１０）は５０個の粒子と１０個の準粒子を使用したＭＣＰＦである。ＧＳＲＦ（５０，５）は５０個のガウシアンコンポーネントと５個の準ガウシアンコンポーネントを使用したＧＳＲＦである。以下の表１にそれぞれフィルタの精度を示す。表１は各フィルタの平均絶対誤差（ＭＡＥ）の１００回のモンテカルロシミュレーションにおける平均値である。
【００９５】
【表１】
2016218674.tif 000017
表１の結果から、ＭＣＰＦ（５０，１０）とＧＳＲＦ（５０，５）の優位性を確認することができる。まとめると、システムノイズのガウス性に限らず、観測ノイズがガウス分布に従う場合はＧＳＲＦが最も良い精度を出すことがわかった。観測ノイズも非ガウス分布に従う場合はＧＳＲＦは適用不可であるため、ＭＣＰＦを使用する。
【００９６】
以上説明したように、本発明の実施の形態に係る状態推定装置によれば、フィルタリング部３０が、観測値ｙｔのノイズが正規分布に従う場合に、ガウス和リサンプリングフィルタを適用し、観測値ｙｔのノイズが正規分布に従わない場合にモンテカルロ粒子フィルタを適用することにより、精度よく、推定対象の状態を観測値から推定することができる。
【００９７】
なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【符号の説明】
【００９８】
１０   入力部
２０   演算部
３０   フィルタリング部
５０   出力部
１００ 状態推定装置
【図１】
2016218674.tif 000018
【図２】
2016218674.tif 000019
【図３】
2016218674.tif 000020