JP2005242580A - パラメータ推定方法、データ予測方法、パラメータ推定装置、データ予測装置、及びコンピュータプログラム - Google Patents

Abstract

【課題】 入力を伴う系から得られるデータからなる時系列のパラメータを確実に推定する方法、現実に即したデータの予測を行う方法、パラメータ推定装置、データ予測装置及びコンピュータプログラムを提供する。
【解決手段】 データ予測装置（パラメータ推定装置）は、現在の標本時系列｛ｘ_n ｝と未来の入力時系列｛ｕ_n ｝との相互相関関数を０としてデータを計算し（Ｓ２）、白色雑音｛ｅ_n ｝を生成し（Ｓ３）、｛ｅ_n ｝の分散σ_e ²を計算し（Ｓ４）、｛ｕ_n ｝と｛ｅ_n ｝との相互相関関数、現在の｛ｘ_n ｝と未来の｛ｅ_n ｝との相互相関関数を０としてデータを計算し（Ｓ５）、｛ｅ_n ｝の自己相関関数をラグ０でσ_e ²、他のラグで０として、ＡＲＭＡＸモデルのパラメータを計算する（Ｓ６）。パラメータの収束を判定し（Ｓ７）、収束していない場合は新たな白色雑音を計算し（Ｓ８）、パラメータが収束するまで計算を繰り返す。
【選択図】 図２

Description

本発明は、データの時系列からパラメータを推定する方法、推定されたパラメータを用いて新たなデータを予測する方法、パラメータ推定装置、データ予測装置、コンピュータをパラメータ推定装置又はデータ予測装置として実現するためのコンピュータプログラムに関する。

あるシステムから得られるデータの時系列を解析して新たなデータの値を予想する方法は、従来様々な方法が提案されている。その方法の一つとして、データの時系列を離散時間線形モデルに当てはめて解析を行い、データ間の関係を示すパラメータを推定し、推定したパラメータを用いて構成される離散時間線形モデルに基づいて新たなデータの値を予測する方法がある。離散時間線形モデルとしては、ＡＲ（Autoregressive：自己回帰）モデル又はＡＲＭＡ（AutoregressiveMoving Average ：自己回帰移動平均）モデル等が利用されている。

データの時系列｛ｘ_n ：ｎ＝１，２，…｝を、定常エルゴード的正規過程ｘ（ｔ）の標本時系列であるとすると、（ｐ，ｑ）次のＡＲＭＡモデルは下記（１）式で表される。

ここで、｛ｅ_n ｝は、平均値０及び分散σ_e ²の正規性白色雑音である。また、｛ａ_i ：i＝１，２，…，ｐ｝及び｛ｂ_j ：ｊ＝１，２，…，ｑ｝はＡＲＭＡモデルのパラメータである。また、｛ａ_i ｝及び｛ｂ_j ｝は、下記式で与えられるＡ（Ｚ^-1）及びＢ（Ｚ^-1）が規約で夫々定常条件、可逆条件および強正実条件を満たすものとする。
Ａ（Ｚ^-1）＝１＋ａ₁ Ｚ^-1＋…＋ａ_p Ｚ^-p
Ｂ（Ｚ^-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：ｎ＝１，２，…｝を定常エルゴード的正規過程ｘ（ｔ）の標本時系列であるとすると、（ｐ，ｍ，ｑ）次のＡＲＭＡＸモデルは、下記（２）式で表される。

ここで、｛ｅ_n ｝は平均値０及び分散σ_e ²の正規性白色雑音であり、｛ｕ_n ｝はシステムに対する入力時系列である。また｛ａ_i ：i＝１，２，…，ｐ｝、｛ｃ_i ：ｉ＝１，２，…，ｍ｝及び｛ｂ_i：ｉ＝１，２，…，ｑ｝はＡＲＭＡＸモデルのパラメータである。また、（２）式は定常条件を満たすとする。現代制御理論では、ｑ≦ｐ，ｍ≦ｐとなることは周知である。以下、簡単にするためにｑ＝ｍ＝ｐとして議論を進める。
ベクトルＺ_n 及びベクトルθを下記のように定義する。
Ｚ_n ＝［-ｘ_n-1,…,-ｘ_n-p,ｕ_n-1,…,ｕ_n-p,ｅ_n-1,…,ｅ_n-p］^T
θ＝［ａ₁,…,ａ_p,ｃ₁,…,ｃ_p,ｂ₁,…,ｂ_p］^T
（２）式をｘ_n についての式に書き直すと、（２）式は下記（３）式で表すことができる。
ｘ_n ＝θ^T Ｚ_n＋ｅ_n …（３）

