(57)【要約】【課題】非線形モデルを対象として、ノイズの統計量が未知であっても、通常のＵＦＫと同程度の精度で、ロバストな解析結果を得る。【解決手段】状態空間モデリング部３１が、システムノイズを用いた状態更新式、及び観測ノイズを用いた観測方程式によって定められた状態空間モデルを定義し、状態推定部３２が、状態空間モデル及び観測値の時系列データに基づいてＵＦＫを行う際に、システムノイズの共分散行列Ｑの対角成分及び観測ノイズの分散値ｒの総和ｔｒ（Ｑ）＋ｒで正規化された正規化共分散行列＾Ｐ'、及び総和ｔｒ（Ｑ）＋ｒを用いて変更された重みに基づいてシグマ点列を発生させ、＾Ｐ'、ｔｒ（Ｑ）＋ｒで正規化された正規化共分散行列￣Ｐ’、正規化共分散行列＾Ｐ’、正規化システムノイズＱ’、及び正規化観測ノイズｒ’を用いて状態更新及び観測更新を行って、状態ベクトルＸ及び共分散行列Ｐを逐次推定する。【選択図】図１
2014041547.tif 000002
【特許請求の範囲】
【請求項１】
システムノイズを用いて状態ベクトルＸを非線形に時間更新するための状態更新式、及び観測ノイズを用いて状態ベクトルと観測値との関係を示す観測方程式によって定められた状態空間モデル、並びに前記観測値の時系列データに基づいて、シグマ点列の発生、状態ベクトルの状態更新、及び観測更新を行うＵＦＫ（Ｕｎｓｃｅｎｔｅｄ  Ｋａｌｍａｎ  Ｆｉｌｔｅｒ）によって、前記状態ベクトルＸ及び該状態ベクトルの推定誤差分散の共分散行列Ｐを逐次推定する際に、前記状態更新後の前記共分散行列Ｐの予測値￣Ｐをシステムノイズの共分散行列Ｑの対角成分及び観測ノイズの分散値ｒの総和で正規化した正規化共分散行列￣Ｐ’、前記観測更新後の前記共分散行列Ｐの推定値＾Ｐを前記総和で正規化した正規化共分散行列＾Ｐ’、前記システムノイズの共分散行列Ｑを前記総和で正規化した正規化システムノイズＱ’、及び前記観測ノイズの分散値ｒを前記総和で正規化した正規化観測ノイズｒ’を用いると共に、前記総和を用いて変更された重みに基づいて、前記シグマ点列を発生させる状態推定手段と、
前記状態推定手段の推定結果に基づいて、前記正規化システムノイズＱ’と前記正規化観測ノイズｒ’との比率をパラメトライズして得られる尤度が最大となるときのパラメータを、前記状態空間モデルのパラメータとして同定するモデル同定手段と、
を含む時系列データ解析装置。
【請求項２】
前記状態推定手段は、前記総和が収束するまで前記ＵＦＫを繰り返し行い、前記総和が収束した際の前記状態ベクトルＸの推定値＾Ｘ及び前記推定値＾Ｐの各々を平滑化した平滑化推定値~Ｘ及び~Ｐを、最終的な推定結果とする請求項１記載の時系列データ解析装置。
【請求項３】
状態推定手段と、モデル同定手段とを含む時系列データ解析装置における時系列データ解析方法であって、
前記状態推定手段が、システムノイズを用いて状態ベクトルＸを非線形に時間更新するための状態更新式、及び観測ノイズを用いて状態ベクトルと観測値との関係を示す観測方程式によって定められた状態空間モデル、並びに前記観測値の時系列データに基づいて、シグマ点列の発生、状態ベクトルの状態更新、及び観測更新を行うＵＦＫ（Ｕｎｓｃｅｎｔｅｄ  Ｋａｌｍａｎ  Ｆｉｌｔｅｒ）によって、前記状態ベクトルＸ及び該状態ベクトルの推定誤差分散の共分散行列Ｐを逐次推定する際に、前記状態更新後の前記共分散行列Ｐの予測値￣Ｐをシステムノイズの共分散行列Ｑの対角成分及び観測ノイズの分散値ｒの総和で正規化した正規化共分散行列￣Ｐ’、前記観測更新後の前記共分散行列Ｐの推定値＾Ｐを前記総和で正規化した正規化共分散行列＾Ｐ’、前記システムノイズの共分散行列Ｑを前記総和で正規化した正規化システムノイズＱ’、及び前記観測ノイズの分散値ｒを前記総和で正規化した正規化観測ノイズｒ’を用いると共に、前記総和を用いて変更された重みに基づいて、前記シグマ点列を発生させ、
前記モデル同定手段が、前記状態推定手段の推定結果に基づいて、前記正規化システムノイズＱ’と前記正規化観測ノイズｒ’との比率をパラメトライズして得られる尤度が最大となるときのパラメータを、前記状態空間モデルのパラメータとして同定する
時系列データ解析方法。
【請求項４】
コンピュータを、請求項１または請求項２記載の時系列データ解析装置を構成する各手段として機能させるための時系列データ解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、時系列データ解析装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
時系列データの解析において、データの生成プロセスをシステムモデルと観測モデルとで記述する状態空間モデリングがよく用いられる。例えば、線形モデルの残差（観測値と予測値との差）とその共分散との関係を線形式で書き下し、最小二乗法でシステムノイズと観測ノイズとを推定する手法が提案されている（例えば、非特許文献１参照）。
【０００３】
しかし、実際の問題ではモデルが非線形になることが多い。非線形モデルに対する有効な状態推定アルゴリズムにＵｎｓｃｅｎｔｅｄ  Ｋａｌｍａｎ  Ｆｉｌｔｅｒ（ＵＫＦ）が存在する。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】BrianJ.Odelsonetal., “A New Autocovariance Least-Squares Method for Estimating Noise Covariances”, Technical report, TWMCC, 2003.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、従来手法であるＵＦＫでは、システムノイズ及び観測ノイズの分散値の未知性及び初期値任意性により、解析結果がこれらノイズの設定値に大きく依存してしまい、ロバストな解析結果を得られない場合がある、という問題があった。
【０００６】
本発明は、上記の事情を鑑みてなされたもので、非線形モデルを対象として、ノイズの統計量が未知であっても、通常のＵＦＫと同程度の精度で、ロバストな解析結果を得ることができる時系列データ解析装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００７】
上記目的を達成するために、本発明の時系列データ解析装置は、システムノイズを用いて状態ベクトルＸを非線形に時間更新するための状態更新式、及び観測ノイズを用いて状態ベクトルと観測値との関係を示す観測方程式によって定められた状態空間モデル、並びに前記観測値の時系列データに基づいて、シグマ点列の発生、状態ベクトルの状態更新、及び観測更新を行うＵＦＫ（ＵｎｓｃｅｎｔｅｄＫａｌｍａｎＦｉｌｔｅｒ）によって、前記状態ベクトルＸ及び該状態ベクトルの推定誤差分散の共分散行列Ｐを逐次推定する際に、前記状態更新後の前記共分散行列Ｐの予測値￣Ｐをシステムノイズの共分散行列Ｑの対角成分及び観測ノイズの分散値ｒの総和で正規化した正規化共分散行列￣Ｐ’、前記観測更新後の前記共分散行列Ｐの推定値＾Ｐを前記総和で正規化した正規化共分散行列＾Ｐ’、前記システムノイズの共分散行列Ｑを前記総和で正規化した正規化システムノイズＱ’、及び前記観測ノイズの分散値ｒを前記総和で正規化した正規化観測ノイズｒ’を用いると共に、前記総和を用いて変更された重みに基づいて、前記シグマ点列を発生させる状態推定手段と、前記状態推定手段の推定結果に基づいて、前記正規化システムノイズＱ’と前記正規化観測ノイズｒ’との比率をパラメトライズして得られる尤度が最大となるときのパラメータを、前記状態空間モデルのパラメータとして同定するモデル同定手段と、を含んで構成されている。
【０００８】
本発明の時系列データ解析装置によれば、状態空間モデルとして、システムノイズを用いて状態ベクトルＸを非線形に時間更新するための状態更新式、及び観測ノイズを用いて状態ベクトルと観測値との関係を示す観測方程式によって定められたモデルを定義する。そして、状態推定手段が、この状態空間モデルと観測値の時系列データとに基づいて、シグマ点列の発生、状態ベクトルの状態更新、及び観測更新を行うＵＦＫ（ＵｎｓｃｅｎｔｅｄＫａｌｍａｎＦｉｌｔｅｒ、無香カルマンフィルタ）によって、状態ベクトルＸ及び状態ベクトルの推定誤差分散の共分散行列Ｐを逐次推定する。この際、状態推定手段は、状態更新後の共分散行列Ｐの予測値￣Ｐをシステムノイズの共分散行列Ｑの対角成分及び観測ノイズの分散値ｒの総和で正規化した正規化共分散行列￣Ｐ’、観測更新後の共分散行列Ｐの推定値＾Ｐを総和で正規化した正規化共分散行列＾Ｐ’、システムノイズの共分散行列Ｑを総和で正規化した正規化システムノイズＱ’、及び観測ノイズの分散値ｒを総和で正規化した正規化観測ノイズｒ’を用いる。これにより、ＵＦＫアルゴリズム内の全ての共分散行列Ｐの推定値＾Ｐが正規化される。さらに、状態推定手段は、システムノイズの共分散行列Ｑの対角成分及び観測ノイズの分散値ｒの総和を用いて変更された重みに基づいて、シグマ点列を発生させる。これにより、推定精度をテイラー展開の意味で２次まで保証し、ＵＫＦと同等の精度を得ることができる。
【０００９】
また、モデル同定手段が、状態推定手段の推定結果に基づいて、正規化システムノイズＱ’と正規化観測ノイズｒ’との比率を、例えば０〜１で振るなどしてパラメトライズして得られる尤度が最大となるときのパラメータを、状態空間モデルのパラメータとして同定する。
【００１０】
このように、システムノイズ及び観測ノイズを、システムノイズの共分散行列Ｑの対角成分及び観測ノイズの分散値ｒの総和で正規化し、さらに、総和を用いて変更された重みに基づいてシグマ点列を発生させることにより、ノイズの統計量が未知であっても、通常のＵＦＫと同程度の精度で、ロバストな解析結果を得ることができる。
【００１１】
また、前記状態推定手段は、前記総和が収束するまで前記ＵＦＫを繰り返し行い、前記総和が収束した際の前記状態ベクトルＸの推定値＾Ｘ及び前記推定値＾Ｐの各々を平滑化した平滑化推定値~Ｘ及び~Ｐを、最終的な推定結果とすることができる。
【００１２】
また、本発明の時系列データ解析方法は、状態推定手段と、モデル同定手段とを含む時系列データ解析装置における時系列データ解析方法であって、前記状態推定手段は、システムノイズを用いて状態ベクトルＸを非線形に時間更新するための状態更新式、及び観測ノイズを用いて状態ベクトルと観測値との関係を示す観測方程式によって定められた状態空間モデル、並びに前記観測値の時系列データに基づいて、シグマ点列の発生、状態ベクトルの状態更新、及び観測更新を行うＵＦＫ（ＵｎｓｃｅｎｔｅｄＫａｌｍａｎＦｉｌｔｅｒ）によって、前記状態ベクトルＸ及び該状態ベクトルの推定誤差分散の共分散行列Ｐを逐次推定する際に、前記状態更新後の前記共分散行列Ｐの予測値￣Ｐをシステムノイズの共分散行列Ｑの対角成分及び観測ノイズの分散値ｒの総和で正規化した正規化共分散行列￣Ｐ’、前記観測更新後の前記共分散行列Ｐの推定値＾Ｐを前記総和で正規化した正規化共分散行列＾Ｐ’、前記システムノイズの共分散行列Ｑを前記総和で正規化した正規化システムノイズＱ’、及び前記観測ノイズの分散値ｒを前記総和で正規化した正規化観測ノイズｒ’を用いると共に、前記総和を用いて変更された重みに基づいて、前記シグマ点列を発生させ、前記モデル同定手段は、前記状態推定手段の推定結果に基づいて、前記正規化システムノイズＱ’と前記正規化観測ノイズｒ’との比率をパラメトライズして得られる尤度が最大となるときのパラメータを、前記状態空間モデルのパラメータとして同定する方法である。
【００１３】
また、本発明の時系列データ解析プログラムは、コンピュータを、上記の時系列データ解析装置を構成する各手段として機能させるためのプログラムである。
【発明の効果】
【００１４】
以上説明したように、本発明の時系列データ解析装置、方法、及びプログラムによれば、システムノイズ及び観測ノイズを、システムノイズの共分散行列Ｑの対角成分及び観測ノイズの分散値ｒの総和で正規化し、さらに、総和を用いて変更された重みに基づいてシグマ点列を発生させることにより、ノイズの統計量が未知であっても、通常のＵＦＫと同程度の精度で、ロバストな解析結果を得ることができる、という効果が得られる。
【図面の簡単な説明】
【００１５】
【図１】本実施の形態に係る時系列データ解析装置の構成を示す概略図である。
【図２】観測値の時系列データのデータ構造の一例を示す図である。
【図３】本実施の形態における時系列データ解析処理ルーチンの内容を示すフローチャートである。
【図４】実験１における人工時系列データを示すグラフである。
【図５】実験１における１期先の予測結果を示すグラフである。
【図６】実験１におけるモデルの同定結果を示す表である。
【図７】実験２における１期先の予測結果を示すグラフである。
【図８】実験２におけるＡＲ係数の推定結果を示すグラフである。
【図９】実験２におけるモデルの同定結果を示す表である。
【発明を実施するための形態】
【００１６】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００１７】
＜本実施の形態の概要＞
本実施の形態に係る時系列データ解析装置では、まず、通常のＵＫＦアルゴリズムに存在する状態の共分散行列の状態更新後の予測値、及び観測更新後の推定値の各々を、システムノイズの共分散行列の対角成分及び観測ノイズの分散値の総和で割る。すなわち、システムノイズの共分散行列の対角成分及び観測ノイズの分散値の総和を１に正規化する。この操作によりアルゴリズムに存在する全ての共分散行列の推定値は正規化されたものになる。次に、正規化された共分散行列推定値から発生されたシグマ点列を使い時間更新及び観測更新を行って状態の期待値及び共分散行列の推定を行うが、ここで各シグマ点列の重みを適切に変更することで、推定精度をテイラー展開の意味で２次まで保証する。これはＵＫＦと同等の精度である。本実施の形態では説明を簡潔にするため観測モデルは線形を仮定するが、非線形観測モデルに対しても同様なスケーリング処理を行うことでＵＫＦと同等な精度の状態推定を行うことができる。最終的なモデルの同定には対数尤度を使用する。
【００１８】
上記概要に従った本実施の形態における手法を、以下では、Ｎｏｉｓｅｓ  Ｎｏｒｍａｌｉｚｅｄ  ＵＫＦ＆ＵＲＴＳ  Ｓｍｏｏｔｈｅｒ（ＮＮ−ＵＫＦ＆ＵＲＴＳＳ）と呼ぶ。
【００１９】
＜システム構成＞
本実施の形態に係る時系列データ解析装置１０は、ＣＰＵと、ＲＡＭと、後述する時系列データ解析処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成されている。このコンピュータは、機能的には、図１に示すように、入力部２０と、演算部３０と、出力部４０とを含んだ構成で表すことができる。
【００２０】
入力部２０は、入力された観測値の時系列データを受け付ける。観測値は、スカラー値であり、すなわち、時系列データは、１次元時系列である。観測値は、例えば、図２に示すように、時刻ｔと観測値ｙｔとの組（ｔ＝１，・・・，Ｎ）からなる時系列データである。
【００２１】
また、入力部２０は、オペレータから、後述する状態空間モデルを定義するための設定値の入力を受け付ける。
【００２２】
演算部３０は、状態空間モデリング部３１と、状態推定部３２と、モデル同定部３３とを含んだ構成で表すことができる。
【００２３】
状態空間モデリング部３１は、オペレータにより入力部２０から入力された設定値を受け付けて、下記（１）式及び（２）式に示す状態空間モデルを定義する。
【００２４】
【数１】
2014041547.tif 000003
【００２５】
ただし、Ｘｔは時刻ｔにおける状態ベクトルで、通常、非観測量が取られる。状態ベクトルＸｔを時系列データｙｔを用いて推定することを目的とする。μはシステムノイズベクトル、ｗは観測ノイズベクトルであり、共に正規分布に従う白色ノイズベクトルである。また、Ｎ（ｍ、ｖ）は平均ｍ、分散ｖの正規分布を表す。したがって、Ｑはシステムノイズの共分散行列に相当し、ｒは観測ノイズの分散値に相当する。上記（２）式に示されるように、ｙｔは時刻ｔにおける観測値で、状態Ｘｔの線形変換Ｈにより生成される。上記（２）式は、観測ノイズを用いて状態ベクトルＸｔと観測値ｙｔとの関係を示している。ｆ（・）が非線形箇所であり、これによりシステムモデルは非線形になる。ｆ（・）の型はモデルにより異なるが、簡単な循環変動モデルでは、下記（３）式に示すように状態が定義され、状態更新式が組まれる。
【００２６】
【数２】
2014041547.tif 000004
【００２７】
ここで、ｘｔ，・・・，ｘｔ−Ｍ＋１は自己回帰過程で表現された循環変動成分、ａ１，・・・，ａＭは各循環変動成分にかかる係数である。以下、この係数をＡＲ係数と呼ぶ。この状態の定義を使用したシステムモデルにおける状態更新式のｆ（Ｘｔ−１）及びｗは、下記（４）式で表わされる。
【００２８】
【数３】
2014041547.tif 000005
【００２９】
ここでｗは観測ノイズベクトルであるが、通常ｗ１以外の要素は非常に小さな値を設定する（≒０）。この循環変動モデルでは、観測方程式（２）式のＨは１行２×Ｍ列のベクトルで、１番目の要素が１でその他は０になる。
【００３０】
このように、状態空間モデリング部３１は、オペレータからの入力を受け付けて、上記（４）式に示すようなｆ（・）の型を定義すると共に、上記（２）式のＨを定義する。
【００３１】
状態推定部３２は、入力された時系列データと、状態空間モデリング部３１でモデリングされた状態空間モデルとに基づき、状態ベクトルの逐次推定を行う。モデルは非線形であるので、上記概要で説明したＮＮ−ＵＫＦ＆ＵＲＴＳＳによって、状態ベクトルＸｔの逐次状態推定（オンライン状態推定）を行う。ＮＮ−ＵＫＦのアルゴリズムは、下記（５）式〜（１７）式で表されるシグマ点列の発生、状態更新（１期先予測）、及び観測更新（フィルタ）からなる。
【００３２】
【数４】
2014041547.tif 000006
【００３３】
ただし、κは予め定められた定数であり、ｎは状態ベクトルＸｔの次元数である。ｖは予測誤差であり、Ｐｔは状態ベクトルＸｔの推定誤差分散の共分散行列である。また、￣Ａ（数式中では、記号Ａ上に“￣”）は状態更新により予測されたＡの予測値を示し、＾Ａ（数式中では、記号Ａ上に“＾”）は観測更新により推定されたＡの推定値を示す。なお、￣Ｐ'ｔ、＾Ｐ'ｔ、Ｑ’、ｒ’は下記（１８）式〜（２１）式に示すように、共分散行列の予測値￣Ｐｔ、共分散行列の推定値＾Ｐｔ、システムノイズの共分散行列Ｑ、及び観測ノイズの分散値ｒを、ｔｒ（Ｑ）＋ｒで割って正規化したものである。
【００３４】
【数５】
2014041547.tif 000007
【００３５】
ここで、ｔｒ（Ｑ）＋ｒはシステムノイズμの共分散行列Ｑの対角成分及び観測ノイズｗの分散値ｒの総和であり、アルゴリズム中のαとの関係はｔｒ（Ｑ）＋ｒ＝１／α２である。このスケーリング操作により、アルゴリズムに存在するシステムノイズの共分散行列の対角成分及び観測ノイズの分散値は正規化され、全て１以下になる。また、（６）式及び（９）式は、通常のＵＫＦのシグマ点列の重みであるが、これを（７）式及び（１０）式のように変更することで、状態の推定精度をテイラー展開の意味で２次まで保証している。多くの場合、ｔｒ（Ｑ）＋ｒの値は未知であるので最初はｔｒ（Ｑ）＋ｒ＝１、すなわちα＝１に設定してＮＮ−ＵＫＦを実行する。そして、ＮＮ−ＵＫＦの実行毎に、下記（２２）式及び（２３）式に基づいて、ｔｒ（Ｑ）＋ｒを推定する。
【００３６】
【数６】
2014041547.tif 000008
【００３７】
状態推定部３２は、ＮＮ−ＵＫＦ実行時に使用した（ｔｒ（Ｑ）＋ｒ）ｎｏｗの値と、ＮＮ−ＵＫＦ実行終了時に推定した（ｔｒ（Ｑ）＋ｒ）ｎｅｗの値との差が所定の閾値以下であれば、ｔｒ（Ｑ）＋ｒの値が収束したとみなし、この実行によって推定された状態の結果を採用する手法を取る。最後の観測値ｙｔ＝Ｎまでの処理が完了すると、（１３）式を使用したＮｏｉｓｅｓ  Ｎｏｒｍａｌｉｚｅｄ  Ｕｎｓｃｅｎｔｅｄ  ＲＴＳ  Ｓｍｏｏｔｈｅｒ（ＮＮ−ＵＲＴＳＳ）により、下記（２４）式〜（２６）式に示す平滑化を行う。
【００３８】
【数７】
2014041547.tif 000009
【００３９】
ここで、Ｓｔ−１は平滑化ゲイン、~Ｘｔ−１及び~Ｐｔ−１（数式中では、記号ＸまたはＰ上に“~”）はそれぞれ状態の期待値及び共分散行列の平滑化推定値であり、この推定値が最終的にアウトプットされる。
【００４０】
モデル同定部３３は、状態推定部３２により正規化されたシステムノイズの共分散行列Ｑ’及び観測ノイズの分散値ｒ’を利用して、Ｑ’とｒ’との比率γをパラメトライズして振って、γ毎に下記（２７）式〜（３０）式に示す対数尤度ＬＬ（γ）を求める。そして、対数尤度ＬＬ（γ）が最大になったγから得られるＱ’及びｒ’を、状態推定モデルのパラメータとして推定することにより、モデルの同定を行う。
【００４１】
【数８】
2014041547.tif 000010
【００４２】
ここで、Ｍ（γ）がシステムノイズμの共分散行列の対角成分及び観測ノイズｗの分散値の総和の推定値となる。
【００４３】
出力部４０は、状態推定部３２により推定された状態ベクトルの平滑化推定値~Ｘｔ−１及びその共分散行列の平滑化推定値~Ｐｔ−１、並びにモデル同定部３３により同定されたモデルを示すパラメータの推定値を出力する。
【００４４】
＜時系列データ解析装置の作用＞
次に、本実施の形態に係る時系列データ解析装置１０の作用について説明する。まず、オペレータにより、状態空間モデルを定義するための設定値が時系列データ解析装置１０に入力されると、時系列データ解析装置１０によって、入力された設定値が、メモリ（図示省略）へ格納される。また、時刻ｔ１〜ｔＮの観測値からなる時系列データが、時系列データ解析装置１０に入力されると、時系列データ解析装置１０によって、入力された時系列データが、メモリへ格納される。そして、時系列データ解析装置１０によって、図３に示す時系列データ解析処理ルーチンが実行される。
【００４５】
まず、ステップ１００で、メモリに格納された状態空間モデルを定義するための設定値及び観測値の時系列データを取得する。そして、ステップ１０２で、状態空間モデリング部３１が、取得した設定値を用いて、上記（１）式及び（２）式に示す状態空間モデルを定義する。
【００４６】
次に、ステップ１０４で、状態推定部３２が、ＮＮ−ＵＫＦ＆ＵＲＴＳＳを実行するための初期設定として、システムノイズの共分散行列の対角成分及び観測ノイズの分散値の総和（ｔｒ（Ｑ）＋ｒ）、並びに時系列データの時刻を示すｔをそれぞれ１にセットする。
【００４７】
次に、ステップ１０６で、状態推定部３２が、ＮＮ−ＵＫＦのアルゴリズムに従って、上記ステップ１００で取得した観測値の時系列データ、及び上記ステップ１０４で設定された初期設定を用いて、シグマ点列の発生、状態更新、及び観測更新を逐次行って、状態ベクトル及びその共分散行列を推定する。
【００４８】
次に、ステップ１０８で、状態推定部３２が、上記（２２）式及び（２３）式によりｔｒ（Ｑ）＋ｒを推定し、上記ステップ１０６で使用した（ｔｒ（Ｑ）＋ｒ）ｎｏｗの値と、本ステップで推定した（ｔｒ（Ｑ）＋ｒ）ｎｅｗの値との差が所定の閾値以下か否かを判定することにより、ｔｒ（Ｑ）＋ｒの値が収束したか否かを判定する。ｔｒ（Ｑ）＋ｒの値が収束していない場合には、ステップ１０６へ戻り、推定された状態ベクトル＾Ｘｔ、推定され正規化された共分散行列＾Ｐ'ｔ、及び推定されたｔｒ（Ｑ）＋ｒを用いて、ＮＮ−ＵＫＦによる状態ベクトル及びその共分散行列の推定を繰り返す。
【００４９】
一方、ｔｒ（Ｑ）＋ｒの値が収束した場合には、ステップ１１０へ移行し、ｔがＮになったか否かを判定することにより、最後の観測値ｙＮまで処理が完了したか否かを判定する。ｔがＮに達していない場合には、ステップ１１２へ移行して、ｔを１インクリメントして、ステップ１０６へ戻る。ｔ＝Ｎとなった場合には、ステップ１１４へ移行する。
【００５０】
ステップ１１４では、状態推定部３２が、上記（１３）式を使用したＮＮ−ＵＲＴＳＳにより、上記（２４）式〜（２６）式に示す平滑化を行い、状態ベクトルの平滑化推定値~Ｘｔ−１及びその共分散行列の平滑化推定値~Ｐｔ−１を得る。
【００５１】
次に、ステップ１１６で、モデル同定部３３が、Ｑ’とｒ’との比率γをパラメトライズして振って、γ毎に上記（２７）式〜（３０）式に示す対数尤度ＬＬ（γ）を求め、対数尤度ＬＬ（γ）が最大になったγから得られるＱ’及びｒ’を、状態推定モデルのパラメータとして推定することにより、モデルの同定を行う。
【００５２】
次に、ステップ１１８で、出力部４０が、上記ステップ１１４で得られた状態ベクトルの平滑化推定値~Ｘｔ−１及びその共分散行列の平滑化推定値~Ｐｔ−１、並びに上記ステップ１０６で同定されたモデルを示すパラメータの推定値を出力し、時系列データ解析処理ルーチンを終了する。
【００５３】
＜実験結果＞
ここで、下記（３１）式（３次のＡＲ過程、ＡＲ係数は1.5，-0.9，0.3、観測ノイズの分散値は１）で生成された人工時系列データを使用したＮＮ−ＵＫＦ＆ＵＲＴＳＳの実験結果について説明する。
【００５４】
【数９】
2014041547.tif 000011
【００５５】
まず、実験１ではＡＲ係数を既知として循環変動成分のみを推定した。図４に上記（３１）式に基づいて生成された人工時系列データ（ＡＲ（３））、図５に実験１による１期先の予測結果、図６にシステムノイズと観測ノイズの値をλ（＝システムノイズ／（システムノイズ＋観測ノイズ））とした各モデルの同定結果を示す。図４及び図５から１期先の予測が精度良く行われていることがわかり、図６から観測ノイズの分散値（図中、ノイズ推定値）の真値１を推定することに成功していることがわかる。
【００５６】
次に、実験２ではＡＲ係数も未知とし、非線形システムモデル（（１）式〜（４）式）を使用して状態推定を行った。図７に実験２による１期先の予測結果、図８に３つのＡＲ係数の逐次推定（オンライン推定）結果を示す。図８からＡＲ係数のそれぞれの真値（1.5，-0.9，0.3）を適切に推定していることがわかる。また、図９はモデルの同定結果であり、本実施の形態の手法であるＮＮ−ＵＫＦ＆ＵＲＴＳＳの有効性を示している。
【００５７】
以上説明したように、本実施の形態に係る時系列データ解析装置によれば、システムノイズの共分散行列の対角成分及び観測ノイズの分散値の総和を１で正規化して、ＵＫＦアルゴリズムに存在する状態ベクトルの誤差分散の共分散行列の状態更新後の予測値、及び観測更新後の推定値の各々を正規化し、正規化された共分散行列の推定値から発生されたシグマ点列を使って時間更新及び観測更新を行うことにより、ノイズの統計量が未知であっても、ロバストな推定結果を提供することが可能になる。また、システムノイズの共分散行列の対角成分及び観測ノイズの分散値の総和に基づいて適切に変更された重みに基づいてシグマ点列を発生させることにより、推定精度をテイラー展開の意味で２次まで保証することができ、従来のＵＫＦと同等の精度を得ることができる。
【００５８】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【００５９】
例えば、状態空間モデルの定義を予め設定しておき、状態空間モデルを定義するための設定値の入力を不要としてもよい。
【００６０】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【００６１】
１０   時系列データ解析装置
２０   入力部
３０   演算部
３１   状態空間モデリング部
３２   状態推定部
３３   モデル同定部
４０   出力部
【図１】
2014041547.tif 000012
【図２】
2014041547.tif 000013
【図３】
2014041547.tif 000014
【図４】
2014041547.tif 000015
【図５】
2014041547.tif 000016
【図６】
2014041547.tif 000017
【図７】
2014041547.tif 000018
【図８】
2014041547.tif 000019
【図９】
2014041547.tif 000020