(57)【要約】  【課題】  入力を伴う系から得られるデータからなる時系列のパラメータを確実に推定する方法、現実に即したデータの予測を行う方法、パラメータ推定装置、データ予測装置及びコンピュータプログラムを提供する。  【解決手段】  データ予測装置（パラメータ推定装置）は、現在の標本時系列｛ｘn ｝と未来の入力時系列｛ｕn ｝との相互相関関数を０としてデータを計算し（Ｓ２）、白色雑音｛ｅn ｝を生成し（Ｓ３）、｛ｅn ｝の分散σe2を計算し（Ｓ４）、｛ｕn ｝と｛ｅn ｝との相互相関関数、現在の｛ｘn ｝と未来の｛ｅn ｝との相互相関関数を０としてデータを計算し（Ｓ５）、｛ｅn ｝の自己相関関数をラグ０でσe2、他のラグで０として、ＡＲＭＡＸモデルのパラメータを計算する（Ｓ６）。パラメータの収束を判定し（Ｓ７）、収束していない場合は新たな白色雑音を計算し（Ｓ８）、パラメータが収束するまで計算を繰り返す。【選択図】        図２
2005242580.tif 000002
【特許請求の範囲】
【請求項１】
記憶部及び演算部を備えたコンピュータを用いて、複数のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…｝及び入力時系列｛ｕn ：ｎ＝１，２，…｝を（ｐ，ｐ，ｐ）次（但し、ｐは自然数）のＡＲＭＡＸモデルに当てはめたときのパラメータを推定する方法において、
前記標本時系列及び前記入力時系列を前記記憶部で記憶し、
下記式で定義される行列Ｒ，Ｗ，Ｖ，ｒ，ｖを前記記憶部で記憶し、
【数１】
2005242580.tif 000031
平均値が０である複数のランダムなデータからなる白色雑音｛ｅn ：ｎ＝１，２，…｝を前記演算部で生成し、
前記演算部で生成した白色雑音の分散σe2を前記演算部で計算し、
下記式で定義される行列Ｔ及びｔを前記演算部で計算し、
【数２】
2005242580.tif 000032
前記ＡＲＭＡＸモデルにおける前記標本時系列に係るｐ個のパラメータを行列要素とするベクトルａ＝［ａ1 ，ａ2 ，…，ａp ］T 、前記入力時系列に係るｐ個のパラメータを行列要素とする行列ｃ＝［ｃ1，ｃ2 ，…，ｃp ］T 及び前記白色雑音に係るｐ個のパラメータを行列要素とする行列ｂ＝［ｂ1，ｂ2 ，…，ｂp ］T を、下記式に基づいて前記演算部で計算し、
【数３】
2005242580.tif 000033
前記演算部で計算したベクトルａ、ｂ及びｃに含まれる前記ＡＲＭＡＸモデルのパラメータ、並びに前記標本時系列及び前記入力時系列値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を前記演算部で計算し、
前記演算部で計算した新たな白色雑音を用いて新たなベクトルａ、ｂ及びｃを前記演算部で計算し、
ベクトルａ、ｂ及びｃが収束するまでベクトルａ、ｂ及びｃの計算を前記演算部で繰り返すこと
を特徴とするパラメータ推定方法。
【請求項２】
記憶部及び演算部を備えたコンピュータを用いて、入力を伴う系から得られる複数のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…｝、及び前記系へ入力される複数のデータからなる入力時系列｛ｕn：ｎ＝１，２，…｝から、前記標本時系列に加わるべき新たなデータを予測する方法において、
前記標本時系列及び前記入力時系列を前記記憶部で記憶し、
請求項１に記載のパラメータ推定方法を用いて、前記標本時系列及び前記入力時系列をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を前記演算部で計算し、
前記演算部で計算した前記パラメータの推定値、並びに前記標本時系列及び前記入力時系列を前記ＡＲＭＡモデルの式に代入して、前記標本時系列に加わるべき新たなデータの予測値を前記演算部で計算すること
を特徴とするデータ予測方法。
【請求項３】
複数のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…｝及び入力時系列｛ｕn ：ｎ＝１，２，…｝を（ｐ，ｐ，ｐ）次（但し、ｐは自然数）のＡＲＭＡＸモデルに当てはめたときのパラメータを推定する装置において、
前記標本時系列及び前記入力時系列を記憶する手段と、
下記式で定義される行列Ｒ，Ｗ，Ｖ，ｒ，ｖを記憶する手段と、
【数４】
2005242580.tif 000034
平均値が０である複数のランダムなデータからなる白色雑音｛ｅn ：ｎ＝１，２，…｝を生成する手段と、
該手段が生成した白色雑音の分散σe2を計算する手段と、
下記式で定義される行列Ｔ及びｔを計算する手段と、
【数５】
2005242580.tif 000035
前記ＡＲＭＡＸモデルにおける前記標本時系列に係るｐ個のパラメータの推定値を行列要素とするベクトルａ＝［ａ1 ，ａ2 ，…，ａp ］T、前記入力時系列に係るｐ個のパラメータの推定値を行列要素とする行列ｃ＝［ｃ1 ，ｃ2 ，…，ｃp］T 及び前記白色雑音に係るｐ個のパラメータの推定値を行列要素とする行列ｂ＝［ｂ1 ，ｂ2 ，…，ｂp］T を、下記式に基づいて計算する手段と、
【数６】
2005242580.tif 000036
該手段が計算したベクトルａ、ｂ及びｃに含まれる前記ＡＲＭＡＸモデルのパラメータの推定値、並びに前記標本時系列及び前記入力時系列値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算する手段と、
該手段が計算した新たな白色雑音を用いて新たなベクトルａ、ｂ及びｃを計算する手段と、
ベクトルａ、ｂ及びｃが収束するまでベクトルａ、ｂ及びｃの計算を繰り返す手段と
を備えることを特徴とするパラメータ推定装置。
【請求項４】
入力を伴う系から得られる複数のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…｝、及び前記系へ入力される複数のデータからなる入力時系列｛ｕn：ｎ＝１，２，…｝から、前記標本時系列に加わるべき新たなデータを予測する装置において、
前記標本時系列及び前記入力時系列を記憶する手段と、
請求項３に記載のパラメータ推定装置と、
該パラメータ推定装置を用いて、前記標本時系列及び前記入力時系列をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を計算する手段と、
該手段が計算した前記パラメータの推定値、並びに前記標本時系列及び前記入力時系列を前記ＡＲＭＡモデルの式に代入して、前記標本時系列に加わるべき新たなデータの予測値を計算する手段と
を備えることを特徴とするデータ予測装置。
【請求項５】
複数のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…｝及び入力時系列｛ｕn ：ｎ＝１，２，…｝、並びに下記式で定義される行列Ｒ，Ｗ，Ｖ，ｒ，ｖを記憶するコンピュータに、
【数７】
2005242580.tif 000037
前記標本時系列及び前記入力時系列を（ｐ，ｐ，ｐ）次（但し、ｐは自然数）のＡＲＭＡＸモデルに当てはめたときのパラメータを推定させるコンピュータプログラムであって、
コンピュータに、平均値が０である複数のランダムなデータからなる白色雑音｛ｅn ：ｎ＝１，２，…｝を生成させる手順と、
コンピュータに、生成した白色雑音の分散σe2を計算させる手順と、
コンピュータに、下記式で定義される行列Ｔ及びｔを計算させる手順と、
【数８】
2005242580.tif 000038
コンピュータに、前記ＡＲＭＡＸモデルにおける前記標本時系列に係るｐ個のパラメータの推定値を行列要素とするベクトルａ＝［ａ1 ，ａ2 ，…，ａp ］T、前記入力時系列に係るｐ個のパラメータの推定値を行列要素とする行列ｃ＝［ｃ1 ，ｃ2 ，…，ｃp］T 及び前記白色雑音に係るｐ個のパラメータの推定値を行列要素とする行列ｂ＝［ｂ1 ，ｂ2 ，…，ｂp］T を、下記式に基づいて計算させる手順と、
【数９】
2005242580.tif 000039
コンピュータに、計算したベクトルａ、ｂ及びｃに含まれる前記ＡＲＭＡＸモデルのパラメータの推定値、並びに前記標本時系列及び前記入力時系列値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算させる手順と、
コンピュータに、計算した新たな白色雑音を用いて新たなベクトルａ、ｂ及びｃを計算させる手順と、
コンピュータに、ベクトルａ、ｂ及びｃが収束するまでベクトルａ、ｂ及びｃの計算を繰り返させる手順と
を含むことを特徴とするコンピュータプログラム。
【請求項６】
入力を伴う系から得られる複数のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…｝、及び前記系へ入力される複数のデータからなる入力時系列｛ｕn：ｎ＝１，２，…｝を記憶するコンピュータに、前記標本時系列に加わるべき新たなデータを予測させるコンピュータプログラムであって、
コンピュータに、請求項１に記載のパラメータ推定方法を用いて、前記標本時系列及び前記入力時系列をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を計算させる手順と、
コンピュータに、計算した前記パラメータの推定値、並びに前記標本時系列及び前記入力時系列を前記ＡＲＭＡモデルの式に代入して、前記標本時系列に加わるべき新たなデータの予測値を計算させる手順と
を含むことを特徴とするコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、データの時系列からパラメータを推定する方法、推定されたパラメータを用いて新たなデータを予測する方法、パラメータ推定装置、データ予測装置、コンピュータをパラメータ推定装置又はデータ予測装置として実現するためのコンピュータプログラムに関する。
【背景技術】
【０００２】
あるシステムから得られるデータの時系列を解析して新たなデータの値を予想する方法は、従来様々な方法が提案されている。その方法の一つとして、データの時系列を離散時間線形モデルに当てはめて解析を行い、データ間の関係を示すパラメータを推定し、推定したパラメータを用いて構成される離散時間線形モデルに基づいて新たなデータの値を予測する方法がある。離散時間線形モデルとしては、ＡＲ（Autoregressive：自己回帰）モデル又はＡＲＭＡ（AutoregressiveMoving Average ：自己回帰移動平均）モデル等が利用されている。
【０００３】
データの時系列｛ｘn ：ｎ＝１，２，…｝を、定常エルゴード的正規過程ｘ（ｔ）の標本時系列であるとすると、（ｐ，ｑ）次のＡＲＭＡモデルは下記（１）式で表される。
【０００４】
【数１】
2005242580.tif 000003
【０００５】
ここで、｛ｅn ｝は、平均値０及び分散σe2の正規性白色雑音である。また、｛ａi ：i＝１，２，…，ｐ｝及び｛ｂj ：ｊ＝１，２，…，ｑ｝はＡＲＭＡモデルのパラメータである。また、｛ａi ｝及び｛ｂj ｝は、下記式で与えられるＡ（Ｚ-1）及びＢ（Ｚ-1）が規約で夫々定常条件、可逆条件および強正実条件を満たすものとする。
Ａ（Ｚ-1）＝１＋ａ1 Ｚ-1＋…＋ａp Ｚ-p 
Ｂ（Ｚ-1）＝１＋ｂ1 Ｚ-1＋…＋ｂq Ｚ-q
【０００６】
標本の時系列｛ｘn ｝から（１）式のパラメータ｛ａi｝及び｛ｂj ｝を推定するための方法として、疑似線形回帰法が従来用いられている。従来の疑似線形回帰法では、白色雑音｛ｅn｝を乱数などを用いて生成し、生成した｛ｅn ｝と時系列｛ｘn ｝とから、｛ａi｝及び｛ｂj ｝を計算し、計算した｛ａi ｝及び｛ｂj｝から（１）式を用いて白色雑音｛ｅn ｝を計算し、パラメータ｛ａi｝及び｛ｂj ｝が収束するまで計算を繰り返す。この方法によりＡＲＭＡモデルのパラメータが推定される。非特許文献１には、ＡＲＭＡモデルのパラメータ｛ａi ｝及び｛ｂj ｝の値の収束を容易にする技術が提案されている。
【０００７】
以上の如きＡＲＭＡモデルでは、データの時系列に含まれる一のデータの値は過去のデータの値と過去の雑音とから決定される。しかし現実には、商品の販売量のデータが気温に影響される場合などのように、システムから得られるデータがシステムへの入力に影響されていることも多い。従って、ＡＲＭＡモデルでは実際の問題に対処しきれないという問題がある。
【０００８】
システムへの入力を考慮した離散時間線形モデルとして、ＡＲＭＡＸ（AutoregressiveMoving  Average with exogenous input：外生入力自己回帰移動平均）モデルが従来用いられている。システムから得られた複数のデータからなる標本時系列｛ｘn：ｎ＝１，２，…｝を定常エルゴード的正規過程ｘ（ｔ）の標本時系列であるとすると、（ｐ，ｍ，ｑ）次のＡＲＭＡＸモデルは、下記（２）式で表される。
【０００９】
【数２】
2005242580.tif 000004
【００１０】
ここで、｛ｅn ｝は平均値０及び分散σe2の正規性白色雑音であり、｛ｕn ｝はシステムに対する入力時系列である。また｛ａi ：i＝１，２，…，ｐ｝、｛ｃi ：ｉ＝１，２，…，ｍ｝及び｛ｂi：ｉ＝１，２，…，ｑ｝はＡＲＭＡＸモデルのパラメータである。また、（２）式は定常条件を満たすとする。現代制御理論では、ｑ≦ｐ，ｍ≦ｐとなることは周知である。以下、簡単にするためにｑ＝ｍ＝ｐとして議論を進める。
ベクトルＺn 及びベクトルθを下記のように定義する。
Ｚn ＝［-ｘn-1,…,-ｘn-p,ｕn-1,…,ｕn-p,ｅn-1,…,ｅn-p］T
θ＝［ａ1,…,ａp,ｃ1,…,ｃp,ｂ1,…,ｂp］T
（２）式をｘn についての式に書き直すと、（２）式は下記（３）式で表すことができる。
ｘn ＝θT Ｚn＋ｅn   …（３）
【００１１】
Ｎ個のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…，Ｎ｝及び入力時系列｛ｕn：ｎ＝１，２，…，Ｎ｝が得られたとき、ベクトルθを時系列から最小二乗法により求めるためには、（３）式から得られる下記（４）式の右辺を最小にするθを求めればよい。
【００１２】
【数３】
2005242580.tif 000005
【００１３】
従って、（４）式の右辺をθで偏微分して０に等しいとすることにより、Ｎ個のデータからなる標本時系列｛ｘn ｝及び入力時系列｛ｕn ｝から推定したベクトルθの推定値θN は、下記（５）式で与えられる。
【００１４】
【数４】
2005242580.tif 000006
【００１５】
ｋ≧０としてラグがｋである｛ｘn ｝の自己相関関数をＲkとすると、自己相関関数の定義によりＲk 及びＲ-kは下記式で表すことができる。ここで、Ｅ［ｙn］は、データ列｛ｙn ：ｎ＝１，２，…｝の期待値である。
Ｒk ＝Ｅ［ｘn ｘn+k ］，  Ｒ-k＝Ｒk
【００１６】
同様に、｛ｕn ｝の自己相関関数をＷk とすると、Ｗk 及びＷ-kは下記式で表すことができる。
Ｗk ＝Ｅ［ｕn ｕn+k ］，  Ｗ-k＝Ｗk
【００１７】
同様に、｛ｅn ｝の自己相関関数をＳｋとすると、Ｓｋ及びＳ-kは下記式で表すことができる。
Ｓｋ＝Ｅ［ｅn ｅn+k ］，  Ｓ-k＝Ｓk
【００１８】
｛ｘn ｝と｛ｅn｝との相互相関関数については、Ｔk 及びＴ-kを下記式で定義する。
Ｔk ＝Ｅ［ｅn ｘn+k ］，  Ｔ-k＝Ｅ［ｅn ｘn-k］
このとき、｛ｘn ｝と｛ｅn｝との夫々の相互相関関数は下記式で表される。
Ｅ［ｘn ｅn-k ］＝Ｅ［ｅnｘn+k ］＝Ｔk 
Ｅ［ｘn ｅn+k ］＝Ｅ［ｅn ｘn-k ］＝Ｔ-k
【００１９】
｛ｘn ｝と｛ｕn｝との相互相関関数についても、Ｖk 及びＶ-kを下記式で定義する。
Ｖk ＝Ｅ［ｕnｘn+k ］，  Ｖ-k＝Ｅ［ｕn ｘn-k］
このとき、｛ｘn ｝と｛ｕn｝との夫々の相互相関関数は下記式で表される。
Ｅ［ｘn ｕn-k ］＝Ｅ［ｕnｘn+k ］＝Ｖk 
Ｅ［ｘn ｕn+k ］＝Ｅ［ｕn ｘn-k ］＝Ｖ-k
【００２０】
｛ｕn ｝と｛ｅn｝との相互相関関数についても、Ｑk 及びＱ-kを下記式で定義する。
Ｑk ＝Ｅ［ｅn ｕn+k ］，  Ｑ-k＝Ｅ［ｕn ｅn+k ］
このとき、｛ｕn ｝と｛ｅn｝との夫々の相互相関関数は下記式で表される。
Ｅ［ｕn ｅn-k ］＝Ｅ［ｅnｕn+k ］＝Ｑk 
Ｅ［ｕn ｅn+k ］＝Ｅ［ｅn ｕn-k ］＝Ｑ-k
【００２１】
ベクトルａ，ｂ，ｃを下記式で定義する。
ａ＝［ａ1 ，ａ2 ，…，ａp］T
ｂ＝［ｂ1 ，ｂ2 ，…，ｂp］T
ｃ＝［ｃ1 ，ｃ2 ，…，ｃp］T
【００２２】
（５）式においてＮ→∞とすると、（５）式の左辺はθN →θとなり、（５）式の右辺は｛ｘn ｝，｛ｕn ｝，｛ｅn ｝の自己相関関数及び相互相関関数で記述できる。従って、（５）式は下記（６）式で表すことができる。
【００２３】
【数５】
2005242580.tif 000007
【００２４】
但し、Ｒ，Ｖ’，Ｔ’，Ｗ，Ｓ’，Ｑ’はｐ行ｐ列の行列であって、夫々下記（７），，（８），（９），（１０），（１１），（１２）式で表される。また、ｒ，ｖ，ｔはｐ行一列の行列であって、下記（１３）式で表される。
【００２５】
【数６】
2005242580.tif 000008
【００２６】
【数７】
2005242580.tif 000009
【００２７】
【数８】
2005242580.tif 000010
【００２８】
【数９】
2005242580.tif 000011
【００２９】
【数１０】
2005242580.tif 000012
【００３０】
【数１１】
2005242580.tif 000013
【００３１】
【数１２】
2005242580.tif 000014
【００３２】
（６）式を書き直すと、下記の連立方程式が得られる。
Ｒａ−Ｖ’ｃ−Ｔ’ｂ＝−ｒ              …（１４）
−Ｖ’T ａ＋Ｗｃ＋Ｑ’ｂ＝ｖ     …（１５）
−Ｔ’T ａ＋Ｑ’T ｃ＋Ｓ’ｂ＝ｔ …（１６）
この連立方程式（１４）〜（１６）を解くことにより、ベクトルａ，ｂ，ｃを計算することができる。
【００３３】
疑似線形回帰法を用いてＡＲＭＡＸモデルのパラメータの推定値であるベクトルａ，ｂ，ｃを求めるためには、白色雑音｛ｅn ｝を乱数などを用いて生成し、生成した｛ｅn｝と｛ｘn ｝及び｛ｕn｝とから、ベクトルａ，ｂ，ｃを（１４）〜（１６）式により計算し、計算したベクトルａ，ｂ，ｃから（２）式を用いて新たな白色雑音｛ｅn ｝を計算し、ベクトルａ，ｂ，ｃが収束するまで計算を繰り返す。
【非特許文献１】竹安数博、外３名，「ロバストなブートストラップ型システム同定アルゴリズム」，Journal of Signal Processing，信号処理学会，２００３年３月，第７巻，第２号，ｐ．１６７−１７６
【発明の開示】
【発明が解決しようとする課題】
【００３４】
前述の如き疑似線形回帰法は、問題の性質によっては、求めるべきパラメータであるベクトルａ，ｂ，ｃが収束せず、パラメータの推定ができない場合があるという問題がある。また、前述の方法では、ベクトルａ，ｂ，ｃを計算するためにｐ行ｐ列の行列Ｒ，Ｖ’，Ｔ’，Ｗ，Ｓ’，Ｑ’を用いた行列計算を行う必要があるため、計算時間が長いという問題があった。
【００３５】
本発明は、斯かる事情に鑑みてなされたものであって、その目的とするところは、外部入力を考慮に入れたＡＲＭＡＸモデルのパラメータの値の収束を容易にし、同時に計算時間の短縮を図ることができるパラメータ推定方法、データ予測方法、パラメータ推定装置、データ予測装置、及びコンピュータプログラムを提供することにある。
【課題を解決するための手段】
【００３６】
第１発明に係るパラメータ推定方法は、記憶部及び演算部を備えたコンピュータを用いて、複数のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…｝及び入力時系列｛ｕn：ｎ＝１，２，…｝を（ｐ，ｐ，ｐ）次（但し、ｐは自然数）のＡＲＭＡＸモデルに当てはめたときのパラメータを推定する方法において、前記標本時系列及び前記入力時系列を前記記憶部で記憶し、下記式で定義される行列Ｒ，Ｗ，Ｖ，ｒ，ｖを前記記憶部で記憶し、
【００３７】
【数１３】
2005242580.tif 000015
【００３８】
平均値が０である複数のランダムなデータからなる白色雑音｛ｅn ：ｎ＝１，２，…｝を前記演算部で生成し、前記演算部で生成した白色雑音の分散σe2を前記演算部で計算し、下記式で定義される行列Ｔ及びｔを前記演算部で計算し、
【００３９】
【数１４】
2005242580.tif 000016
【００４０】
前記ＡＲＭＡＸモデルにおける前記標本時系列に係るｐ個のパラメータを行列要素とするベクトルａ＝［ａ1 ，ａ2 ，…，ａp ］T、前記入力時系列に係るｐ個のパラメータを行列要素とする行列ｃ＝［ｃ1 ，ｃ2 ，…，ｃp ］T及び前記白色雑音に係るｐ個のパラメータを行列要素とする行列ｂ＝［ｂ1 ，ｂ2 ，…，ｂp ］Tを、下記式に基づいて前記演算部で計算し、
【００４１】
【数１５】
2005242580.tif 000017
【００４２】
前記演算部で計算したベクトルａ、ｂ及びｃに含まれる前記ＡＲＭＡＸモデルのパラメータ、並びに前記標本時系列及び前記入力時系列値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を前記演算部で計算し、前記演算部で計算した新たな白色雑音を用いて新たなベクトルａ、ｂ及びｃを前記演算部で計算し、ベクトルａ、ｂ及びｃが収束するまでベクトルａ、ｂ及びｃの計算を前記演算部で繰り返すことを特徴とする。
【００４３】
第２発明に係るデータ予測方法は、記憶部及び演算部を備えたコンピュータを用いて、入力を伴う系から得られる複数のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…｝、及び前記系へ入力される複数のデータからなる入力時系列｛ｕn：ｎ＝１，２，…｝から、前記標本時系列に加わるべき新たなデータを予測する方法において、前記標本時系列及び前記入力時系列を前記記憶部で記憶し、第１発明に係るパラメータ推定方法を用いて、前記標本時系列及び前記入力時系列をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を前記演算部で計算し、前記演算部で計算した前記パラメータの推定値、並びに前記標本時系列及び前記入力時系列を前記ＡＲＭＡモデルの式に代入して、前記標本時系列に加わるべき新たなデータの予測値を前記演算部で計算することを特徴とする。
【００４４】
第３発明に係るパラメータ推定装置は、複数のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…｝及び入力時系列｛ｕn：ｎ＝１，２，…｝を（ｐ，ｐ，ｐ）次（但し、ｐは自然数）のＡＲＭＡＸモデルに当てはめたときのパラメータを推定する装置において、前記標本時系列及び前記入力時系列を記憶する手段と、下記式で定義される行列Ｒ，Ｗ，Ｖ，ｒ，ｖを記憶する手段と、
【００４５】
【数１６】
2005242580.tif 000018
【００４６】
平均値が０である複数のランダムなデータからなる白色雑音｛ｅn ：ｎ＝１，２，…｝を生成する手段と、該手段が生成した白色雑音の分散σe2を計算する手段と、下記式で定義される行列Ｔ及びｔを計算する手段と、
【００４７】
【数１７】
2005242580.tif 000019
【００４８】
前記ＡＲＭＡＸモデルにおける前記標本時系列に係るｐ個のパラメータの推定値を行列要素とするベクトルａ＝［ａ1 ，ａ2 ，…，ａp ］T、前記入力時系列に係るｐ個のパラメータの推定値を行列要素とする行列ｃ＝［ｃ1 ，ｃ2 ，…，ｃp］T 及び前記白色雑音に係るｐ個のパラメータの推定値を行列要素とする行列ｂ＝［ｂ1 ，ｂ2 ，…，ｂp］T を、下記式に基づいて計算する手段と、
【００４９】
【数１８】
2005242580.tif 000020
【００５０】
該手段が計算したベクトルａ、ｂ及びｃに含まれる前記ＡＲＭＡＸモデルのパラメータの推定値、並びに前記標本時系列及び前記入力時系列値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算する手段と、該手段が計算した新たな白色雑音を用いて新たなベクトルａ、ｂ及びｃを計算する手段と、ベクトルａ、ｂ及びｃが収束するまでベクトルａ、ｂ及びｃの計算を繰り返す手段とを備えることを特徴とする。
【００５１】
第４発明に係るデータ予測装置は、入力を伴う系から得られる複数のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…｝、及び前記系へ入力される複数のデータからなる入力時系列｛ｕn：ｎ＝１，２，…｝から、前記標本時系列に加わるべき新たなデータを予測する装置において、前記標本時系列及び前記入力時系列を記憶する手段と、第３発明に係るパラメータ推定装置と、該パラメータ推定装置を用いて、前記標本時系列及び前記入力時系列をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を計算する手段と、該手段が計算した前記パラメータの推定値、並びに前記標本時系列及び前記入力時系列を前記ＡＲＭＡモデルの式に代入して、前記標本時系列に加わるべき新たなデータの予測値を計算する手段とを備えることを特徴とする。
【００５２】
第５発明に係るコンピュータプログラムは、複数のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…｝及び入力時系列｛ｕn：ｎ＝１，２，…｝、並びに下記式で定義される行列Ｒ，Ｗ，Ｖ，ｒ，ｖを記憶するコンピュータに、
【００５３】
【数１９】
2005242580.tif 000021
【００５４】
前記標本時系列及び前記入力時系列を（ｐ，ｐ，ｐ）次（但し、ｐは自然数）のＡＲＭＡＸモデルに当てはめたときのパラメータを推定させるコンピュータプログラムであって、コンピュータに、平均値が０である複数のランダムなデータからなる白色雑音｛ｅn ：ｎ＝１，２，…｝を生成させる手順と、コンピュータに、生成した白色雑音の分散σe2を計算させる手順と、コンピュータに、下記式で定義される行列Ｔ及びｔを計算させる手順と、
【００５５】
【数２０】
2005242580.tif 000022
【００５６】
コンピュータに、前記ＡＲＭＡＸモデルにおける前記標本時系列に係るｐ個のパラメータの推定値を行列要素とするベクトルａ＝［ａ1 ，ａ2 ，…，ａp ］T、前記入力時系列に係るｐ個のパラメータの推定値を行列要素とする行列ｃ＝［ｃ1 ，ｃ2 ，…，ｃp］T 及び前記白色雑音に係るｐ個のパラメータの推定値を行列要素とする行列ｂ＝［ｂ1 ，ｂ2 ，…，ｂp］T を、下記式に基づいて計算させる手順と、
【００５７】
【数２１】
2005242580.tif 000023
【００５８】
コンピュータに、計算したベクトルａ、ｂ及びｃに含まれる前記ＡＲＭＡＸモデルのパラメータの推定値、並びに前記標本時系列及び前記入力時系列値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算させる手順と、コンピュータに、計算した新たな白色雑音を用いて新たなベクトルａ、ｂ及びｃを計算させる手順と、コンピュータに、ベクトルａ、ｂ及びｃが収束するまでベクトルａ、ｂ及びｃの計算を繰り返させる手順とを含むことを特徴とする。
【００５９】
第６発明に係るコンピュータプログラムは、入力を伴う系から得られる複数のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…｝、及び前記系へ入力される複数のデータからなる入力時系列｛ｕn：ｎ＝１，２，…｝を記憶するコンピュータに、前記標本時系列に加わるべき新たなデータを予測させるコンピュータプログラムであって、コンピュータに、第１発明に係るパラメータ推定方法を用いて、前記標本時系列及び前記入力時系列をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を計算させる手順と、コンピュータに、計算した前記パラメータの推定値、並びに前記標本時系列及び前記入力時系列を前記ＡＲＭＡモデルの式に代入して、前記標本時系列に加わるべき新たなデータの予測値を計算させる手順とを含むことを特徴とする。
【００６０】
本発明者は、ＡＲＭＡＸモデルの式の構成から導かれるa priori knowledgeを導入することにより、従来の疑似線形回帰法の改良を試みた。まず正規性雑音｛ｅn ｝の自己相関関数を考える。正規性白色雑音｛ｅn｝に含まれる各要素は互いに関連性がないので、｛ｅn ｝の自己相関関数Ｓｋは理論的に下記式で表すことができる。
【００６１】
【数２２】
2005242580.tif 000024
【００６２】
従って、（１１）式で定義される行列Ｓ’に理論的な値を導入することによって行列Ｓ’を改良した行列Ｓは、単位行列をＩとして下記式で表される。
【００６３】
【数２３】
2005242580.tif 000025
【００６４】
次に、｛ｘn ｝と｛ｅn｝との相互相関関数を考える。（２）式で示されるＡＲＭＡＸモデルでは、現在のｘn の値は過去のｅn の値に関連づけられている。しかし、現在のｘnの値は未来のｅn の値とは関連性がないので、ｘn とこのデータよりも未来のｅn との相互相関関数は下記式で表すことができる。
Ｅ［ｘn ｅn+k ］＝０
同様に下記式が成り立つ。
Ｅ［ｅn ｘn-k ］＝０
従って、Ｔ-k＝０となる。また（２）式より現在のｘn の値と現在のｅn の値とは関連性があるので、Ｔ0≠０である。従って、（９）式で定義される行列Ｔ’に理論的な値を導入することによって行列Ｔ’を改良した行列Ｔは下記式で表される。
【００６５】
【数２４】
2005242580.tif 000026
【００６６】
｛ｘn ｝と｛ｕn ｝との相互相関関数についても、同様に、現在のｘnの値は過去のｕn の値に関連づけられているが、現在のｘn の値は未来のｕn の値とは関連性がないので、下記式が成り立つ。
Ｅ［ｘn ｕn+k ］＝Ｅ［ｕn ｘn-k ］＝Ｖ-k＝０
また（２）式より現在のｘn の値と現在のｕn の値とは関連性がないので、Ｖ0 ＝０である。従って、（８）式で定義される行列Ｖ’に理論的な値を導入することによって行列Ｖ’を改良した行列Ｖは下記式で表される。
【００６７】
【数２５】
2005242580.tif 000027
【００６８】
更に、入力時系列｛ｕn ｝と正規性白色雑音｛ｅn ｝とは無相関であるはずであるので、下記式に示す如く｛ｕn ｝と｛ｅn｝との相互相関関数は全て０となる。
Ｅ［ｅn ｕn+k ］＝Ｅ［ｕn ｅn-k ］＝Ｑk＝０
Ｅ［ｅn ｕn-k ］＝Ｅ［ｕnｅn+k ］＝Ｑ-k＝０
従って、（１２）式で定義される行列Ｑ’に理論的な値を導入することによって行列Ｑ’を改良した行列Ｑは、ゼロ行列をＯとして下記式で表される。
Ｑ＝Ｏ      …（２１）
【００６９】
Ｎ個のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…，Ｎ｝及び入力時系列｛ｕn：ｎ＝１，２，…，Ｎ｝が得られたときに各自己相関関数及び各相互相関関数を計算した場合は、Ｎの値が有限であるので相互相関関数及び自己相関関数の値は理論的な値とは違った値となる。この実際に計算される値の内、理論的な値が得られる部分は理論的な値と置き換えることにより、ＡＲＭＡＸモデルのパラメータを求めるための（６）式は下記式で表すことができるようになる。
【００７０】
【数２６】
2005242580.tif 000028
【００７１】
（２２）式を書き直すと、下記の連立方程式が得られる。
Ｒａ−Ｖｃ−Ｔｂ＝−ｒ    …（２３）
−ＶT ａ＋Ｗｃ＝ｖ        …（２４）
−ＴT ａ＋σe2ｂ＝ｔ             …（２５）
（２５）式及び（２４）より、下記の（２６）式及び（２７）式が得られる。
ｂ＝（ｔ＋ＴT ａ）／σe2  …（２６）
ｃ＝Ｗ-1（ｖ＋ＶT ａ）       …（２７）
（２６）式及び（２７）式を（２３）式に代入することにより、下記（２８）式が得られる。
【００７２】
【数２７】
2005242580.tif 000029
【００７３】
（２８）式よりベクトルａを計算することができる。計算したベクトルａを（２６）式及び（２７）式へ代入することにより、ベクトルｂ及びベクトルｃを計算することができる。ＡＲＭＡＸモデルのパラメータの推定値であるベクトルａ，ｂ，ｃを求めるためには、白色雑音｛ｅn ｝を乱数などを用いて生成し、生成した｛ｅn｝と｛ｘn ｝及び｛ｕn｝とから、ベクトルａ，ｂ，ｃを（２６）〜（２８）式により計算し、計算したベクトルａ，ｂ，ｃから（２）式を用いて新たな白色雑音｛ｅn ｝を計算し、ベクトルａ，ｂ，ｃが収束するまで計算を繰り返す。
【００７４】
第１、第３及び第５発明においては、入力を伴う系から得られた複数のデータからなる標本時系列｛ｘn ｝及び系へ入力された複数のデータからなる入力時系列｛ｕn ｝をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を、入力時系列の自己相関関数、白色雑音の自己相関関数、標本時系列と入力時系列との相互相関関数、標本時系列と白色雑音の相互相関関数、及び入力時系列と白色雑音との相互相関関数に理論的な値を用いた（２６）〜（２８）式により計算する。
【００７５】
第２、第４及び第６発明においては、（２６）〜（２８）式により計算したＡＲＭＡＸモデルのパラメータの推定値を用い、標本時系列｛ｘn ｝に加わるべき新たなデータの予測値を計算する。
【発明の効果】
【００７６】
第１、第３及び第５発明にあっては、入力を伴う系から得られた複数のデータからなる標本時系列｛ｘn ｝及び系へ入力された複数のデータからなる入力時系列｛ｕn ｝をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を計算する際に、白色雑音の自己相関関数、標本時系列と入力時系列との相互相関関数、標本時系列と白色雑音の相互相関関数、及び入力時系列と白色雑音との相互相関関数に理論的な値を用いて計算することにより、従来の疑似線形回帰法でパラメータの推定値の収束が困難であった場合についても、パラメータの推定値の収束が容易となる。またパラメータが収束する確率が向上するので、パラメータ推定の計算時間を短縮する事ができる。
【００７７】
第２、第４及び第６発明にあっては、系から得られるデータの値が過去のデータの値と雑音と更に系への入力とから決定されるＡＲＭＡＸモデルのパラメータの推定値を用いて新たなデータの予測値を計算することにより、株の取引量に影響される株価の値、又は商品の入荷量に影響される商品の販売量など、外部入力に影響される系から得られるデータの値を推定することができる。従って、株価又は商品の販売量などの予測をより現実に即して行うことができる等、本発明は優れた効果を奏する。
【発明を実施するための最良の形態】
【００７８】
以下本発明をその実施の形態を示す図面に基づき具体的に説明する。
図１は、本発明のデータ予測装置の構成を示すブロック図である。データ予測装置１は、本発明のパラメータ推定装置としての機能を兼ね備えており、汎用コンピュータを用いて構成されている。データ予測装置１は、演算を行うＣＰＵ（演算部）１１と、演算に伴って発生する一時的な情報を記憶するＲＡＭ（記憶部）１２と、ＣＤ−ＲＯＭドライブ等の外部記憶装置１３と、ハードディスク等の内部記憶装置１４とを備えている。ＣＰＵ１１は、ＣＤ−ＲＯＭ等の記録媒体２から本発明のコンピュータプログラム２０を外部記憶装置１３にて読み取り、読み取ったコンピュータプログラム２０を内部記憶装置１４に記憶させる。コンピュータプログラム２０は必要に応じて内部記憶装置１４からＲＡＭ１２へロードされ、ロードされたコンピュータプログラム２０に基づいてＣＰＵ１１はデータ予測装置１に必要な処理を実行する。データ予測装置１は、キーボード又はマウス等の入力装置１５と、液晶ディスプレイ又はＣＲＴディスプレイ等の出力装置１６とを備えており、データの入力を初めとするオペレータからの操作を受け付ける構成となっている。
【００７９】
内部記憶装置１４は、株価又は商品の販売量などの複数のデータからなる標本時系列｛ｘn ：ｎ＝１，２，…｝と、株の取引量又は商品の入荷量など、標本時系列が得られる系に対して入力される複数のデータからなる入力時系列｛ｕn：ｎ＝１，２，…｝とを記憶している。
【００８０】
なお、コンピュータプログラム２０は、通信ネットワークＮに接続されている図示しない外部のサーバ装置から本発明に係る状態監視装置１へロードされて内部記憶装置１４に記憶される形態であってもよい。
【００８１】
図２は、本発明のデータ予測装置１が行うパラメータ推定の処理およびデータ予測の処理を示すフローチャートである。データ予測装置１のＣＰＵ１１は、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って以下の処理を行う。データ予測装置１のＣＰＵ１１は、内部記憶装置１４が記憶している標本時系列｛ｘn ｝及び入力時系列｛ｕn ｝をＲＡＭ１２に読み出す（Ｓ１）。ＣＰＵ１１は、次に、｛ｘn ｝及び｛ｕn｝に含まれるデータの値を用いて｛ｘn ｝の自己相関関数Ｒk ＝Ｅ［ｘn ｘn+k ］，｛ｕn ｝の自己相関関数Ｗk ＝Ｅ［ｕnｕn+k ］，｛ｘn ｝と｛ｕn ｝との相互相関関数Ｖk＝Ｅ［ｕn ｘn+k ］を計算して、（７）式，（１０）式，（２０）式，（１３）式で定義される行列Ｒ，Ｗ，Ｖ，ｒ，ｖを計算する（Ｓ２）。ＣＰＵ１１は、次に、乱数生成のルーチンを用いて時刻の値を初期値として乱数を発生させる等の方法により、平均値が０である複数のランダムなデータからなる白色雑音｛ｅn：ｎ＝１，２，…｝を生成し（Ｓ３）、生成した白色雑音｛ｅn ｝の分散σe2を計算する（Ｓ４）。
【００８２】
ＣＰＵ１１は、次に、｛ｘn ｝及び｛ｅn ｝に含まれるデータの値を用いて｛ｘn ｝と｛ｅn ｝との相互相関関数Ｔk＝Ｅ［ｅn ｘn+k ］を計算して、（１９）式，（１３）式で定義される行列Ｔ，ｔを計算する（Ｓ５）。ＣＰＵ１１は、次に、（２８）式によりベクトルａを計算し、ベクトルａを（２６）式及び（２７）式へ代入することによりベクトルｂ及びベクトルｃを計算することによって、｛ｘn ｝及び｛ｕn ｝をＡＲＭＡＸモデルに当てはめたときのＡＲＭＡＸモデルのパラメータの推定値からなるベクトルａ，ｂ，ｃを計算する（Ｓ６）。
【００８３】
ＣＰＵ１１は、次に、計算したベクトルａ，ｂ，ｃが収束しているか否かを判定する（Ｓ７）。このときＣＰＵ１１は、例えば、最新のベクトルａ，ｂ，ｃと前回に計算したベクトルａ，ｂ，ｃであるベクトルａ’＝［ａ’1 ，ａ’2 ，…，ａ’p ］T，ｂ’＝［ｂ’1 ，ｂ’2 ，…，ｂ’p ］T ，ｃ’＝［ｃ’1，ｃ’2 ，…，ｃ’p ］T を用いて下記（２９）式で定義されるＤを計算し、計算したＤが所定の微少量εに対してＤ＜εである状態が所定回数連続で成立している場合にベクトルａ，ｂ，ｃが収束していると判定する。なお、その他の方法を用いてベクトルａ，ｂ，ｃの収束を判定する処理を用いてもよい。
【００８４】
【数２８】
2005242580.tif 000030
【００８５】
ステップＳ７にて、ベクトルａ，ｂ，ｃが収束していないと判定された場合は（Ｓ７：ＮＯ）、ベクトルａ，ｂ，ｃと｛ｘn ｝及び｛ｕn ｝とをｍ＝ｑ＝ｐとした（２）式へ代入して、新たな白色雑音｛ｅn：ｎ＝１，２，…｝を計算し（Ｓ８）、処理をステップＳ４へ戻す。ＣＰＵ１１は、ステップＳ７にてベクトルａ，ｂ，ｃが収束したと判定されるまでステップＳ４〜Ｓ８の処理を繰り返す。
【００８６】
ステップＳ７にて、ベクトルａ，ｂ，ｃが収束していると判定された場合は（Ｓ７：ＹＥＳ）、ＣＰＵ１１は、計算したベクトルａ，ｂ，ｃ、｛ｘn ｝、｛ｕn ｝、｛ｅn ｝を（２）式へ代入することにより、標本時系列｛ｘn｝に加わるべきデータの予測値、即ち株価又は商品の販売量などの予測データを計算する（Ｓ９）。なお、ステップＳ９では、入力時系列に加わる新たなデータ又は当該データの予測値を必要に応じて用いてもよい。以上の処理により、ＡＲＭＡＸモデルのパラメータが推定され、標本時系列に加わるべき新たなデータが予測できる。
【００８７】
以上詳述した如く、本発明においては、入力を伴う系から得られた複数のデータからなる標本時系列｛ｘn ｝及び系へ入力された複数のデータからなる入力時系列｛ｕn ｝をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を計算する際に、白色雑音の自己相関関数、標本時系列と入力時系列との相互相関関数、標本時系列と白色雑音の相互相関関数、及び入力時系列と白色雑音との相互相関関数に理論的な値を用いて計算する。これにより、従来の疑似線形回帰法でパラメータの推定値の収束が困難であった場合についても、パラメータの推定値の収束が容易となる。またパラメータが収束する確率が向上するので、パラメータ推定の計算時間を短縮する事ができる。
【００８８】
また、本発明においては、系から得られるデータの値が過去のデータの値と雑音と更に系への入力とから決定されるＡＲＭＡＸモデルを離散時間線形モデルとして用いることにより、株の取引量に影響される株価の値、又は商品の入荷量に影響される商品の販売量など、外部入力に影響される系から得られるデータの値を推定することができる。従って、株価又は商品の販売量などの予測をより現実に即して行うことができる。
【００８９】
なお、本実施の形態においては、本発明のデータ予測装置は、内部で標本時系列｛ｘn ｝及び入力時系列｛ｕn ｝を記憶している形態を示しているが、これに限るものではなく、標本時系列｛ｘn ｝及び入力時系列｛ｕn ｝に含まれるデータを外部から受け付ける手段を備え、受け付けたデータについて予測データを計算する処理を行う形態であってもよい。
【００９０】
また、本実施の形態においては、本発明のデータ予測装置１は、本発明のパラメータ推定装置の機能を兼ね備えている形態を示したが、これに限るものではなく、ＡＲＭＡＸモデルのパラメータの予測値を計算するパラメータ推定装置とＡＲＭＡＸモデルによりデータを予測するデータ予測装置とを夫々別のコンピュータにて構成し、互いに情報を入出力してパラメータの推定とデータの予測とを行う形態としても良い。
【図面の簡単な説明】
【００９１】
【図１】本発明のデータ予測装置の構成を示すブロック図である。
【図２】本発明のデータ予測装置が行うパラメータ推定の処理およびデータ予測の処理を示すフローチャートである。
【符号の説明】
【００９２】
１  データ予測装置（パラメータ推定装置）
１１  ＣＰＵ（演算部）
１２  ＲＡＭ（記憶部）
２  記録媒体
２０  コンピュータプログラム
【図１】
2005242580.tif 000040
【図２】
2005242580.tif 000041