Ｎ個のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…，Ｎ｝及び入力時系列｛ｕ_n：ｎ＝１，２，…，Ｎ｝が得られたとき、ベクトルθを時系列から最小二乗法により求めるためには、（３）式から得られる下記（４）式の右辺を最小にするθを求めればよい。

従って、（４）式の右辺をθで偏微分して０に等しいとすることにより、Ｎ個のデータからなる標本時系列｛ｘ_n ｝及び入力時系列｛ｕ_n ｝から推定したベクトルθの推定値θ_N は、下記（５）式で与えられる。

ｋ≧０としてラグがｋである｛ｘ_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

ベクトルａ，ｂ，ｃを下記式で定義する。
ａ＝［ａ₁ ，ａ₂ ，…，ａ_p］^T
ｂ＝［ｂ₁ ，ｂ₂ ，…，ｂ_p］^T
ｃ＝［ｃ₁ ，ｃ₂ ，…，ｃ_p］^T

（５）式においてＮ→∞とすると、（５）式の左辺はθ_N →θとなり、（５）式の右辺は｛ｘ_n ｝，｛ｕ_n ｝，｛ｅ_n ｝の自己相関関数及び相互相関関数で記述できる。従って、（５）式は下記（６）式で表すことができる。

但し、Ｒ，Ｖ’，Ｔ’，Ｗ，Ｓ’，Ｑ’はｐ行ｐ列の行列であって、夫々下記（７），，（８），（９），（１０），（１１），（１２）式で表される。また、ｒ，ｖ，ｔはｐ行一列の行列であって、下記（１３）式で表される。

（６）式を書き直すと、下記の連立方程式が得られる。
Ｒａ−Ｖ’ｃ−Ｔ’ｂ＝−ｒ …（１４）
−Ｖ’^T ａ＋Ｗｃ＋Ｑ’ｂ＝ｖ …（１５）
−Ｔ’^T ａ＋Ｑ’^T ｃ＋Ｓ’ｂ＝ｔ …（１６）
この連立方程式（１４）〜（１６）を解くことにより、ベクトルａ，ｂ，ｃを計算することができる。

疑似線形回帰法を用いてＡＲＭＡＸモデルのパラメータの推定値であるベクトルａ，ｂ，ｃを求めるためには、白色雑音｛ｅ_n ｝を乱数などを用いて生成し、生成した｛ｅ_n｝と｛ｘ_n ｝及び｛ｕ_n｝とから、ベクトルａ，ｂ，ｃを（１４）〜（１６）式により計算し、計算したベクトルａ，ｂ，ｃから（２）式を用いて新たな白色雑音｛ｅ_n ｝を計算し、ベクトルａ，ｂ，ｃが収束するまで計算を繰り返す。
竹安数博、外３名，「ロバストなブートストラップ型システム同定アルゴリズム」，Journal of Signal Processing，信号処理学会，２００３年３月，第７巻，第２号，ｐ．１６７−１７６

前述の如き疑似線形回帰法は、問題の性質によっては、求めるべきパラメータであるベクトルａ，ｂ，ｃが収束せず、パラメータの推定ができない場合があるという問題がある。また、前述の方法では、ベクトルａ，ｂ，ｃを計算するためにｐ行ｐ列の行列Ｒ，Ｖ’，Ｔ’，Ｗ，Ｓ’，Ｑ’を用いた行列計算を行う必要があるため、計算時間が長いという問題があった。

本発明は、斯かる事情に鑑みてなされたものであって、その目的とするところは、外部入力を考慮に入れたＡＲＭＡＸモデルのパラメータの値の収束を容易にし、同時に計算時間の短縮を図ることができるパラメータ推定方法、データ予測方法、パラメータ推定装置、データ予測装置、及びコンピュータプログラムを提供することにある。

第１発明に係るパラメータ推定方法は、記憶部及び演算部を備えたコンピュータを用いて、複数のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…｝及び入力時系列｛ｕ_n：ｎ＝１，２，…｝を（ｐ，ｐ，ｐ）次（但し、ｐは自然数）のＡＲＭＡＸモデルに当てはめたときのパラメータを推定する方法において、前記標本時系列及び前記入力時系列を前記記憶部で記憶し、下記式で定義される行列Ｒ，Ｗ，Ｖ，ｒ，ｖを前記記憶部で記憶し、

