(57)【要約】【課題】単位時間に対して取得されていないデータを有する時系列データであっても、時系列モデルパラメータの推定ができる時系列モデルパラメータの推定方法を提供する。【解決手段】時間順で取得された複数系列の時系列データにおいて、取得されたデータ毎の当該データと一つ前に取得されたデータとの差分からなる差分時系列データ、及び、これらの時系列データに共通する最小単位時間に基づいて、時系列モデルパラメータを推定する時系列モデルパラメータの推定方法であって、系列数に対応する多変量正規分布の平均及び分散共分散行列、並びに、差分毎に、当該差分を求めたときのデータと一つ前に取得されたデータとの最小単位時間に基づく時間間隔、及び、複数系列のうち、異なる２つの系列毎の差分の重複する最小単位時間に基づく時間間隔により、差分時系列データの尤度関数を決定し、この尤度関数の値を最大化する平均及び分散共分散行列を求める。【選択図】図４
2017151497.tif 000002
【特許請求の範囲】
【請求項１】
時間順で取得された複数系列の時系列データにおいて、取得されたデータ毎に、当該データと一つ前に取得されたデータとの差分からなる差分時系列データ、及び、前記複数系列の時系列データに共通する最小単位時間に基づいて、前記時系列データの時系列モデルパラメータを推定する時系列モデルパラメータの推定方法であって、
前記複数系列の系列数に対応する多変量正規分布の平均及び分散共分散行列、並びに、前記時系列データの差分毎に、当該差分を求めたときの前記データと前記一つ前に取得されたデータとの前記最小単位時間に基づく時間間隔、及び、前記複数系列のうち、異なる２つの系列毎の前記差分の重複する前記最小単位時間に基づく時間間隔により、前記差分時系列データの尤度関数を決定し、
前記尤度関数の値を最大化する前記差分時系列データの平均及び分散共分散行列を求めることを特徴とする時系列モデルパラメータの推定方法。
【請求項２】
時間順で取得された複数系列の時系列データにおいて、取得されたデータ毎に、当該データと一つ前に取得されたデータとの差分からなる差分時系列データ、及び、前記複数系列の時系列データに共通する最小単位時間に基づいて、前記時系列データの時系列モデルパラメータを推定する時系列モデルパラメータの推定方法であって、
前記複数系列の系列数に対応する多変量正規分布の平均及び分散共分散行列、並びに、前記時系列データの差分毎に、当該差分を求めたときの前記データと前記一つ前に取得されたデータとの前記最小単位時間に基づく時間間隔、及び、前記複数系列のうち、異なる２つの系列毎の前記差分の重複する前記最小単位時間に基づく時間間隔により、前記差分時系列データの尤度関数を決定するステップと、
前記複数系列の時系列データにおいて、前記最小単位時間で前記データが取得されたと仮定したときと、取得された前記時系列データとを比較して、前記データが取得されていない部分にデータを補填して、前記尤度関数の値を最大化する前記差分時系列データの平均及び分散共分散行列を求めるステップと、を有することを特徴とする時系列モデルパラメータの推定方法。
【請求項３】
時間順で取得された複数系列の時系列データにおいて、取得されたデータ毎に、当該データと一つ前に取得されたデータとの差分からなる差分時系列データ、及び、前記複数系列の時系列データに共通する最小単位時間に基づいて、前記時系列データの時系列モデルパラメータを推定する時系列モデルパラメータの推定方法であって、
前記複数系列のうち、一部の系列の差分時系列データを取り出し、前記一部の系列の系列数に対応する多変量正規分布の平均及び分散共分散行列、並びに、前記差分時系列データの前記差分毎に、当該差分を求めたときの前記データと前記一つ前に取得されたデータとの前記最小単位時間に基づく時間間隔、及び、前記複数系列のうち、異なる２つの系列毎の前記差分の重複する前記最小単位時間に基づく時間間隔により、前記一部の系列の尤度関数を決定し、当該尤度関数の値を最大化する前記差分時系列データの平均及び分散共分散行列を求める第１のステップと、
前記複数の系列のうち、残りの未知の平均、分散共分散行列の要素が１つ以上含まれるように、一部の系列の差分時系列データを前記第１のステップと同様に取り出し、前記第１のステップで取得された前記平均及び前記分散共分散行列を既知として尤度関数を決定し、当該尤度関数の値を最大化する前記差分時系列データの平均及び分散共分散行列を求める第２のステップと、を有し、
全ての系列の差分時系列データに対する平均及び分散共分散行列が求まるまで、前記第２のステップを繰り返すことを特徴とする時系列モデルパラメータの推定方法。
【請求項４】
時間順で取得された複数系列の時系列データにおいて、取得されたデータ毎に、当該データと一つ前に取得されたデータとの差分からなる差分時系列データ、及び、前記複数系列の時系列データに共通する最小単位時間に基づいて、前記時系列データの時系列モデルパラメータを推定する時系列モデルパラメータの推定方法であって、
前記複数系列のうち、一部の系列の差分時系列データを取り出し、前記一部の系列の系列数に対応する多変量正規分布の平均及び分散共分散行列、並びに、前記差分時系列データの前記差分毎に、当該差分を求めたときの前記データと前記一つ前に取得されたデータとの前記最小単位時間に基づく時間間隔、及び、前記複数系列のうち、異なる２つの系列毎の前記差分の重複する前記最小単位時間に基づく時間間隔により、前記一部の系列の尤度関数を決定し、前記一部の系列の時系列データにおいて、前記最小単位時間で前記データが取得されたと仮定したときと、取得された前記時系列データとを比較して、前記データが取得されていない部分にデータを補填して、前記尤度関数の値を最大化する前記差分時系列データの平均及び分散共分散行列を求める第１のステップと、
前記複数の系列のうち、残りの未知の平均、分散共分散行列の要素が１つ以上含まれるように、一部の系列の差分時系列データを前記第１のステップと同様に取り出し、前記第１のステップで取得された前記平均及び前記分散共分散行列を既知として尤度関数を決定し、前記残りの系列の一部の系列の時系列データにおいて、前記最小単位時間で前記データが取得されたと仮定したときと、取得された前記時系列データとを比較して、前記データが取得されていない部分にデータを補填して、前記尤度関数の値を最大化する前記差分時系列データの平均及び分散共分散行列を求める第２のステップと、を有し、
全ての系列の差分時系列データに対する平均及び分散共分散行列が求まるまで、前記第２のステップを繰り返すことを特徴とする時系列モデルパラメータの推定方法。
【請求項５】
前記全ての系列の差分時系列データに対する前記差分時系列データの平均及び分散共分散行列が半正定値条件を満たしていないときは、半正定値条件を満たすように、前記平均及び前記分散共分散行列を補正することを特徴する請求項３または４に記載の時系列モデルパラメータの推定方法。
【請求項６】
前記差分時系列データの平均は、前記系列毎の平均と前記差分を求めたときの前記データと前記一つ前に取得されたデータとの前記最小単位時間に基づく時間間隔との積であり、
前記差分時系列データの分散は、前記系列数に対応する多変量正規分布の分散と当該差分を求めたときの前記データと前記一つ前に取得されたデータとの前記最小単位時間に基づく時間間隔との積であり、
前記差分時系列データの共分散のうち、前記異なる２つの系列の同一の時刻の値は、前記系列数に対応する多変量正規分布の共分散と前記異なる２つの系列毎の前記差分の重複する前記最小単位時間に基づく時間間隔との積であり、
前記差分時系列データの共分散のうち、前記異なる２つの系列の異なる時刻の値は０であることを特徴とする請求項１〜５のいずれか一項に記載の時系列モデルパラメータの推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、時系列モデルパラメータの推定方法に関する。
【背景技術】
【０００２】
多次元時系列データ（系列数＝ｎ、時間数＝Ｔ）から時系列モデルパラメータを推定する手法については、欠損が無いケースでは、広く知られている時系列分析手法（Auto Regressive Integrated Moving Average モデル）がある（例えば、特許文献１参照）。この手法は、時系列データの差分（系列数＝ｎ、時間数＝Ｔ−１）が多変量正規分布に従うことにより、標本である差分データから多変量正規分布の母数（平均、分散、共分散）を推定するものである。具体的には、
１．ｎ次元データが多変量正規分布に従うことから尤度を計算
２．時間毎のｎ次元データが独立で同一な分布に従うことから、時間毎に計算された（Ｔ−１）個の尤度の積を計算
３．この積を最大化するような多変量正規分布の母数（平均、分散、共分散）を計算し、推定値とする
の手順を経る。平均、分散、共分散の推定値は、それぞれ標本平均、標本分散、標本共分散に一致し、数式で表せることがわかっている。標本平均がドリフト項、標本分散及び標本共分散が拡散項を表すパラメータの推定値となる。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開平１１−２０３３６１号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
多次元時系列データに欠損がないときは上述した方法によりパラメータの推定値を求めることができるが、実際のデータには、系列毎に様々な欠損が発生する場合や、サンプリング間隔が系列毎やサンプリング毎に異なる場合がある。このような場合、各系列で存在するデータの時間が不揃いとなり、ｎ個のデータの組みの間だけの相関を考慮すればよいという仮定が成立しないため、欠損を伴う（Ｔ−１）個のｎ次元データから多変量正規分布の母数の推定に上述した１〜３の従来技術の手順をそのまま適用することができないという課題があった。なお、以降の説明において、サンプリング間隔が異なる場合であって、最小単位時間で計測した場合に、データが存在しない場合も「欠損」と呼ぶこととする。
【０００５】
本発明はこのような課題に鑑みてなされたものであり、欠損を伴う多次元時系列データであっても、時系列モデルパラメータの推定ができる時系列モデルパラメータ推定方法を提供することを目的とする。
【課題を解決するための手段】
【０００６】
前記課題を解決するために、第１の本発明に係る時系列モデルパラメータの推定方法は、時間順で取得された複数系列の時系列データにおいて、取得されたデータ毎に、当該データと一つ前に取得されたデータとの差分からなる差分時系列データ、及び、前記複数系列の時系列データに共通する最小単位時間に基づいて、前記時系列データの時系列モデルパラメータを推定する時系列モデルパラメータの推定方法であって、前記複数系列の系列数に対応する多変量正規分布の平均及び分散共分散行列、並びに、前記時系列データの差分毎に、当該差分を求めたときの前記データと前記一つ前に取得されたデータとの前記最小単位時間に基づく時間間隔、及び、前記複数系列のうち、異なる２つの系列毎の前記差分の重複する前記最小単位時間に基づく時間間隔により、前記差分時系列データの尤度関数を決定し、前記尤度関数の値を最大化する前記差分時系列データの平均及び分散共分散行列を求めることを特徴とする。
【０００７】
また、第２の本発明に係る時系列モデルパラメータの推定方法は、時間順で取得された複数系列の時系列データにおいて、取得されたデータ毎に、当該データと一つ前に取得されたデータとの差分からなる差分時系列データ、及び、前記複数系列の時系列データに共通する最小単位時間に基づいて、前記時系列データの時系列モデルパラメータを推定する時系列モデルパラメータの推定方法であって、前記複数系列の系列数に対応する多変量正規分布の平均及び分散共分散行列、並びに、前記時系列データの差分毎に、当該差分を求めたときの前記データと前記一つ前に取得されたデータとの前記最小単位時間に基づく時間間隔、及び、前記複数系列のうち、異なる２つの系列毎の前記差分の重複する前記最小単位時間に基づく時間間隔により、前記差分時系列データの尤度関数を決定するステップと、前記複数系列の時系列データにおいて、前記最小単位時間で前記データが取得されたと仮定したときと、取得された前記時系列データとを比較して、前記データが取得されていない欠損部分にデータを補填して、前記尤度関数の値を最大化する前記差分時系列データの平均及び分散共分散行列を求めるステップと、を有することを特徴とする。
【０００８】
また、第３の本発明に係る時系列モデルパラメータの推定方法は、時間順で取得された複数系列の時系列データにおいて、取得されたデータ毎に、当該データと一つ前に取得されたデータとの差分からなる差分時系列データ、及び、前記複数系列の時系列データに共通する最小単位時間に基づいて、前記時系列データの時系列モデルパラメータを推定する時系列モデルパラメータの推定方法であって、前記複数系列のうち、一部の系列の差分時系列データを取り出し、前記一部の系列の系列数に対応する多変量正規分布の平均及び分散共分散行列、並びに、前記差分時系列データの前記差分毎に、当該差分を求めたときの前記データと前記一つ前に取得されたデータとの前記最小単位時間に基づく時間間隔、及び、前記複数系列のうち、異なる２つの系列毎の前記差分の重複する前記最小単位時間に基づく時間間隔により、前記一部の系列の尤度関数を決定し、当該尤度関数の値を最大化する前記差分時系列データの平均及び分散共分散行列を求める第１のステップと、前記複数の系列のうち、残りの未知の平均、分散共分散行列の要素が１つ以上含まれるように、一部の系列の差分時系列データを前記第１のステップと同様に取り出し、前記第１のステップで取得された前記平均及び前記分散共分散行列を既知として尤度関数を決定し、当該尤度関数の値を最大化する前記差分時系列データの平均及び分散共分散行列を求める第２のステップと、を有し、全ての系列の差分時系列データに対する平均及び分散共分散行列が求まるまで、前記第２のステップを繰り返すことを特徴とする。
【０００９】
また、第４の本発明に係る時系列モデルパラメータの推定方法は、時間順で取得された複数系列の時系列データにおいて、取得されたデータ毎に、当該データと一つ前に取得されたデータとの差分からなる差分時系列データ、及び、前記複数系列の時系列データに共通する最小単位時間に基づいて、前記時系列データの時系列モデルパラメータを推定する時系列モデルパラメータの推定方法であって、前記複数系列のうち、一部の系列の差分時系列データを取り出し、前記一部の系列の系列数に対応する多変量正規分布の平均及び分散共分散行列、並びに、前記差分時系列データの前記差分毎に、当該差分を求めたときの前記データと前記一つ前に取得されたデータとの前記最小単位時間に基づく時間間隔、及び、前記複数系列のうち、異なる２つの系列毎の前記差分の重複する前記最小単位時間に基づく時間間隔により、前記一部の系列の尤度関数を決定し、前記一部の系列の時系列データにおいて、前記最小単位時間で前記データが取得されたと仮定したときと、取得された前記時系列データとを比較して、前記データが取得されていない欠損部分にデータを補填して、前記尤度関数の値を最大化する前記差分時系列データの平均及び分散共分散行列を求める第１のステップと、前記複数の系列のうち、残りの未知の平均、分散共分散行列の要素が１つ以上含まれるように、一部の系列の差分時系列データを前記第１のステップと同様に取り出し、前記第１のステップで取得された前記平均及び前記分散共分散行列を既知として尤度関数を決定し、前記残りの系列の一部の系列の時系列データにおいて、前記最小単位時間で前記データが取得されたと仮定したときと、取得された前記時系列データとを比較して、前記データが取得されていない欠損部分にデータを補填して、前記尤度関数の値を最大化する前記差分時系列データの平均及び分散共分散行列を求める第２のステップと、を有し、全ての系列の差分時系列データに対する平均及び分散共分散行列が求まるまで、前記第２のステップを繰り返すことを特徴とする。
【発明の効果】
【００１０】
本発明に係る時系列モデルパラメータの推定方法を以上のように構成すると、単位時間に対して取得されていないデータを有する多次元時系列データであっても、時系列モデルパラメータの推定ができる。
【図面の簡単な説明】
【００１１】
【図１】欠損があるときの時系列データに対する差分時系列データを説明するための説明図である。
【図２】３次元の４時刻分の差分時系列データを説明するための説明図である。
【図３】時系列モデルパラメータ推定システムの構成を示す説明図である。
【図４】時系列モデルパラメータの推定方法のフローチャートである。
【発明を実施するための形態】
【００１２】
以下、本発明の好ましい実施形態について図面を参照して説明する。
【００１３】
１．時系列データ｛ｙi,t｝の定義
まず、系列ｎ、時間数Ｔとするときの時系列データを｛ｙi,t｝（ｉ＝１〜ｎ、ｔ＝１〜Ｔ）とすると、この時系列データ｛ｙi,t｝の隣接するデータ同士の差分データ（差分時系列データ）である｛ｘi,t｝（ｉ＝１〜ｎ、ｔ＝１〜（Ｔ−１））を、ｘi,t＝ｙi,(t+1)−ｙi,tと定義する。ここで、ｔは最小単位時間に基づいて振り出されている。時系列データが、対数を取るなど所定の変換をした結果として「拡散項が正規分布に従うドリフト（定数）項付きランダムウォークモデル」に従う場合には、以下に示す式（１）が成立する。ここで、ｘは差分時系列データであり、Ｎは正規分布であり、μは平均であり、Σは分散共分散行列である。
【００１４】
【数１】
2017151497.tif 000003
【００１５】
このような差分時系列データから、全てのｉ，ｊ，ｔ，ｕに対して、平均Ｅ、分散Ｖ及び共分散Ｃｏｖについては、以下に示す式（２）〜（５）が成立する。なお、ｊはｉとは異なる系列（ｉ≠ｊ）を示している。
【００１６】
【数２】
2017151497.tif 000004
【００１７】
２．欠損データがある場合について
次に、ｉ番目（ｉ＝１〜ｎ）のデータの一部に欠損があるときの差分時系列データについて説明する。図１は、ｙi,t(s)とｙi,t(s+1)の間のデータが欠損している場合を示している。ここで、ｓは非欠損データだけを並べた際に、何番目であるかを示し、ｔ（ｓ）はｓの関数でｓ番目の時刻を表している。例えば、この図１においては、ｙi,t(s)の次にｙi,t(s+1)が取得できたことを示している。このとき、欠損データ間の差分データは取得することができず、ｙi,t(s)とｙi,t(s+1)の差分データ（ｘi,s*とする）だけを取得することができる。すなわち、欠損データがあるときの差分時系列データ（ｘi,s*）は、以下に示す式（６）で表される。
【００１８】
【数３】
2017151497.tif 000005
【００１９】
ここで、式（７）に示すように、差分データの和の期待値Ｅは、各々の差分データの期待値の和で表され、また、差分データの和の分散Ｖは、各々の差分データの和とそれらの差分データの共分散で表される。上述した式（２）〜（５）の性質を用いると、以下に示す式（７）、（８）が成立する。
【００２０】
【数４】
2017151497.tif 000006
【００２１】
この式（７）、（８）の関係より、上述した式（６）で表される差分時系列データの平均Ｅ及び分散Ｖはそれぞれ、式（９）、（１０）で表される。
【００２２】
【数５】
2017151497.tif 000007
【００２３】
また、欠損を伴う２つの時系列データにおいて、差分時系列データの共分散は、以下に示す式（１１）で表される。ここで、２つの時系列データをｉ，ｊとし、非欠損データの順番をｓ，ｕで表すものとする。なお、式（１１）において、時刻ｔのときが重複しており、（ｔ−１），（ｔ＋１）のときは重複していないことを示している。また、時刻の異なるデータの共分散の値は０となる。
【００２４】
【数６】
2017151497.tif 000008
【００２５】
この式（１１）より、欠損を伴う２つの時系列データそれぞれの差分時系列データ（ｘi,s*とｘj,u*）間の共分散は、以下に示す式（１２）で表される。
【００２６】
【数７】
2017151497.tif 000009
【００２７】
したがって、式（９）、（１０）、（１２）より、欠損を伴う２つの時系列データの差分時系列データは、式（１３）に示す正規分布に従うことになる。
【００２８】
【数８】
2017151497.tif 000010
【００２９】
３．欠損がある差分時系列データからのモデルパラメータの最尤推定
上述したように、本実施形態においては、差分時系列データが多変量正規分布に従うことから、当該正規分布の確率密度関数は、当該データを既知とした場合のモデルパラメータの関数（尤度関数）とみることができる。したがって、これを最大化するモデルパラメータ（μi，σi2，σij）の値を求めることで、最尤推定量が得られる。なお、差分データの系列間は、系列数を次元とする多変量正規分布に従うことから、式（１）に示す分散共分散行列Σ（＝｛σi2，σij｝）が半正定値行列でなければならない。つまり、Σの全固有値が非負となる条件下で上記尤度関数を最大化することになる。
【００３０】
上述した欠損がある差分時系列データについて、図２に示すように、３次元の４時刻分を例にして説明する。ここでは、ｘ、ｙ、ｚの３つの系列があり、上段に示すように、それぞれ４時刻分のデータが発生するが、下段に示すように、一部のデータに欠損が発生している。具体的には、系列ｘでは、ｘ1が取得できずにｘ1及びｘ2の合計であるｘ1*が取得され、ｘ3が取得できずにｘ3及びｘ4の合計であるｘ2*が取得されている。また、系列ｙでは、ｙ1はｙ1*として取得されているが、ｙ2が取得できずにｙ2及びｙ3の合計であるｙ2*が取得され、ｙ4がｙ3*として取得されている。また、系列ｚでは、ｚ1はｚ1*として取得されているが、ｚ2〜ｚ4が取得できずにｚ2〜ｚ4の合計であるｚ2*が取得されている。
【００３１】
上述した欠損を含む差分時系列データをＸとし、平均をμとし、分散共分散行列をΣとすると、差分時系列データＸは、平均μ及び分散共分散行列Σの正規分布に従い、以下に示す式（１４）で表される。
【００３２】
【数９】
2017151497.tif 000011
【００３３】
ここで、差分時系列データＸは、取得された全てのデータのベクトルである。そして、平均μは、式（９）を用いて説明したように、各系列の平均（例えば、系列ｘにおけるμx）に各時刻の長さ（個数）をかけた値である。具体的には、図２に示すように、データｘ1*は、ｘ11とｘ12に対応しているため、その時刻の長さ（個数）は「２」となり、その平均は２μxとなる。他のデータについても同様である。
【００３４】
次に、分散共分散行列Σは、式（１３）を用いて説明したように、分散は、平均と同様に、各系列の分散（例えば、系列ｘにおけるσx2）に各時刻の長さ（個数）をかけた値である。具体的には、図２に示すように、ｘ1*の分散は２σx2となる。また、共分散は、重複する時刻がないときは「０」となり、重複する時刻があるときは、共分散の値に重複する時刻の長さ（個数）をかけたものである。具体的には、図２に示すように、ｘ1*とｙ3*とは重複する時刻がないため「０」となり、ｘ1*とｙ2*とは２番目の時刻のデータが重複するため、系列ｘと系列ｙとの共分散σxyにその重複する時刻の長さ（個数）である「１」をけたσxyとなる。
【００３５】
以上より求められた平均μ及び分散共分散行列Σより、尤度関数Ｌを求めると、以下に示す式（１５）で表され、また、対数尤度ＬＬは式（１６）で表され、この対数尤度ＬＬを最大化する平均μ及び分散共分散行列Σを求めることで時系列モデルパラメータの値を決定することができる。なお、式（１５）及び式（１６）の第１項の「７」は、図２に示すように、取得されたデータが、ｘ1*〜ｚ2*の数が７個であることが理由である。
【００３６】
【数１０】
2017151497.tif 000012
【００３７】
なお、系列数３を次元とする多変量正規分布の分散共分散行列は、以下に示す式（１７）で表され、この分散共分散行列の固有値が非負となる条件下での尤度関数Ｌの最大化を行うことが必要である。
【００３８】
【数１１】
2017151497.tif 000013
【００３９】
４．欠損が全くない場合
ここまで説明したように、時系列データに欠損がある場合でも、得られた差分データから時系列モデルパラメータを推定することができる。そこで、この方法が、欠損が全くない場合でも適用可能なことについて説明する。
【００４０】
欠損を考慮した場合の差分時系列データｘi,t(s)において、欠損がないときはｔ（ｓ）＝ｓとなる。したがって、欠損がないときの差分時系列データをｘi,s（ｉ＝１〜Ｎ，ｓ＝１〜Ｓ（＝Ｔ−１））とおくと、差分時系列データＸ、平均μは以下に示す式（１８）、（１９）で表される。なお、ｍは推定された平均を示す。
【００４１】
【数１２】
2017151497.tif 000014
【００４２】
また、推定されるべき分散共分散行列を以下の式（２０）に示すＶとすると、全体の分散共分散行列Σは式（２１）で表される。なお、Ｏは０行列である。
【００４３】
【数１３】
2017151497.tif 000015
【００４４】
以上より、対数尤度ＬＬは以下に示す式（２２）で表される。
【００４５】
【数１４】
2017151497.tif 000016
【００４６】
この式（２２）は、多変量正規分布に従うＮ次元ベクトルがＳ個あって、それらが独立である場合の対数尤度ＬＬにほかならない。この場合の平均ｍ及び推定された分散共分散行列Ｖにかかる最尤推定量については、解析解が知られており、以下の式（２３）、（２４）で示される通り、標本平均、標本分散共分散行列に一致する。また、この標本分散共分散行列が半正定値行列となることも知られている。
【００４７】
【数１５】
2017151497.tif 000017
【００４８】
５．一次元の場合
ここまでの説明では複数の系列の時系列データについて説明してきたが、ここでは１次元の場合について説明する。上述した式（６）より、１次元の差分データｘis*は、以下に示す式（２５）のように正規分布に従う。なお、１次元であるため、ｉ＝１となる。また、差分時系列データＸ、平均μ及び分散共分散行列Σも示す。
【００４９】
【数１６】
2017151497.tif 000018
【００５０】
そして、この式（２５）より、尤度関数Ｌ及び対数尤度ＬＬは式（２６）、（２７）で表される。
【００５１】
【数１７】
2017151497.tif 000019
【００５２】
ここで、式（２７）に示される対数尤度ＬＬを最大化する平均μi及びσi2は、解析的に求めることができ、最尤推定量は、式（２８）、（２９）となる。なお、もし欠損が全くない場合には、以下の式（２８），（２９）は標本平均、標本分散に一致する。
【００５３】
【数１８】
2017151497.tif 000020
【００５４】
６．導出した尤度関数を最大化するモデルパラメータの計算手法（一括使用による解法）
以上のようにして導出された尤度関数に基づいて、この尤度関数を最大化するモデルパラメータ（平均、分散、共分散）を求める方法として、まず、一括使用による解法を説明する。この一括使用による解法では、尤度関数から数理最適化手法（例えば、準ニュートン法等）を用いてモデルパラメータを求める方法である。但し、制約条件（半正定値条件）及び目的関数（対数尤度最大化）がともに非線形となるため、モデルパラメータを直接変えながらモデルパラメータの最尤推定量を探索するのではなく、欠損部分に入れるデータを変化させながら、都度見かけ上欠損を無くし、その上で、式（２３），（２４）に示した解析解に従って計算し、しかる後、欠損を考慮した尤度関数の値を計算することを繰り返しながら探索する。なお、この一括使用による解法では、式（２４）が半正定値行列であることから、半正定値条件を充足しながら、対数尤度を最大化するモデルパラメータを探索することができる。
【００５５】
７．導出した尤度関数を最大化するモデルパラメータの計算手法（逐次使用による解法）
ここでは、ｎ個の系列の中から一部の系列を抜き出して、モデルパラメータを計算し、ついで、残りの未知のパラメータが１つ以上含まれるように一部の系列を前回同様に抜き出し、既に計算されているモデルパラメータを固定値として未知のパラメータを計算し、上記の処理を繰り返して全てのモデルパラメータを計算するという方法である。なお、この方法の場合、全てのモデルパラメータが得られた後、そのモデルパラメータにより半正定値条件を満足するか否かを確認し、半正定値条件を満たしていない場合は、モデルパラメータから定まる分散共分散行列が半正定値条件を満たすように変形し補正解を求めることが必要である。
【００５６】
例えば、まず、ｎ個の系列を個々に１次元系列ととらえ、上述した式（２８），（２９）により、平均と分散の最尤推定量を計算し、ついで、ｎ個の系列から２系列を取り出す全ての組み合わせを考え、それぞれの組み合わせを２次元系列ととらえ、１次元系列として得た平均と分散の推定値を固定した対数尤度から、共分散の最尤推定量を計算し、最後に、半正定値条件を確認し、この条件を満たしていない場合には適宜変形するという方法がある。なお、このような解析的な計算方法の代わりに、上述した数理最適化手法による計算方法を用いることも可能である。
【００５７】
なお、半正定値条件を満たさない場合の変形方法としては、上述の方法により得られた分散共分散行列において、条件を満たさない推定分散共分散行列に対し、条件を満たす分散共分散行列を別に考え、両者の差（例えば、何らかの行列ノルム）が最小となるものを見いだす方法がある。あるいは、上述の方法により得られた分散共分散行列において、条件を満たさない推定分散共分散行列を固有値分解し、対角行列の対角成分の値を非負となるよう調整し、分解式に当てはめて分散共分散行列を修正するという方法がある。
【００５８】
８．時系列モデルパラメータ推定システム
次に、これまで説明した時系列モデルパラメータの推定を行う時系列モデルパラメータ推定システム１００について図３及び図４を用いて説明する。この時系列モデルパラメータ推定システム１００は、図３に示すように、ＣＰＵやＲＡＭ、ＲＯＭ等を有し、プログラムを実行することにより時系列モデルパラメータの推定を行う処理部１１０と、所定の情報を入力する入力部１２０と、推定結果等を記憶するためのメモリーやハードディスク等からなる記憶部１３０と、推定結果等を出力する出力部１４０と、を有して構成されている。
【００５９】
図４を用いて、処理部１１０で実行される処理について説明する。処理部１１０は、時系列モデルパラメータ推定処理が実行されると、まず、差分時系列データ（欠損を含む場合がある）を読み込む（ステップＳ２００）。なお、予め時系列データ（例えば、取得時刻と取得されたデータの組み合わせ）を記憶部１３０に記憶しておいてこれを読み込むように構成してもよいし、時系列データを読み込んで、この時系列データから差分時系列データを生成するように構成してもよい。そして、一括使用か逐次使用かを判断する（ステップＳ２０２）。どちらの方法を選択するかを予め記憶部１３０に設定しておいてもよいし、入力部１２０からの入力により選択するように構成してもよい。
【００６０】
ステップＳ２０２において、一括使用が選択された場合には、上述したように、欠損が無い状態での時系列データから求められる差分時系列データを生成し（ステップＳ２０４）、この欠損のない時系列データから求められる差分時系列データから平均、分散、共分散を計算し（ステップＳ２０６）、その結果から尤度（対数尤度）を計算する（ステップＳ２０８）。このステップＳ２０４〜Ｓ２０８の処理を、欠損データを代えて繰り返し実行し、尤度が最大になったか否かを判断し（ステップＳ２１０）、最大になったと判断したときは、そのときの平均、分散、共分散を時系列モデルパラメータの推定値として出力する（ステップＳ２２４）。出力方法として、ディスプレイやプリンタ等からなる出力部１４０に出力（表示）してもよいし、記憶部１３０に記憶させてもよい。
【００６１】
一方、ステップＳ２０２において、逐次使用が選択された場合には、全系列データから一部の系列データを抜き出し、その系列に対して、上述したように、時系列モデルパラメータ（平均、分散共分散行列）を計算する（ステップＳ２１２）。ついで、残りの未知のパラメータが１つ以上含まれるように一部の系列を前回同様に抜き出し（ステップＳ２１４）、既に計算されている系列の時系列モデルパラメータを既知の値として未知の時系列モデルパラメータを計算し（ステップＳ２１６）、全ての系列の計算が終了するまでステップＳ２１４〜Ｓ２１６を繰り返す（ステップＳ２１８）。そして、全ての系列が計算されたと判断すると、上述した処理により算出された最終的な時系列モデルパラメータが半正定値条件を満たすか否かを判断し（ステップＳ２２０）、満たさないときは補正解を算出して補正し（ステップＳ２２２）、結果を平均、分散、共分散を時系列モデルパラメータの推定値として出力する（ステップＳ２２４）。
【実施例】
【００６２】
それでは、具体的な実施例として、下記の表１に示すように、ｘ，ｙ，ｚの３つの系統における７時刻分の時系列データに対し、上述した一括使用による解法及び逐次使用による解法に基づいて時系列モデルパラメータを求める場合について説明する。なお、以下の表１において、「？」はそのデータが欠損していることを示している。また、各系列における「差」の値は、当該時刻におけるデータと、一つ前に取得されたデータとの差分を示している。
【００６３】
【表１】
2017151497.tif 000021
【００６４】
また、この表１に示す時系列データ（差分時系列データ）に対して、全ての差分を、ベクトルにした差分時系列データＸ、並びに、この差分時系列データＸから求められる平均μ及び分散共分散行列Σは以下のように表される。なお、Ｄは差分時系列データＸのデータ数を示している（ｘtに対して４個の差分データが取得され、ｙtに対して６個の差分データが取得され、ｚtに対して５個の差分データが取得されているため、合計でＤ＝１５となる）。
【００６５】
【数１９】
2017151497.tif 000022
【００６６】
また、以降の説明において、対数尤度ＬＬを、データ数Ｄと差分時系列データＸをパラメータに有する、平均μ及び分散共分散行列Σの関数として以下のように定義し、この対数尤度ＬＬが最大になる時系列モデルパラメータを求めるものとする。
【００６７】
【数２０】
2017151497.tif 000023
【００６８】
但し、この対数尤度ＬＬにおいて、系列xt〜ztの3系列（3次元）の多変量正規分布の分散共分散行列をＶとしたとき、この分散共分散行列Ｖは、次式のように表され、その固有値が非負となる（すなわち、半正定値行列となる）条件下での対数尤度ＬＬの最大化を行う。
【００６９】
【数２１】
2017151497.tif 000024
【００７０】
（第１の実施例）
まず、第１の実施例として、上述した差分時系列データＸに対し、一括解法（一括使用による解法）で、時系列モデルパラメータを推定する場合について説明する。
【００７１】
対数尤度ＬＬの最大化にあたって、上述した分散共分散行列Ｖの半正定値条件を充足させるために、欠損した時系列データに対して何らかの値を与えて欠損無しの状態を生成し、この場合に理論的に得られる平均μ及び分散共分散行列Ｖの推定値を用いて対数尤度ＬＬを求める。ここで、欠損無しの場合の理論値は必ず半正定値条件を満たすことから、欠損を補填する値を変えて計算を行うことで、対数尤度ＬＬを最大とする平均μ及び分散共分散行列Ｖを見いだすことにすれば、半正定値条件の下で対数尤度ＬＬの最大化を行うこととなる。
【００７２】
まず、表１の値より、ｘ、ｙ、ｚ各系列の平均μx、μy、μz、並びに、分散共分散行列Ｖの値を、式（２３）、（２４）を用いて求める。
【００７３】
【数２２】
2017151497.tif 000025
【００７４】
上述したように、欠損データであるｘ３，ｘ4，ｚ4を変化させながら、対数尤度ＬＬが最大となる平均μ及び分散共分散行列Ｖを求める。
【００７５】
（第２の実施例）
次に、第２の実施例として、逐次解法（逐次使用による解法）で、時系列モデルパラメータを推定する場合について説明する。ここでは、ｘt、ｙt、ｚtからなる系列数３の時系列データに対し、それぞれ１次元データとしてモデルパラメータを算出し、次に、３つの系列から２つの系列を選択し、１次元の結果を用いることで、２次元データとしてモデルパラメータを算出し、それらの結果から、最終的なモデルパラメータの推定値を求める場合について説明する。
【００７６】
［１］ｘtだけの一次元データとしたときのモデルパラメータ
以下の表３は、表１から系列ｘtだけの時系列データ及びその差分を抽出したものである。
【００７７】
【表２】
2017151497.tif 000026
【００７８】
この表２から、式（２８）（２９）により系列ｘtの平均μx及び分散σ2xを求める。これをｍx及びｓx2（但しｓx≧０）とする。
【００７９】
［２］ｙtだけの一次元データとしたときのモデルパラメータ
表１から系列ｙtだけの時系列データ及びその差分を抽出し、［１］と同様の手法により平均μy及び分散σy2を求める。これをｍy及びｓy2（但しｓy≧０）とする。
【００８０】
［３］ｚtだけの一次元データとしたときのモデルパラメータ
表１から系列ｚtだけの時系列データ及びその差分を抽出し、［１］と同様の手法により平均μz及び分散σz2を求める。これをｍz及びｓz2（但しｓz≧０）とする。
【００８１】
［４］ｘt，ｙtの二次元データとしたときのモデルパラメータ
以下の表３は、表１から系列ｘt，ｙtの時系列データ及びその差分を抽出したものである。
【００８２】
【表３】
2017151497.tif 000027
【００８３】
系列ｘt，ｙtにおいて、差分データをベクトル化したものをＸxyとし、そのときの平均をμxy、分散共分散行列をΣxyとすると、以下のように表される。ここで、μx、μy、σx2、σy2は、一次元のデータとして求められた値を使用する。また、系列ｘt，ｙtにおいては、５個の差分データが取得されているため、Ｄ＝１０となる。
【００８４】
【数２３】
2017151497.tif 000028
【００８５】
そして、ＬＬ＝ＬＬ（μxy，Σxy；１０，Ｘxy）が最大となるσxyを求める。なお、系列
ｘt及び系列ｙtの２系列（２次元）多変量正規分布の分散共分散行列Ｖxyは、以下のように表される。
【００８６】
【数２４】
2017151497.tif 000029
【００８７】
ここで、分散共分散行列Ｖxyが半正定値行列、すなわち、σxyが次式の範囲であることが必要である。
【００８８】
【数２５】
2017151497.tif 000030
【００８９】
以上より、（μx、μy、σx2、σy2）を（ｍx、ｍy、ｓx2、ｓy2）に固定し、σxyを前式の範囲で変化させながら、対数尤度ＬＬが最大となるσxyを求める。これをsxyとする。
【００９０】
［５］ｘt，ｚtの二次元データとしたときのモデルパラメータ
表１から系列ｘt，ｚtの時系列データ及びその差分を抽出し、［４］と同様の手法により共分散σxzを求める。これをｓxzとする。
【００９１】
［６］ｙt，ｚtの二次元データとしたときのモデルパラメータ
表１から系列ｙt，ｚtの時系列データ及びその差分を抽出し、［４］と同様の手法により共分散σyzを求める。これをｓyzとする。
【００９２】
［７］推定値の補正
上述した（１）〜（６）から、本実施例における時系列モデルパラメータの推定値は以下のようになる。
【００９３】
【数２６】
2017151497.tif 000031
【００９４】
ここで、分散共分散行列Ｖを固有値分解すると以下のようになる。
【００９５】
【数２７】
2017151497.tif 000032
【００９６】
この分散共分散行列Ｖにおいて、右辺の２番目の行列には、対角成分に固有値が並んでいるが、負となる成分が１つ以上ある場合には、この分散共分散行列Ｖは半正定値行列となっていない。このため、半正定値行列となるように、負の成分を補正する。以下の例は、正の微小数ε（＝０．０００１）に補正した場合を示している（補正後の分散共分散行列をＶ*とする）。
【００９７】
【数２８】
2017151497.tif 000033
【００９８】
このが最終的な推定値となる。
【００９９】
以上のように、本実施形態に係る時系列モデルパラメータの推定方法によれば、多次元時系列データの一部に欠損があったとしても、差分データのパターンに応じた同時分布を生成できるようになったため、尤度さらにはこれを最大化する母数の値を計算することが可能になった。これにより、欠損を伴う多次元時系列データであっても、欠損していない残りの全てのデータを使用した時系列モデルパラメータの推定が可能となる。なお、欠損を伴う多次元時系列データだけでなく、最小単位時間は共通するが、この最小単位時間の整数倍（系列毎に倍数が異なる場合や、１つの系列内でも倍数が一定でない場合、及びその組み合わせも含む）でデータが取得される時系列データに対しても、単位時間に対して取得されていないデータを上記欠損と考えれば、本発明を適用することができる。
【符号の説明】
【０１００】
１００  時系列モデルパラメータ推定システム
１１０  処理部
１２０  入力部
１３０  記憶部
１４０  表示部
【図１】
2017151497.tif 000034
【図２】
2017151497.tif 000035
【図３】
2017151497.tif 000036
【図４】
2017151497.tif 000037