平均値が０である複数のランダムなデータからなる白色雑音｛ｅ_n ：ｎ＝１，２，…｝を前記演算部で生成し、前記演算部で生成した白色雑音の分散σ_e ²を前記演算部で計算し、下記式で定義される行列Ｔ及びｔを前記演算部で計算し、

前記ＡＲＭＡＸモデルにおける前記標本時系列に係るｐ個のパラメータを行列要素とするベクトルａ＝［ａ₁ ，ａ₂ ，…，ａ_p ］^T、前記入力時系列に係るｐ個のパラメータを行列要素とする行列ｃ＝［ｃ₁ ，ｃ₂ ，…，ｃ_p ］^T及び前記白色雑音に係るｐ個のパラメータを行列要素とする行列ｂ＝［ｂ₁ ，ｂ₂ ，…，ｂ_p ］^Tを、下記式に基づいて前記演算部で計算し、

前記演算部で計算したベクトルａ、ｂ及びｃに含まれる前記ＡＲＭＡＸモデルのパラメータ、並びに前記標本時系列及び前記入力時系列値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を前記演算部で計算し、前記演算部で計算した新たな白色雑音を用いて新たなベクトルａ、ｂ及びｃを前記演算部で計算し、ベクトルａ、ｂ及びｃが収束するまでベクトルａ、ｂ及びｃの計算を前記演算部で繰り返すことを特徴とする。

第２発明に係るデータ予測方法は、記憶部及び演算部を備えたコンピュータを用いて、入力を伴う系から得られる複数のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…｝、及び前記系へ入力される複数のデータからなる入力時系列｛ｕ_n：ｎ＝１，２，…｝から、前記標本時系列に加わるべき新たなデータを予測する方法において、前記標本時系列及び前記入力時系列を前記記憶部で記憶し、第１発明に係るパラメータ推定方法を用いて、前記標本時系列及び前記入力時系列をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を前記演算部で計算し、前記演算部で計算した前記パラメータの推定値、並びに前記標本時系列及び前記入力時系列を前記ＡＲＭＡモデルの式に代入して、前記標本時系列に加わるべき新たなデータの予測値を前記演算部で計算することを特徴とする。

第３発明に係るパラメータ推定装置は、複数のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…｝及び入力時系列｛ｕ_n：ｎ＝１，２，…｝を（ｐ，ｐ，ｐ）次（但し、ｐは自然数）のＡＲＭＡＸモデルに当てはめたときのパラメータを推定する装置において、前記標本時系列及び前記入力時系列を記憶する手段と、下記式で定義される行列Ｒ，Ｗ，Ｖ，ｒ，ｖを記憶する手段と、

平均値が０である複数のランダムなデータからなる白色雑音｛ｅ_n ：ｎ＝１，２，…｝を生成する手段と、該手段が生成した白色雑音の分散σ_e ²を計算する手段と、下記式で定義される行列Ｔ及びｔを計算する手段と、

前記ＡＲＭＡＸモデルにおける前記標本時系列に係るｐ個のパラメータの推定値を行列要素とするベクトルａ＝［ａ₁ ，ａ₂ ，…，ａ_p ］^T、前記入力時系列に係るｐ個のパラメータの推定値を行列要素とする行列ｃ＝［ｃ₁ ，ｃ₂ ，…，ｃ_p］^T 及び前記白色雑音に係るｐ個のパラメータの推定値を行列要素とする行列ｂ＝［ｂ₁ ，ｂ₂ ，…，ｂ_p］^T を、下記式に基づいて計算する手段と、

該手段が計算したベクトルａ、ｂ及びｃに含まれる前記ＡＲＭＡＸモデルのパラメータの推定値、並びに前記標本時系列及び前記入力時系列値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算する手段と、該手段が計算した新たな白色雑音を用いて新たなベクトルａ、ｂ及びｃを計算する手段と、ベクトルａ、ｂ及びｃが収束するまでベクトルａ、ｂ及びｃの計算を繰り返す手段とを備えることを特徴とする。

第４発明に係るデータ予測装置は、入力を伴う系から得られる複数のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…｝、及び前記系へ入力される複数のデータからなる入力時系列｛ｕ_n：ｎ＝１，２，…｝から、前記標本時系列に加わるべき新たなデータを予測する装置において、前記標本時系列及び前記入力時系列を記憶する手段と、第３発明に係るパラメータ推定装置と、該パラメータ推定装置を用いて、前記標本時系列及び前記入力時系列をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を計算する手段と、該手段が計算した前記パラメータの推定値、並びに前記標本時系列及び前記入力時系列を前記ＡＲＭＡモデルの式に代入して、前記標本時系列に加わるべき新たなデータの予測値を計算する手段とを備えることを特徴とする。

第５発明に係るコンピュータプログラムは、複数のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…｝及び入力時系列｛ｕ_n：ｎ＝１，２，…｝、並びに下記式で定義される行列Ｒ，Ｗ，Ｖ，ｒ，ｖを記憶するコンピュータに、

前記標本時系列及び前記入力時系列を（ｐ，ｐ，ｐ）次（但し、ｐは自然数）のＡＲＭＡＸモデルに当てはめたときのパラメータを推定させるコンピュータプログラムであって、コンピュータに、平均値が０である複数のランダムなデータからなる白色雑音｛ｅ_n ：ｎ＝１，２，…｝を生成させる手順と、コンピュータに、生成した白色雑音の分散σ_e ²を計算させる手順と、コンピュータに、下記式で定義される行列Ｔ及びｔを計算させる手順と、

コンピュータに、前記ＡＲＭＡＸモデルにおける前記標本時系列に係るｐ個のパラメータの推定値を行列要素とするベクトルａ＝［ａ₁ ，ａ₂ ，…，ａ_p ］^T、前記入力時系列に係るｐ個のパラメータの推定値を行列要素とする行列ｃ＝［ｃ₁ ，ｃ₂ ，…，ｃ_p］^T 及び前記白色雑音に係るｐ個のパラメータの推定値を行列要素とする行列ｂ＝［ｂ₁ ，ｂ₂ ，…，ｂ_p］^T を、下記式に基づいて計算させる手順と、

コンピュータに、計算したベクトルａ、ｂ及びｃに含まれる前記ＡＲＭＡＸモデルのパラメータの推定値、並びに前記標本時系列及び前記入力時系列値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算させる手順と、コンピュータに、計算した新たな白色雑音を用いて新たなベクトルａ、ｂ及びｃを計算させる手順と、コンピュータに、ベクトルａ、ｂ及びｃが収束するまでベクトルａ、ｂ及びｃの計算を繰り返させる手順とを含むことを特徴とする。

第６発明に係るコンピュータプログラムは、入力を伴う系から得られる複数のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…｝、及び前記系へ入力される複数のデータからなる入力時系列｛ｕ_n：ｎ＝１，２，…｝を記憶するコンピュータに、前記標本時系列に加わるべき新たなデータを予測させるコンピュータプログラムであって、コンピュータに、第１発明に係るパラメータ推定方法を用いて、前記標本時系列及び前記入力時系列をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を計算させる手順と、コンピュータに、計算した前記パラメータの推定値、並びに前記標本時系列及び前記入力時系列を前記ＡＲＭＡモデルの式に代入して、前記標本時系列に加わるべき新たなデータの予測値を計算させる手順とを含むことを特徴とする。

本発明者は、ＡＲＭＡＸモデルの式の構成から導かれるa priori knowledgeを導入することにより、従来の疑似線形回帰法の改良を試みた。まず正規性雑音｛ｅ_n ｝の自己相関関数を考える。正規性白色雑音｛ｅ_n｝に含まれる各要素は互いに関連性がないので、｛ｅ_n ｝の自己相関関数Ｓ_ｋは理論的に下記式で表すことができる。

従って、（１１）式で定義される行列Ｓ’に理論的な値を導入することによって行列Ｓ’を改良した行列Ｓは、単位行列をＩとして下記式で表される。

次に、｛ｘ_n ｝と｛ｅ_n｝との相互相関関数を考える。（２）式で示されるＡＲＭＡＸモデルでは、現在のｘ_n の値は過去のｅ_n の値に関連づけられている。しかし、現在のｘ_nの値は未来のｅ_n の値とは関連性がないので、ｘ_n とこのデータよりも未来のｅ_n との相互相関関数は下記式で表すことができる。
Ｅ［ｘ_n ｅ_n+k ］＝０
同様に下記式が成り立つ。
Ｅ［ｅ_n ｘ_n-k ］＝０
従って、Ｔ_-k＝０となる。また（２）式より現在のｘ_n の値と現在のｅ_n の値とは関連性があるので、Ｔ₀≠０である。従って、（９）式で定義される行列Ｔ’に理論的な値を導入することによって行列Ｔ’を改良した行列Ｔは下記式で表される。

｛ｘ_n ｝と｛ｕ_n ｝との相互相関関数についても、同様に、現在のｘ_nの値は過去のｕ_n の値に関連づけられているが、現在のｘ_n の値は未来のｕ_n の値とは関連性がないので、下記式が成り立つ。
Ｅ［ｘ_n ｕ_n+k ］＝Ｅ［ｕ_n ｘ_n-k ］＝Ｖ_-k＝０
また（２）式より現在のｘ_n の値と現在のｕ_n の値とは関連性がないので、Ｖ₀ ＝０である。従って、（８）式で定義される行列Ｖ’に理論的な値を導入することによって行列Ｖ’を改良した行列Ｖは下記式で表される。

更に、入力時系列｛ｕ_n ｝と正規性白色雑音｛ｅ_n ｝とは無相関であるはずであるので、下記式に示す如く｛ｕ_n ｝と｛ｅ_n｝との相互相関関数は全て０となる。
Ｅ［ｅ_n ｕ_n+k ］＝Ｅ［ｕ_n ｅ_n-k ］＝Ｑ_k＝０
Ｅ［ｅ_n ｕ_n-k ］＝Ｅ［ｕ_nｅ_n+k ］＝Ｑ_-k＝０
従って、（１２）式で定義される行列Ｑ’に理論的な値を導入することによって行列Ｑ’を改良した行列Ｑは、ゼロ行列をＯとして下記式で表される。
Ｑ＝Ｏ …（２１）

Ｎ個のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…，Ｎ｝及び入力時系列｛ｕ_n：ｎ＝１，２，…，Ｎ｝が得られたときに各自己相関関数及び各相互相関関数を計算した場合は、Ｎの値が有限であるので相互相関関数及び自己相関関数の値は理論的な値とは違った値となる。この実際に計算される値の内、理論的な値が得られる部分は理論的な値と置き換えることにより、ＡＲＭＡＸモデルのパラメータを求めるための（６）式は下記式で表すことができるようになる。

（２２）式を書き直すと、下記の連立方程式が得られる。
Ｒａ−Ｖｃ−Ｔｂ＝−ｒ …（２３）
−Ｖ^T ａ＋Ｗｃ＝ｖ …（２４）
−Ｔ^T ａ＋σ_e ²ｂ＝ｔ …（２５）
（２５）式及び（２４）より、下記の（２６）式及び（２７）式が得られる。
ｂ＝（ｔ＋Ｔ^T ａ）／σ_e ² …（２６）
ｃ＝Ｗ^-1（ｖ＋Ｖ^T ａ） …（２７）
（２６）式及び（２７）式を（２３）式に代入することにより、下記（２８）式が得られる。

（２８）式よりベクトルａを計算することができる。計算したベクトルａを（２６）式及び（２７）式へ代入することにより、ベクトルｂ及びベクトルｃを計算することができる。ＡＲＭＡＸモデルのパラメータの推定値であるベクトルａ，ｂ，ｃを求めるためには、白色雑音｛ｅ_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 ｝の分散σ_e ²を計算する（Ｓ４）。

ＣＰＵ１１は、次に、｛ｘ_n ｝及び｛ｅ_n ｝に含まれるデータの値を用いて｛ｘ_n ｝と｛ｅ_n ｝との相互相関関数Ｔ_k＝Ｅ［ｅ_n ｘ_n+k ］を計算して、（１９）式，（１３）式で定義される行列Ｔ，ｔを計算する（Ｓ５）。ＣＰＵ１１は、次に、（２８）式によりベクトルａを計算し、ベクトルａを（２６）式及び（２７）式へ代入することによりベクトルｂ及びベクトルｃを計算することによって、｛ｘ_n ｝及び｛ｕ_n ｝をＡＲＭＡＸモデルに当てはめたときのＡＲＭＡＸモデルのパラメータの推定値からなるベクトルａ，ｂ，ｃを計算する（Ｓ６）。

ＣＰＵ１１は、次に、計算したベクトルａ，ｂ，ｃが収束しているか否かを判定する（Ｓ７）。このときＣＰＵ１１は、例えば、最新のベクトルａ，ｂ，ｃと前回に計算したベクトルａ，ｂ，ｃであるベクトルａ’＝［ａ’₁ ，ａ’₂ ，…，ａ’_p ］^T，ｂ’＝［ｂ’₁ ，ｂ’₂ ，…，ｂ’_p ］^T ，ｃ’＝［ｃ’₁，ｃ’₂ ，…，ｃ’_p ］^T を用いて下記（２９）式で定義されるＤを計算し、計算したＤが所定の微少量εに対してＤ＜εである状態が所定回数連続で成立している場合にベクトルａ，ｂ，ｃが収束していると判定する。なお、その他の方法を用いてベクトルａ，ｂ，ｃの収束を判定する処理を用いてもよい。

ステップＳ７にて、ベクトルａ，ｂ，ｃが収束していないと判定された場合は（Ｓ７：ＮＯ）、ベクトルａ，ｂ，ｃと｛ｘ_n ｝及び｛ｕ_n ｝とをｍ＝ｑ＝ｐとした（２）式へ代入して、新たな白色雑音｛ｅ_n：ｎ＝１，２，…｝を計算し（Ｓ８）、処理をステップＳ４へ戻す。ＣＰＵ１１は、ステップＳ７にてベクトルａ，ｂ，ｃが収束したと判定されるまでステップＳ４〜Ｓ８の処理を繰り返す。

ステップＳ７にて、ベクトルａ，ｂ，ｃが収束していると判定された場合は（Ｓ７：ＹＥＳ）、ＣＰＵ１１は、計算したベクトルａ，ｂ，ｃ、｛ｘ_n ｝、｛ｕ_n ｝、｛ｅ_n ｝を（２）式へ代入することにより、標本時系列｛ｘ_n｝に加わるべきデータの予測値、即ち株価又は商品の販売量などの予測データを計算する（Ｓ９）。なお、ステップＳ９では、入力時系列に加わる新たなデータ又は当該データの予測値を必要に応じて用いてもよい。以上の処理により、ＡＲＭＡＸモデルのパラメータが推定され、標本時系列に加わるべき新たなデータが予測できる。

以上詳述した如く、本発明においては、入力を伴う系から得られた複数のデータからなる標本時系列｛ｘ_n ｝及び系へ入力された複数のデータからなる入力時系列｛ｕ_n ｝をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を計算する際に、白色雑音の自己相関関数、標本時系列と入力時系列との相互相関関数、標本時系列と白色雑音の相互相関関数、及び入力時系列と白色雑音との相互相関関数に理論的な値を用いて計算する。これにより、従来の疑似線形回帰法でパラメータの推定値の収束が困難であった場合についても、パラメータの推定値の収束が容易となる。またパラメータが収束する確率が向上するので、パラメータ推定の計算時間を短縮する事ができる。

また、本発明においては、系から得られるデータの値が過去のデータの値と雑音と更に系への入力とから決定されるＡＲＭＡＸモデルを離散時間線形モデルとして用いることにより、株の取引量に影響される株価の値、又は商品の入荷量に影響される商品の販売量など、外部入力に影響される系から得られるデータの値を推定することができる。従って、株価又は商品の販売量などの予測をより現実に即して行うことができる。

なお、本実施の形態においては、本発明のデータ予測装置は、内部で標本時系列｛ｘ_n ｝及び入力時系列｛ｕ_n ｝を記憶している形態を示しているが、これに限るものではなく、標本時系列｛ｘ_n ｝及び入力時系列｛ｕ_n ｝に含まれるデータを外部から受け付ける手段を備え、受け付けたデータについて予測データを計算する処理を行う形態であってもよい。

また、本実施の形態においては、本発明のデータ予測装置１は、本発明のパラメータ推定装置の機能を兼ね備えている形態を示したが、これに限るものではなく、ＡＲＭＡＸモデルのパラメータの予測値を計算するパラメータ推定装置とＡＲＭＡＸモデルによりデータを予測するデータ予測装置とを夫々別のコンピュータにて構成し、互いに情報を入出力してパラメータの推定とデータの予測とを行う形態としても良い。

本発明のデータ予測装置の構成を示すブロック図である。 本発明のデータ予測装置が行うパラメータ推定の処理およびデータ予測の処理を示すフローチャートである。

符号の説明

１ データ予測装置（パラメータ推定装置）
１１ ＣＰＵ（演算部）
１２ ＲＡＭ（記憶部）
２ 記録媒体
２０ コンピュータプログラム

Claims (6)

1. 記憶部及び演算部を備えたコンピュータを用いて、複数のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…｝及び入力時系列｛ｕ_n ：ｎ＝１，２，…｝を（ｐ，ｐ，ｐ）次（但し、ｐは自然数）のＡＲＭＡＸモデルに当てはめたときのパラメータを推定する方法において、
   前記標本時系列及び前記入力時系列を前記記憶部で記憶し、
   下記式で定義される行列Ｒ，Ｗ，Ｖ，ｒ，ｖを前記記憶部で記憶し、
   

   

   平均値が０である複数のランダムなデータからなる白色雑音｛ｅ_n ：ｎ＝１，２，…｝を前記演算部で生成し、
   前記演算部で生成した白色雑音の分散σ_e ²を前記演算部で計算し、
   下記式で定義される行列Ｔ及びｔを前記演算部で計算し、
   

   

   前記ＡＲＭＡＸモデルにおける前記標本時系列に係るｐ個のパラメータを行列要素とするベクトルａ＝［ａ₁ ，ａ₂ ，…，ａ_p ］^T 、前記入力時系列に係るｐ個のパラメータを行列要素とする行列ｃ＝［ｃ₁，ｃ₂ ，…，ｃ_p ］^T 及び前記白色雑音に係るｐ個のパラメータを行列要素とする行列ｂ＝［ｂ₁，ｂ₂ ，…，ｂ_p ］^T を、下記式に基づいて前記演算部で計算し、
   

   

   前記演算部で計算したベクトルａ、ｂ及びｃに含まれる前記ＡＲＭＡＸモデルのパラメータ、並びに前記標本時系列及び前記入力時系列値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を前記演算部で計算し、
   前記演算部で計算した新たな白色雑音を用いて新たなベクトルａ、ｂ及びｃを前記演算部で計算し、
   ベクトルａ、ｂ及びｃが収束するまでベクトルａ、ｂ及びｃの計算を前記演算部で繰り返すこと
   を特徴とするパラメータ推定方法。
2. 記憶部及び演算部を備えたコンピュータを用いて、入力を伴う系から得られる複数のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…｝、及び前記系へ入力される複数のデータからなる入力時系列｛ｕ_n：ｎ＝１，２，…｝から、前記標本時系列に加わるべき新たなデータを予測する方法において、
   前記標本時系列及び前記入力時系列を前記記憶部で記憶し、
   請求項１に記載のパラメータ推定方法を用いて、前記標本時系列及び前記入力時系列をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を前記演算部で計算し、
   前記演算部で計算した前記パラメータの推定値、並びに前記標本時系列及び前記入力時系列を前記ＡＲＭＡモデルの式に代入して、前記標本時系列に加わるべき新たなデータの予測値を前記演算部で計算すること
   を特徴とするデータ予測方法。
3. 複数のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…｝及び入力時系列｛ｕ_n ：ｎ＝１，２，…｝を（ｐ，ｐ，ｐ）次（但し、ｐは自然数）のＡＲＭＡＸモデルに当てはめたときのパラメータを推定する装置において、
   前記標本時系列及び前記入力時系列を記憶する手段と、
   下記式で定義される行列Ｒ，Ｗ，Ｖ，ｒ，ｖを記憶する手段と、
   

   

   平均値が０である複数のランダムなデータからなる白色雑音｛ｅ_n ：ｎ＝１，２，…｝を生成する手段と、
   該手段が生成した白色雑音の分散σ_e ²を計算する手段と、
   下記式で定義される行列Ｔ及びｔを計算する手段と、
   

   

   前記ＡＲＭＡＸモデルにおける前記標本時系列に係るｐ個のパラメータの推定値を行列要素とするベクトルａ＝［ａ₁ ，ａ₂ ，…，ａ_p ］^T、前記入力時系列に係るｐ個のパラメータの推定値を行列要素とする行列ｃ＝［ｃ₁ ，ｃ₂ ，…，ｃ_p］^T 及び前記白色雑音に係るｐ個のパラメータの推定値を行列要素とする行列ｂ＝［ｂ₁ ，ｂ₂ ，…，ｂ_p］^T を、下記式に基づいて計算する手段と、
   

   

   該手段が計算したベクトルａ、ｂ及びｃに含まれる前記ＡＲＭＡＸモデルのパラメータの推定値、並びに前記標本時系列及び前記入力時系列値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算する手段と、
   該手段が計算した新たな白色雑音を用いて新たなベクトルａ、ｂ及びｃを計算する手段と、
   ベクトルａ、ｂ及びｃが収束するまでベクトルａ、ｂ及びｃの計算を繰り返す手段と
   を備えることを特徴とするパラメータ推定装置。
4. 入力を伴う系から得られる複数のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…｝、及び前記系へ入力される複数のデータからなる入力時系列｛ｕ_n：ｎ＝１，２，…｝から、前記標本時系列に加わるべき新たなデータを予測する装置において、
   前記標本時系列及び前記入力時系列を記憶する手段と、
   請求項３に記載のパラメータ推定装置と、
   該パラメータ推定装置を用いて、前記標本時系列及び前記入力時系列をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を計算する手段と、
   該手段が計算した前記パラメータの推定値、並びに前記標本時系列及び前記入力時系列を前記ＡＲＭＡモデルの式に代入して、前記標本時系列に加わるべき新たなデータの予測値を計算する手段と
   を備えることを特徴とするデータ予測装置。
5. 複数のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…｝及び入力時系列｛ｕ_n ：ｎ＝１，２，…｝、並びに下記式で定義される行列Ｒ，Ｗ，Ｖ，ｒ，ｖを記憶するコンピュータに、
   

   

   前記標本時系列及び前記入力時系列を（ｐ，ｐ，ｐ）次（但し、ｐは自然数）のＡＲＭＡＸモデルに当てはめたときのパラメータを推定させるコンピュータプログラムであって、
   コンピュータに、平均値が０である複数のランダムなデータからなる白色雑音｛ｅ_n ：ｎ＝１，２，…｝を生成させる手順と、
   コンピュータに、生成した白色雑音の分散σ_e ²を計算させる手順と、
   コンピュータに、下記式で定義される行列Ｔ及びｔを計算させる手順と、
   

   

   コンピュータに、前記ＡＲＭＡＸモデルにおける前記標本時系列に係るｐ個のパラメータの推定値を行列要素とするベクトルａ＝［ａ₁ ，ａ₂ ，…，ａ_p ］^T、前記入力時系列に係るｐ個のパラメータの推定値を行列要素とする行列ｃ＝［ｃ₁ ，ｃ₂ ，…，ｃ_p］^T 及び前記白色雑音に係るｐ個のパラメータの推定値を行列要素とする行列ｂ＝［ｂ₁ ，ｂ₂ ，…，ｂ_p］^T を、下記式に基づいて計算させる手順と、
   

   

   コンピュータに、計算したベクトルａ、ｂ及びｃに含まれる前記ＡＲＭＡＸモデルのパラメータの推定値、並びに前記標本時系列及び前記入力時系列値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算させる手順と、
   コンピュータに、計算した新たな白色雑音を用いて新たなベクトルａ、ｂ及びｃを計算させる手順と、
   コンピュータに、ベクトルａ、ｂ及びｃが収束するまでベクトルａ、ｂ及びｃの計算を繰り返させる手順と
   を含むことを特徴とするコンピュータプログラム。
6. 入力を伴う系から得られる複数のデータからなる標本時系列｛ｘ_n ：ｎ＝１，２，…｝、及び前記系へ入力される複数のデータからなる入力時系列｛ｕ_n：ｎ＝１，２，…｝を記憶するコンピュータに、前記標本時系列に加わるべき新たなデータを予測させるコンピュータプログラムであって、
   コンピュータに、請求項１に記載のパラメータ推定方法を用いて、前記標本時系列及び前記入力時系列をＡＲＭＡＸモデルに当てはめたときのパラメータの推定値を計算させる手順と、
   コンピュータに、計算した前記パラメータの推定値、並びに前記標本時系列及び前記入力時系列を前記ＡＲＭＡモデルの式に代入して、前記標本時系列に加わるべき新たなデータの予測値を計算させる手順と
   を含むことを特徴とするコンピュータプログラム。

Images