(57)【要約】【課題】二つの異なる時系列データ間で因果関係が存在しているのか否かを判定できるようにする。【解決手段】ＲＡＭに記憶された二つの異なる時系列データＸｄａｔａ，Ｙｄａｔａに対して、ＹｄａｔａからＸｄａｔａへの各時刻ｔの因果関係の程度とその因果関係に対するＹ→Ｘ重みｗy→x（ｔ）とを用いた評価式の値が最大となるよう各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）を演算すると共にこれを用いてＹｄａｔａからＸｄａｔａへの各時刻ｔの因果関係の有無を判定して結果を出力する（Ｓ１１０〜Ｓ１３０）。そして、同様の手法により、各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）を演算すると共にこれをを用いてＸｄａｔａからＹｄａｔａへの各時刻ｔの因果関係の有無を判定して結果を出力する（Ｓ１４０〜Ｓ１６０）。【選択図】図２
2008257595.tif 000002
【特許請求の範囲】
【請求項１】
二つの異なる時系列データ間の各時刻の因果関係の有無を判定する因果関係判定装置であって、
情報を入力する入力手段と、
情報を記憶する記憶手段と、
情報を出力する出力手段と、
前記二つの異なる時系列データを前記入力手段により入力すると共に該入力した二つの異なる時系列データを前記記憶手段に記憶する入力処理手段と、
前記記憶手段に記憶された二つの異なる時系列データに基づく該二つの異なる時系列データ間の各時刻の前記因果関係の程度と該因果関係の程度に対する各時刻の重みとを用いて所定の演算手法により演算される結果が最大となるよう各時刻の前記重みを演算し、該演算結果を前記記憶手段に記憶する重み演算処理手段と、
前記記憶手段に記憶された各時刻の重みを用いて該記憶手段に記憶された二つの異なる時系列データ間の各時刻の前記因果関係の有無を判定し、該判定結果を前記出力手段により出力する判定出力処理手段と、
を備える因果関係判定装置。
【請求項２】
請求項１記載の因果関係判定装置であって、
前記二つの異なる時系列データは、Ｘｄａｔａ，Ｙｄａｔａと置くと、式（１）および式（２）により表わされるものであり、
前記重み演算処理手段は、任意の時刻における前記Ｙｄａｔａから前記Ｘｄａｔａへの各時刻の前記因果関係の程度であるＹ→Ｘ程度と該Ｙｄａｔａから該Ｘｄａｔａへの各時刻の前記重みであるＹ→Ｘ重みとを用いた評価式であるＹ→Ｘ評価式の値が最大となるよう該Ｙｄａｔａから該Ｘｄａｔａへの各時刻の該重みを演算すると共に該演算結果を前記記憶手段に記憶し、任意の時刻における前記Ｘｄａｔａから前記Ｙｄａｔａへの各時刻の前記因果関係の程度であるＸ→Ｙ程度と該Ｘｄａｔａから該Ｙｄａｔａへの各時刻の前記重みであるＸ→Ｙ重みとを用いた評価式であるＸ→Ｙ評価式の値が最大となるよう該Ｘｄａｔａから該Ｙｄａｔａへの各時刻の該重みを演算すると共に該演算結果を前記記憶手段に記憶する手段である
因果関係判定装置。
【数１】
2008257595.tif 000025
ここで、式（１）中、ｘ（１），…，ｘ（Ｎ）は、時刻１，…，ＮのＸｄａｔａの値である。また、式（２）中、ｙ（１），…，ｙ（Ｎ）は、時刻１，…，ＮのＹｄａｔａの値である。
【請求項３】
請求項２記載の因果関係判定装置であって、
前記Ｙ→Ｘ程度は、各時刻に対して、前記Ｘｄａｔａの第１の時刻の入力値と該Ｘｄａｔａの該第１の時刻よりも前の時刻の１以上の入力値を用いて予測される該Ｘｄａｔａの該第１の時刻の予測値との誤差であるＸ→Ｘ誤差と、前記Ｘｄａｔａの前記第１の時刻の入力値と該Ｘｄａｔａおよび前記Ｙｄａｔａの該第１の時刻よりも前の時刻の１以上の入力値を用いて予測される該Ｘｄａｔａの該第１の時刻の予測値との誤差であるＸＹ→Ｘ誤差と、を用いて表わされるものであり、
前記Ｘ→Ｙ程度は、各時刻に対して、前記Ｙｄａｔａの前記第１の時刻の入力値と該Ｙｄａｔａの該第１の時刻よりも前の時刻の１以上の入力値を用いて予測される該Ｙｄａｔａの該第１の時刻の予測値との誤差であるＹ→Ｙ誤差と、前記Ｙｄａｔａの前記第１の時刻の入力値と前記Ｘｄａｔａおよび該Ｙｄａｔａの該第１の時刻よりも前の時刻の１以上の入力値を用いて予測される該Ｙｄａｔａの該第１の時刻の予測値との誤差であるＹＸ→Ｙ誤差と、を用いて表わされるものである
因果関係判定装置。
【請求項４】
請求項３記載の因果関係判定装置であって、
前記Ｙ→Ｘ評価式は、前記Ｙ→Ｘ重みと、前記Ｘ→Ｘ誤差の２乗と前記ＸＹ→Ｘ誤差の２乗との差と、の積の期待値であるＹ→Ｘ期待値または該Ｙ→Ｘ期待値に第１の補正項を加えた値を演算する式であり、
前記Ｘ→Ｙ評価式は、前記Ｘ→Ｙ重みと、前記Ｙ→Ｙ誤差の２乗と前記ＹＸ→Ｙ誤差の２乗との差と、の積の期待値であるＸ→Ｙ期待値または該Ｘ→Ｙ期待値に第２の補正項を加えた値を演算する式である
因果関係判定装置。
【請求項５】
請求項３または４記載の因果関係判定装置であり、
前記Ｘ→Ｘ誤差は、時刻ｔの誤差を「ε(x)（ｔ）」と置くと、式（３）により表わされ、
前記ＸＹ→Ｘ誤差は、時刻ｔの誤差を「ε(x,y)（ｔ）」と置くと、式（４）により表わされ、
前記Ｙ→Ｘ誤差は、時刻ｔの誤差を「ε(y)（ｔ）」と置くと、式（５）により表わされ、
前記ＹＸ→Ｘ誤差は、時刻ｔの誤差を「ε(y,x)（ｔ）」と置くと、式（６）により表わされる
因果関係判定装置。
【数２】
2008257595.tif 000026
ここで、式（３）〜式（６）中、「α」，「Ａ」，「Ｂ」，「β」，「Ｃ」，「Ｄ」は係数ベクトル，「π」，「θ」，「ｐ」，「ｑ」は定数，「Ｘt-1」，「Ｙt-1」は前記Ｘｄａｔａ，前記Ｙｄａｔａの時刻ｔ−１の状態ベクトルである。
【請求項６】
請求項５記載の因果関係判定装置であって、
前記状態ベクトルＸt-1は、前記Ｘｄａｔａの前記時刻ｔ−１以前の複数の時刻の入力値を要素の少なくとも一部とするベクトルであり、
前記状態ベクトルＹt-1は、前記Ｙｄａｔａの前記時刻ｔ−１以前の複数の時刻の入力値を要素の少なくとも一部とするベクトルである
因果関係判定装置。
【請求項７】
請求項６記載の因果関係判定装置であって、
前記状態ベクトルＸt-1は、式（７）により表わされるベクトルであり、
前記状態ベクトルＹt-1は、式（８）により表わされるベクトルである
因果関係判定装置。
【数３】
2008257595.tif 000027
【数４】
2008257595.tif 000028
ここで、前記式（７），（８）中、「ｍ」，「ｎ」は埋め込み次元を示し、「τ」，「ν」は状態ベクトルの各要素間の時間ずれ（ラグ）であり、「ｏ1」，「ｏ2」は項の次数の最大値である。
【請求項８】
請求項５ないし７いずれか記載の因果関係判定装置であって、
前記Ｙ→Ｘ評価式は、「ｃ1」と置くと、式（９）により表わされ、
前記Ｘ→Ｙ評価式は、「ｃ2」と置くと、式（１０）により表わされる
因果関係判定装置。
【数５】
2008257595.tif 000029
ここで、前記式（９）中「ｗy→x（ｔ）」は前記Ｙ→Ｘ重みであり、前記式（１０）中「ｗx→y（ｔ）」は前記Ｘ→Ｙ重みである。
【請求項９】
前記重み演算処理手段は、任意の時刻に対して、前記Ｘ→Ｘ誤差ε(x)（ｔ）の大きさに比して前記ＸＹ→Ｘ誤差ε(x,y)（ｔ）の大きさが小さい時刻ほど前記Ｙ→Ｘ重みｗy→x（ｔ）が大きくなるよう各時刻ｔの前記Ｙ→Ｘ重みｗy→x（ｔ）を演算し、前記Ｙ→Ｙ誤差ε(y)（ｔ）の大きさに比して前記ＹＸ→Ｙ誤差ε(y,x)（ｔ）の大きさが小さい時刻ほど前記Ｘ→Ｙ重みｗx→y（ｔ）が大きくなるよう各時刻ｔの前記Ｘ→Ｙ重みｗx→y（ｔ）を演算する手段である請求項８記載の因果関係判定装置。
【請求項１０】
請求項５ないし７いずれか記載の因果関係判定装置であって、
前記Ｙ→Ｘ評価式は、「ｃ3」と置くと、式（１１）により表わされ、
前記Ｘ→Ｙ評価式は、「ｃ4」と置くと、式（１２）により表わされる
因果関係判定装置。
【数６】
2008257595.tif 000030
ここで、前記式（１１）中、「ｗy→x（ｔ）」は前記Ｙ→Ｘ重みであり、「λ1」はパラメータであり、「Ｅ（ｗy→x）」は式（１３）により表わされる。また、前記式（１２）中、「ｗx→y（ｔ）」は前記Ｘ→Ｙ重みであり、「λ2」はパラメータであり、「Ｅ（ｗx→y）」は式（１４）により表わされる。
【数７】
2008257595.tif 000031
ここで、前記式（１３），（１４）中、「ｅ」はパラメータである。
【請求項１１】
前記重み演算処理手段は、任意の時刻の前記Ｙ→Ｘ重みｗy→x（ｔ）に第１の所定値を設定し、該設定した各時刻ｔの重みｗy→x（ｔ）を用いて式（１５）および式（１６）により前記係数ベクトルα，Ａ，Ｂおよび前記定数π，ｐを演算する処理と、該演算した係数ベクトルα，Ａ，Ｂおよび定数π，ｐを用いて前記式（３）および前記式（４）によりＸ→Ｘ誤差ε(x)（ｔ）およびＸＹ→Ｘ誤差ε(x,y)（ｔ）を演算する処理と、該演算したＸ→Ｘ誤差ε(x)（ｔ）およびＸＹ→Ｘ誤差ε(x,y)（ｔ）を用いて前記Ｙ→Ｘ評価式の値が最大となるよう各時刻ｔの前記Ｙ→Ｘ重みｗy→x（ｔ）を演算する処理と、の３処理を繰り返し実行して各時刻ｔの前記Ｙ→Ｘ重みｗy→x（ｔ）を演算し、任意の時刻の前記Ｘ→Ｙ重みｗx→y（ｔ）に第２の所定値を設定し、該設定した各時刻ｔの重みｗx→y（ｔ）を用いて式（１７）および式（１８）により前記係数ベクトルβ，Ｃ，Ｄおよび前記定数θ，ｑを演算する処理と、該演算した係数ベクトルβ，Ｃ，Ｄおよび定数θ，ｑを用いて前記式（５）および前記式（６）によりＹ→Ｙ誤差ε(y)（ｔ）およびＹＸ→Ｙ誤差ε(y,x)（ｔ）を演算する処理と、該演算したＹ→Ｙ誤差ε(y)（ｔ）およびＹＸ→Ｙ誤差ε(y,x)（ｔ）を用いて前記Ｘ→Ｙ評価式の値が最大となるよう各時刻ｔの前記Ｘ→Ｙ重みｗx→y（ｔ）を演算する処理と、の３処理を繰り返し実行して各時刻ｔの前記Ｘ→Ｙ重みｗx→y（ｔ）を演算する手段である請求項８ないし１０いずれか記載の因果関係判定装置。
【数８】
2008257595.tif 000032
【数９】
2008257595.tif 000033
【請求項１２】
前記二つの異なる時系列データは、気象に関する時系列データである請求項１ないし１１いずれか記載の因果関係判定装置。
【請求項１３】
前記二つの異なる時系列データは、生体に関する時系列データである請求項１ないし１１いずれか記載の因果関係判定装置。
【請求項１４】
情報を入力する入力手段と、情報を記憶する記憶手段と、情報を出力する出力手段と、演算処理を行なう演算処理装置と、を備えるコンピュータを二つの異なる時系列データ間の各時刻の因果関係の有無を判定する因果関係判定装置として機能させるプログラムであって、
前記二つの異なる時系列データを前記入力手段により入力すると共に該入力した二つの異なる時系列データを前記記憶手段に記憶する入力処理モジュールと、
前記記憶手段に記憶された二つの異なる時系列データに基づく該二つの異なる時系列データ間の各時刻の前記因果関係の程度と該因果関係の程度に対する各時刻の重みとを用いて所定の演算手法により演算される結果が最大となるよう各時刻の前記重みを演算し、該演算結果を前記記憶手段に記憶する重み演算処理モジュールと、
前記記憶手段に記憶された各時刻の重みを用いて該記憶手段に記憶された二つの異なる時系列データ間の各時刻の前記因果関係の有無を判定し、該判定結果を前記出力手段により出力する判定出力処理モジュールと、
を備えるプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、因果関係判定装置およびプログラムに関し、詳しくは、二つの異なる時系列データ間の各時刻の因果関係の有無を判定する因果関係判定装置およびコンピュータを因果関係判定装置として機能させるプログラムに関する。
【背景技術】
【０００２】
従来、この種の因果関係判定装置としては、二つの異なる時系列データ間で因果関係が存在し続けるものに対して、線形回帰モデル（例えば、非特許文献１参照）やこれを非線形に拡張したモデル（例えば、非特許文献２参照）を用いて因果関係の程度を求めるものが提案されている。
【非特許文献１】北川源四郎．時系列解析入門．岩波書店，２００５．
【非特許文献２】N. Ancona, D. Marinazzo, and S. Stramaglia. Radial basis function approach to nonlinear Granger causality of time series. Physical Review E, Vol. 70, p. 056221, 2004.
【発明の開示】
【発明が解決しようとする課題】
【０００３】
しかしながら、一般に、現実世界の事象から得られる二つの時系列データ間では、各時刻に因果関係が存在したり消滅したりするものが多いため、こうした二つの時系列データ間の各時刻の因果関係の有無を判定できるようにすることが望まれている。
【０００４】
本発明の因果関係判定装置およびプログラムは、二つの異なる時系列データ間で因果関係が存在しているのか否かを判定できるようにすること、特に、気象に関する時系列データや生体に関する時系列データ等の自然界における二つの異なる時系列データ間で因果関係が存在しているのか否かを判定できるようにすることを目的とする。
【課題を解決するための手段】
【０００５】
本発明の因果関係判定装置およびプログラムは、上述の目的を達成するために以下の手段を採った。
【０００６】
本発明の因果関係判定装置は、
二つの異なる時系列データ間の各時刻の因果関係の有無を判定する因果関係判定装置であって、
情報を入力する入力手段と、
情報を記憶する記憶手段と、
情報を出力する出力手段と、
前記二つの異なる時系列データを前記入力手段により入力すると共に該入力した二つの異なる時系列データを前記記憶手段に記憶する入力処理手段と、
前記記憶手段に記憶された二つの異なる時系列データに基づく該二つの異なる時系列データ間の各時刻の前記因果関係の程度と該因果関係の程度に対する各時刻の重みとを用いて所定の演算手法により演算される結果が最大となるよう各時刻の前記重みを演算し、該演算結果を前記記憶手段に記憶する重み演算処理手段と、
前記記憶手段に記憶された各時刻の重みを用いて該記憶手段に記憶された二つの異なる時系列データ間の各時刻の前記因果関係の有無を判定し、該判定結果を前記出力手段により出力する判定出力処理手段と、
を備えることを要旨とする。
【０００７】
この本発明の因果関係判定装置では、まず、二つの異なる時系列データを入力手段により入力すると共に入力した二つの異なる時系列データを記憶手段に記憶する。続いて、記憶手段に記憶された二つの異なる時系列データに基づく二つの異なる時系列データ間の各時刻の因果関係の程度とその因果関係の程度に対する各時刻の重みとを用いて所定の演算手法により演算される結果が最大となるよう各時刻の重みを演算すると共に演算結果を記憶手段に記憶する。そして、記憶手段に記憶された各時刻の重みを用いて記憶手段に記憶された二つの異なる時系列データ間の各時刻の因果関係の有無を判定すると共に判定結果を出力手段により出力する。これにより、二つの異なる時系列データ間の各時刻の因果関係の程度に対する重みを用いてこれらの間で各時刻に因果関係が存在しているか否かを判定することができる。
【０００８】
こうした本発明の因果関係判定装置において、前記二つの異なる時系列データは、Ｘｄａｔａ，Ｙｄａｔａと置くと、式（１）および式（２）により表わされるものであり、前記重み演算処理手段は、任意の時刻における前記Ｙｄａｔａから前記Ｘｄａｔａへの各時刻の前記因果関係の程度であるＹ→Ｘ程度と該Ｙｄａｔａから該Ｘｄａｔａへの各時刻の前記重みであるＹ→Ｘ重みとを用いた評価式であるＹ→Ｘ評価式の値が最大となるよう該Ｙｄａｔａから該Ｘｄａｔａへの各時刻の該重みを演算すると共に該演算結果を前記記憶手段に記憶し、任意の時刻における前記Ｘｄａｔａから前記Ｙｄａｔａへの各時刻の前記因果関係の程度であるＸ→Ｙ程度と該Ｘｄａｔａから該Ｙｄａｔａへの各時刻の前記重みであるＸ→Ｙ重みとを用いた評価式であるＸ→Ｙ評価式の値が最大となるよう該Ｘｄａｔａから該Ｙｄａｔａへの各時刻の該重みを演算すると共に該演算結果を前記記憶手段に記憶する手段であるものとすることもできる。ここで、式（１）中、ｘ（１），…，ｘ（Ｎ）は、時刻１，…，ＮのＸｄａｔａの値である。また、式（２）中、ｙ（１），…，ｙ（Ｎ）は、時刻１，…，ＮのＹｄａｔａの値である。
【０００９】
【数１】
2008257595.tif 000003
【００１０】
この二つの異なる時系列データをＸｄａｔａ，Ｙｄａｔａと置いた態様の本発明の因果関係判定装置において、前記Ｙ→Ｘ程度は、各時刻に対して、前記Ｘｄａｔａの第１の時刻の入力値と該Ｘｄａｔａの該第１の時刻よりも前の時刻の１以上の入力値を用いて予測される該Ｘｄａｔａの該第１の時刻の予測値との誤差であるＸ→Ｘ誤差と、前記Ｘｄａｔａの前記第１の時刻の入力値と該Ｘｄａｔａおよび前記Ｙｄａｔａの該第１の時刻よりも前の時刻の１以上の入力値を用いて予測される該Ｘｄａｔａの該第１の時刻の予測値との誤差であるＸＹ→Ｘ誤差と、を用いて表わされるものであり、前記Ｘ→Ｙ程度は、各時刻に対して、前記Ｙｄａｔａの前記第１の時刻の入力値と該Ｙｄａｔａの該第１の時刻よりも前の時刻の１以上の入力値を用いて予測される該Ｙｄａｔａの該第１の時刻の予測値との誤差であるＹ→Ｙ誤差と、前記Ｙｄａｔａの前記第１の時刻の入力値と前記Ｘｄａｔａおよび該Ｙｄａｔａの該第１の時刻よりも前の時刻の１以上の入力値を用いて予測される該Ｙｄａｔａの該第１の時刻の予測値との誤差であるＹＸ→Ｙ誤差と、を用いて表わされるものであるものとすることもできる。こうすれば、Ｘ→Ｘ誤差やＸＹ→Ｘ誤差を用いてＹｄａｔａからＸｄａｔａへの各時刻の重みを演算することができ、Ｙ→Ｙ誤差やＹＸ→Ｙ誤差を用いてＸｄａｔａからＹｄａｔａへの各時刻の重みを演算することができる。
【００１１】
このＸ→Ｘ誤差やＸＹ→Ｘ誤差を用いてＹ→Ｘ程度が表わされると共にＹ→Ｙ誤差やＹＸ→Ｙ誤差を用いてＸ→Ｙ程度が表わされる態様の本発明の因果関係判定装置において、前記Ｙ→Ｘ評価式は、前記Ｙ→Ｘ重みと、前記Ｘ→Ｘ誤差の２乗と前記ＸＹ→Ｘ誤差の２乗との差と、の積の期待値であるＹ→Ｘ期待値または該Ｙ→Ｘ期待値に第１の補正項を加えた値を演算する式であり、前記Ｘ→Ｙ評価式は、前記Ｘ→Ｙ重みと、前記Ｙ→Ｙ誤差の２乗と前記ＹＸ→Ｙ誤差の２乗との差と、の積の期待値であるＸ→Ｙ期待値または該Ｘ→Ｙ期待値に第２の補正項を加えた値を演算する式であるものとすることもできる。
【００１２】
また、Ｘ→Ｘ誤差やＸＹ→Ｘ誤差を用いてＹ→Ｘ程度が表わされると共にＹ→Ｙ誤差やＹＸ→Ｙ誤差を用いてＸ→Ｙ程度が表わされる態様の本発明の因果関係判定装置において、前記Ｘ→Ｘ誤差は、時刻ｔの誤差を「ε(x)（ｔ）」と置くと、式（３）により表わされ、前記ＸＹ→Ｘ誤差は、時刻ｔの誤差を「ε(x,y)（ｔ）」と置くと、式（４）により表わされ、前記Ｙ→Ｘ誤差は、時刻ｔの誤差を「ε(y)（ｔ）」と置くと、式（５）により表わされ、前記ＹＸ→Ｘ誤差は、時刻ｔの誤差を「ε(y,x)（ｔ）」と置くと、式（６）により表わされるものとすることもできる。ここで、式（３）〜式（６）中、「α」，「Ａ」，「Ｂ」，「β」，「Ｃ」，「Ｄ」は係数ベクトル，「π」，「θ」，「ｐ」，「ｑ」は定数，「Ｘt-1」，「Ｙt-1」は前記Ｘｄａｔａ，前記Ｙｄａｔａの時刻ｔ−１の状態ベクトルである。
【００１３】
【数２】
2008257595.tif 000004
【００１４】
この各誤差が式（３）〜式（６）により表わされる態様の本発明の因果関係判定装置において、前記状態ベクトルＸt-1は、前記Ｘｄａｔａの前記時刻ｔ−１以前の複数の時刻の入力値を要素の少なくとも一部とするベクトルであり、前記状態ベクトルＹt-1は、前記Ｙｄａｔａの前記時刻ｔ−１以前の複数の時刻の入力値を要素の少なくとも一部とするベクトルであるものとすることもできる。この場合、前記状態ベクトルＸt-1は、式（７）により表わされるベクトルであり、前記状態ベクトルＹt-1は、式（８）により表わされるベクトルであるものとすることもできる。ここで、前記式（７），（８）中、「ｍ」，「ｎ」は埋め込み次元を示し、「τ」，「ν」は状態ベクトルの各要素間の時間ずれ（ラグ）であり、「ｏ1」，「ｏ2」は項の次数の最大値である。
【００１５】
【数３】
2008257595.tif 000005
【００１６】
【数４】
2008257595.tif 000006
【００１７】
また、各誤差が式（３）〜式（６）により表わされる態様の本発明の因果関係判定装置において、前記Ｙ→Ｘ評価式は、「ｃ1」と置くと、式（９）により表わされ、前記Ｘ→Ｙ評価式は、「ｃ2」と置くと、式（１０）により表わされるものとすることもできる。ここで、前記式（９）中「ｗy→x（ｔ）」は前記Ｙ→Ｘ重みであり、前記式（１０）中「ｗx→y（ｔ）」は前記Ｘ→Ｙ重みである。この場合、前記重み演算処理手段は、任意の時刻に対して、前記Ｘ→Ｘ誤差ε(x)（ｔ）の大きさに比して前記ＸＹ→Ｘ誤差ε(x,y)（ｔ）の大きさが小さい時刻ほど前記Ｙ→Ｘ重みｗy→x（ｔ）が大きくなるよう各時刻ｔの前記Ｙ→Ｘ重みｗy→x（ｔ）を演算し、前記Ｙ→Ｙ誤差ε(y)（ｔ）の大きさに比して前記ＹＸ→Ｙ誤差ε(y,x)（ｔ）の大きさが小さい時刻ほど前記Ｘ→Ｙ重みｗx→y（ｔ）が大きくなるよう各時刻ｔの前記Ｘ→Ｙ重みｗx→y（ｔ）を演算する手段であるものとすることもできる。
【００１８】
【数５】
2008257595.tif 000007
【００１９】
さらに、各誤差が式（３）〜式（６）により表わされる態様の本発明の因果関係判定装置において、前記Ｙ→Ｘ評価式は、「ｃ3」と置くと、式（１１）により表わされ、前記Ｘ→Ｙ評価式は、「ｃ4」と置くと、式（１２）により表わされるものとすることもできる。ここで、前記式（１１）中、「ｗy→x（ｔ）」は前記Ｙ→Ｘ重みであり、「λ1」はパラメータであり、「Ｅ（ｗy→x）」は式（１３）により表わされる。また、前記式（１２）中、「ｗx→y（ｔ）」は前記Ｘ→Ｙ重みであり、「λ2」はパラメータであり、「Ｅ（ｗx→y）」は式（１４）により表わされる。ここで、前記式（１３），（１４）中、「ｅ」はパラメータである。
【００２０】
【数６】
2008257595.tif 000008
【００２１】
【数７】
2008257595.tif 000009
【００２２】
このＹ→Ｘ評価式が式（９）または式（１１），Ｘ→Ｙ評価式が式（１０）または（１２）により表わされる態様の本発明の因果関係判定装置において、前記重み演算処理手段は、任意の時刻の前記Ｙ→Ｘ重みｗy→x（ｔ）に第１の所定値を設定し、該設定した各時刻ｔの重みｗy→x（ｔ）を用いて式（１５）および式（１６）により前記係数ベクトルα，Ａ，Ｂおよび前記定数π，ｐを演算する処理と、該演算した係数ベクトルα，Ａ，Ｂおよび定数π，ｐを用いて前記式（３）および前記式（４）によりＸ→Ｘ誤差ε(x)（ｔ）およびＸＹ→Ｘ誤差ε(x,y)（ｔ）を演算する処理と、該演算したＸ→Ｘ誤差ε(x)（ｔ）およびＸＹ→Ｘ誤差ε(x,y)（ｔ）を用いて前記Ｙ→Ｘ評価式の値が最大となるよう各時刻ｔの前記Ｙ→Ｘ重みｗy→x（ｔ）を演算する処理と、の３処理を繰り返し実行して各時刻ｔの前記Ｙ→Ｘ重みｗy→x（ｔ）を演算し、任意の時刻の前記Ｘ→Ｙ重みｗx→y（ｔ）に第２の所定値を設定し、該設定した各時刻ｔの重みｗx→y（ｔ）を用いて式（１７）および式（１８）により前記係数ベクトルβ，Ｃ，Ｄおよび前記定数θ，ｑを演算する処理と、該演算した係数ベクトルβ，Ｃ，Ｄおよび定数θ，ｑを用いて前記式（５）および前記式（６）によりＹ→Ｙ誤差ε(y)（ｔ）およびＹＸ→Ｙ誤差ε(y,x)（ｔ）を演算する処理と、該演算したＹ→Ｙ誤差ε(y)（ｔ）およびＹＸ→Ｙ誤差ε(y,x)（ｔ）を用いて前記Ｘ→Ｙ評価式の値が最大となるよう各時刻ｔの前記Ｘ→Ｙ重みｗx→y（ｔ）を演算する処理と、の３処理を繰り返し実行して各時刻ｔの前記Ｘ→Ｙ重みｗx→y（ｔ）を演算する手段であるものとすることもできる。
【００２３】
【数８】
2008257595.tif 000010
【００２４】
【数９】
2008257595.tif 000011
【００２５】
本発明の因果関係判定装置において、前記二つの異なる時系列データは、気象に関する時系列データであるものとすることもできるし、生体に関する時系列データであるものとすることもできる。ここで、二つの異なる時系列データが気象に関する時系列データである場合には、例えば、一地点の気温と気圧の時系列データや、二地点の気温と気温の時系列データ，二地点の気温と気圧の時系列データであるものとすることができる。また、二つの異なる時系列データが生体に関する時系列データである場合には、例えば、心拍数と呼吸の時系列データであるものとすることができる。さらに、二つの異なる時系列データは、自然界の時系列データであるものとすることもできる。例えば、医療解析に関する時系列データ（例えば、薬剤の投与量と体温の時系列データや、薬剤の投与量と血中のある成分の濃度の時系列データなど）であるものとすることもできるし、通信ネットワークに関する時系列データ（例えば、サーバの負荷と通信速度の時系列データなど）であるものとすることもできるし、環境アセスメントに関する時系列データ（例えば、工事箇所の数と環境指標の値の時系列データなど）であるものとすることもできるし、神経科学的分析に関する時系列データ（例えば、脳の二点の神経信号の時系列データ（電圧など）など）であるものとすることもできるし、遺伝子ネットワーク解析に関する時系列データ（例えば、二つの異なる遺伝子の発現状態を示すＪ系列データなど）であるものとすることもできるし、行動認識に関する時系列データ（例えば、画像データと人の位置の時系列データや、二人の位置の時系列データなど）であるものとすることもできるし、自動車運転行動分析に関する時系列データ（例えば、車間距離とアクセル操作量の時系列データや、歩行者の位置とハンドル操作量の時系列データ，他車の位置と視線方向の時系列データなど）であるものとすることもできる。あるいは、二つの異なる時系列データは、自然界以外の時系列データであるものとすることもできる。
【００２６】
本発明のプログラムは、
情報を入力する入力手段と、情報を記憶する記憶手段と、情報を出力する出力手段と、演算処理を行なう演算処理装置と、を備えるコンピュータを二つの異なる時系列データ間の各時刻の因果関係の有無を判定する因果関係判定装置として機能させるプログラムであって、
前記二つの異なる時系列データを前記入力手段により入力すると共に該入力した二つの異なる時系列データを前記記憶手段に記憶する入力処理モジュールと、
前記記憶手段に記憶された二つの異なる時系列データに基づく該二つの異なる時系列データ間の各時刻の前記因果関係の程度と該因果関係の程度に対する各時刻の重みとを用いて所定の演算手法により演算される結果が最大となるよう各時刻の前記重みを演算し、該演算結果を前記記憶手段に記憶する重み演算処理モジュールと、
前記記憶手段に記憶された各時刻の重みを用いて該記憶手段に記憶された二つの異なる時系列データ間の各時刻の前記因果関係の有無を判定し、該判定結果を前記出力手段により出力する判定出力処理モジュールと、
を備えることを要旨とする。
【００２７】
この本発明のプログラムでは、まず、入力処理モジュールにより、二つの異なる時系列データを入力手段により入力すると共に入力した二つの異なる時系列データを記憶手段に記憶する。続いて、重み演算処理モジュールにより、記憶手段に記憶された二つの異なる時系列データに基づく二つの異なる時系列データ間の各時刻の因果関係の程度とその因果関係の程度に対する各時刻の重みとを用いて所定の演算手法により演算される結果が最大となるよう各時刻の重みを演算すると共に演算結果を記憶手段に記憶する。そして、判定出力処理モジュールにより、記憶手段に記憶された各時刻の重みを用いて記憶手段に記憶された二つの異なる時系列データ間の各時刻の因果関係の有無を判定すると共に判定結果を出力手段により出力する。これにより、二つの異なる時系列データ間の各時刻の因果関係の程度に対する重みを用いてこれらの間で各時刻に因果関係が存在しているか否かを判定することができる。
【発明を実施するための最良の形態】
【００２８】
次に、本発明を実施するための最良の形態を実施例を用いて説明する。
【実施例】
【００２９】
図１は、本発明の一実施例としての因果関係判定装置２０の構成の概略を示す構成図である。実施例の因果関係判定装置２０は、図示するように、一般的な汎用のコンピュータ２２にアプリケーションソフトウエアとしての二つの異なる時系列データ（以下、Ｘｄａｔａ，Ｙｄａｔａという）間の各時刻の因果関係の有無を判定する因果関係判定プログラム３０がインストールされたものとして構成されている。
【００３０】
コンピュータ２２は、周知のＣＰＵ２３，ＲＯＭ２４，ＲＡＭ２５，クロック回路２６，ハードディスク２７，入出力処理回路２８などを備え、これらがバス２９によって接続されている。ハードディスク２７には、アプリケーションソフトウエアとしてインストールされた因果関係判定プログラム３０が記憶されている。この因果関係判定プログラム３０は、二つの異なる時系列データを式（１９），（２０）により表わされるＸｄａｔａ，Ｙｄａｔａとして入力すると共に入力したＸｄａｔａ，ＹｄａｔａをＲＡＭ２５に記憶する入力処理モジュール３２と、ＲＡＭ２５に記憶されたＸｄａｔａ，Ｙｄａｔａに基づくＹｄａｔａからＸｄａｔａへの各時刻ｔ（１≦ｔ≦Ｎ）の因果関係の程度とその因果関係の程度に対する各時刻ｔの重みであるＹ→Ｘ重みｗy→x（ｔ）とを用いて各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）を演算すると共に演算結果をＲＡＭ２５に記憶するＹ→Ｘ重み演算処理モジュール３３と、ＲＡＭ２５に記憶されたＹ→Ｘ重みｗy→x（ｔ）を用いてＹｄａｔａからＸｄａｔａへの各時刻ｔの因果関係の有無を判定すると共に判定結果を出力するＹ→Ｘ判定出力処理モジュール３４と、ＲＡＭ２５に記憶されたＸｄａｔａ，Ｙｄａｔａに基づくＸｄａｔａからＹｄａｔａへの各時刻ｔの因果関係の程度とその因果関係の程度に対する各時刻ｔの重みであるＸ→Ｙ重みｗx→y（ｔ）とを用いて各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）を演算すると共に演算結果をＲＡＭ２５に記憶するＸ→Ｙ重み演算処理モジュール３５と、ＲＡＭ２５に記憶されたＸ→Ｙ重みｗx→y（ｔ）を用いてＸｄａｔａからＹｄａｔａへの各時刻ｔの因果関係の有無を判定すると共に判定結果を出力するＸ→Ｙ判定出力処理モジュール３６と、から構成されている。なお、式（１９），（２０）中、「Ｎ」は、Ｘｄａｔａ，Ｙｄａｔａのデータ数であり、ｘ（１），…，ｘ（Ｎ）は、時刻１，…，ＮのＸｄａｔａの値であり、ｙ（１），…ｙ（Ｎ）は、時刻１，…，ＮのＹｄａｔａの値である。
【００３１】
【数１０】
2008257595.tif 000012
【００３２】
このコンピュータ２２には、ＣＤやＤＶＤに対して読み書き可能なディスクドライブ装置４１が入出力処理回路２８を介してデータを入出力可能に接続され、キーボード４２及びマウス４３が入出力処理回路２８を介してデータを入力可能に接続され、情報を出力表示するディスプレイ４４が入出力処理回路２８を介してデータを出力可能に接続されている。
【００３３】
次に、こうして構成される因果関係判定装置２０の動作、特にユーザによるキーボード４２やマウス４３の操作により因果関係判定プログラム３０の実行が指示されたときの動作について説明する。図２は、ＣＰＵ２３により実行される因果関係判定プログラム３０の一例を示すフローチャートである。この因果関係判定プログラム３０は、ユーザにより実行が指示されたときに、ＣＰＵ２３により、ＲＡＭ２５の所定アドレスに書き込まれる。そして、ＣＰＵ２３は、ＲＡＭ２５に書き込まれた因果関係判定プログラム３０を読み込んで実行する。
【００３４】
因果関係判定プログラム３０が実行されると、ＣＰＵ２３は、まず、前述の式（１９），（２０）により表わされるＸｄａｔａ，Ｙｄａｔａを入力する処理を実行する（ステップＳ１００）。ここで、Ｘｄａｔａ，Ｙｄａｔａの入力は、例えば、ユーザのキーボード４２やマウス４３の操作に応じて、ディスクドライブ装置４１にセットされたＣＤやＤＶＤに書き込まれているデータやウェブサイトからインターネットを介してダウンロードされたデータなどをＲＡＭ２５の所定アドレスに書き込むと共にこれを読み込むことにより入力することができる。
【００３５】
こうしてデータを入力すると、入力したＸｄａｔａ，Ｙｄａｔａを用いて図３に例示するＹ→Ｘ重み演算処理によりＹｄａｔａからＸｄａｔａへの各時刻ｔの因果関係の程度に対する各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）を演算すると共に演算結果をＲＡＭ２５の所定アドレスに記憶する（ステップＳ１１０）。Ｙ→Ｘ重み演算処理は、実施例では、式（２１）および式（２２）により表わされる回帰モデルと各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）とを用いて行なうものとした。ここで、式（２１）および式（２２）中、「α」，「Ａ」，「Ｂ」は係数ベクトル，「π」，「ｐ」は定数，「ε(x)（ｔ）」はＸｄａｔａの時刻ｔの入力値ｘ（ｔ）とＸｄａｔａの時刻ｔ−１の状態ベクトルＸt-1を用いて予測される時刻ｔの予測値との誤差（以下、Ｘ→Ｘ誤差という），「ε(x,y)（ｔ）」はＸｄａｔａの時刻ｔの入力値ｘ（ｔ）とＸｄａｔａおよびＹｄａｔａの時刻ｔ−１の状態ベクトルＸt-1，Ｙt-1を用いて予測される時刻ｔの予測値との誤差（以下、ＸＹ→Ｘ誤差という）である。状態ベクトルＸt-1，Ｙt-1は、式（２３）および式（２４）により表わされる。式（２３），（２４）中、「ｍ」，「ｎ」は埋め込み次元を示し、「τ」，「ν」は状態ベクトルの各要素間の時間ずれ（ラグ）であり、「ｏ1」，「ｏ2」は項の次数の最大値である。この「ｍ」，「ｎ」，「τ」，「ν」，「ｏ1」，「ｏ2」は、ユーザによるキーボード４２やマウス４３の操作に応じてＲＡＭ２５の所定アドレスに書き込まれたものを読み込んで用いるものとしてもよいし、入力されたＸｄａｔａ，Ｙｄａｔａの種類に応じて図示しないパラメータ設定ルーチンにより設定されてＲＡＭ２５の所定アドレスに書き込まれたものを読み込んで用いるものとしてもよいし、固定値として因果関係プログラム３０と共にＲＡＭ２５の所定アドレスに書き込まれたものを読み込んで用いるものとしてもよい。このような状態ベクトルＸt，Ｙtを用いることにより、式（２５）および式（２６）により表わされる状態ベクトルＸt，Ｙtを用いるものに比して考慮される要素が多くなると共に非線形な要素を導入することができるため、時刻ｔの予測値をより適正に演算することができる。以下、図２の因果関係判定ルーチンの説明を一旦中断し、図３のＹ→Ｘ重み演算処理について説明する。
【００３６】
【数１１】
2008257595.tif 000013
【００３７】
【数１２】
2008257595.tif 000014
【００３８】
【数１３】
2008257595.tif 000015
【００３９】
【数１４】
2008257595.tif 000016
【００４０】
Ｙ→Ｘ重み演算処理では、まず、任意の時刻に対して各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）に値１を設定する（ステップＳ２００）。そして、式（２１）および式（２２）により表わされる回帰モデルと各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）とを用いた式（２７）および式（２８）により係数ベクトルα，Ａ，Ｂおよび定数π，ｐを演算する（ステップＳ２１０）。式（２７）は、各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）とＸ→Ｘ誤差ε(x)（ｔ）の２乗との積の期待値が最小となるよう係数ベクトルαおよび定数πを演算する式である。この式（２７）では、Ｙ→Ｘ重みｗy→x（ｔ）とＸｄａｔａの入力値ｘ（ｔ）と状態ベクトルＸt-1が定まっているため、Ｘ→Ｘ誤差ε(x)（ｔ）の２乗の期待値が最小となるよう係数ベクトルαおよび定数πを演算することになる。また、式（２８）は、Ｙ→Ｘ重みｗy→x（ｔ）とＸＹ→Ｘ誤差ε(x,y)（ｔ）の２乗との積の期待値が最小となるよう係数ベクトルＡ，Ｂおよび定数ｐを演算する式である。この式（２８）では、各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）とＸｄａｔａの入力値ｘ（ｔ）と状態ベクトルＸt-1，Ｙt-1が定まっているため、ＸＹ→Ｘ誤差ε(x,y)（ｔ）の期待値が最小となるよう係数ベクトルＡ，Ｂおよび定数ｐを演算することになる。
【００４１】
【数１５】
2008257595.tif 000017
【００４２】
こうして係数ベクトルα，Ａ，Ｂおよび定数π，ｐを演算すると、演算した係数ベクトルαおよび定数πとＸｄａｔａの入力値ｘ（ｔ）とを用いて前述の式（２１）により各時刻ｔのＸ→Ｘ誤差ε(x)（ｔ）を演算すると共に係数ベクトルＡ，Ｂおよび定数ｐとＸｄａｔａの入力値ｘ（ｔ）とを用いて式（２２）により各時刻ｔのＸＹ→Ｘ誤差ε(x,y)（ｔ）を演算する（ステップＳ２２０）。
【００４３】
そして、各時刻ｔのＸ→Ｘ誤差ε(x)（ｔ）およびＸＹ→Ｘ誤差ε(x,y)（ｔ）を用いて式（２９）の評価式ｃ5の値が最大となるよう各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）を演算する（ステップＳ２３０）。ここで、式（２９）中、Ｘ→Ｘ誤差ε(x)（ｔ）やＸＹ→Ｘ誤差ε(x,y)（ｔ）を含む項は、各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）と、Ｘ→Ｘ誤差ε(x)（ｔ）の２乗とＸＹ→Ｘ誤差ε(x,y)（ｔ）の２乗との差と、の積の期待値である。また、式（２９）中、「λ1×Ｅ（ｗy→x）」はＹ→Ｘ重みｗy→x（ｔ）が時刻の変化に対して激しく振動しないようにするための振動抑制項である。そして、「Ｅ（ｗy→x）」は、式（３０）により表わされ、式（３０）中、「ｅ」はパラメータである。また、「λ1」はパラメータであり、ユーザによるキーボード４２やマウス４３の操作に応じてＲＡＭ２５の所定アドレスに書き込まれたものを読み込んで用いるものとしてもよいし、入力されたＸｄａｔａ，Ｙｄａｔａの種類に応じて図示しないパラメータ設定ルーチンにより設定されてＲＡＭ２５の所定アドレスに書き込まれたものを読み込んで用いるものとしてもよいし、固定値として因果関係プログラム３０と共にＲＡＭ２５の所定アドレスに書き込まれたものを読み込んで用いるものとしてもよい。式（２９）より、Ｘ→Ｘ誤差ε(x)（ｔ）の大きさに比してＸＹ→Ｘ誤差ε(x,y)（ｔ）の大きさが小さいほど｛（ε(x)（ｔ））2−（ε(x,y)（ｔ））2｝は大きくなることが解る。ここで、Ｘ→Ｘ誤差ε(x)（ｔ）の大きさに比してＸＹ→Ｘ誤差ε(x,y)（ｔ）の大きさが小さいとは、状態ベクトルＸt-1を考慮して得られる時刻ｔの予測値に比して状態ベクトルＸt-1，Ｙt-1を考慮して得られる時刻ｔの予測値の方がＸｄａｔａの時刻ｔの入力値ｘ（ｔ）に近い、即ち、状態ベクトルＹt-1がＸｄａｔａの時刻ｔの入力値ｘ（ｔ）に影響を与えている、ということを意味する。したがって、Ｘ→Ｘ誤差ε(x)（ｔ）の大きさに比してＸＹ→Ｘ誤差ε(x,y)（ｔ）の大きさが小さいほど、状態ベクトルＹt-1がＸｄａｔａの時刻ｔの入力値ｘ（ｔ）に与える影響が大きい、即ち、時刻ｔでのＹｄａｔａからＸｄａｔａへの因果関係の程度が大きいといえる。このことと、各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）の総和を固定値（値Ｎ）としたこととを踏まえると、ステップＳ２３０の評価式ｃ5の値が最大となるよう各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）を演算する処理は、式（２９）のパラメータλ1が値０のとき、即ち振動抑制項（λ1×Ｅ（ｗy→x））を無視するときを考えれば、ＹｄａｔａからＸｄａｔａへの因果関係の程度が大きい時刻ｔほど大きなＹ→Ｘ重みｗy→x（ｔ）を演算し、ＹｄａｔａからＸｄａｔａへの因果関係の程度が小さい時刻ｔほど小さなＹ→Ｘ重みｗy→x（ｔ）を演算する処理となる。次に、式（２９）において、パラメータλ1が値０でないとき、即ち、振動抑制項（λ1×Ｅ（ｗy→x））を無視しないときを考える。Ｘｄａｔａ，Ｙｄａｔａによっては｛（ε(x)（ｔ））2−（ε(x,y)（ｔ））2｝が時刻の変化に対して激しく振動することがあり、この場合、振動抑制項（λ1×Ｅ（ｗy→x））を考慮しないと、各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）が時刻の変化に対して大きく振動することになる。これに対して、振動抑制項（λ1×Ｅ（ｗy→x））を加味して評価式ｃ5の値が最大となるよう各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）を演算することにより、各時刻ｔの重みｗy→x（ｔ）の時刻の変化に対する変化を平滑化させることができる。なお、実施例では、この各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）を演算する処理は、最急降下法（例えば、引用文献１参照）を用いて行なうものとした。また、式（２９）は、実施例では、各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）の総和が値Ｎになるという条件と、各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）が値０以上となる条件と、を共に満たすものとした。
引用文献１：馬場則夫・坂和正敏．数理計画法入門  ＰＣ９８０１パソコン学習支援システム，共立出版．
【００４４】
【数１６】
2008257595.tif 000018
【００４５】
こうして各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）を演算すると、Ｙ→Ｘ重み演算処理を終了するための演算終了条件が成立したか否かを判定し（ステップＳ２４０）、演算終了条件が成立していないときには、ステップＳ２１０に戻る。こうしてステップＳ２１０〜Ｓ２４０の処理を繰り返し実行することにより各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）を更新していき、演算終了条件が成立したときに、各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）をＲＡＭ２５の所定アドレスに書き込んで（ステップＳ２５０）、Ｙ→Ｘ重み演算処理を終了する。ここで、演算終了条件は、例えば、更新された各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）と更新前の対応する各時刻ｔのＹ→Ｘ重み（更新前ｗy→x（ｔ））との偏差の最大値が所定値以下になる条件や、ステップＳ２１０〜Ｓ２４０の処理を所定回数に亘って繰り返し実行した条件などとすることができる。
【００４６】
図２の因果関係判定ルーチンの説明に戻る。各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）を演算すると（ステップＳ１１０）、演算した各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）を用いてＹｄａｔａからＸｄａｔａへの各時刻ｔの因果関係の有無を判定し（ステップＳ１２０）、判定結果をユーザが指定した出力ファイルに出力したりディスプレイ４４に出力したりディスクドライブ装置４１にセットされたＣＤやＤＶＤに出力したりして（ステップＳ１３０）、因果関係判定プログラム３０を終了する。ステップＳ１２０の判定は、実施例では、ステップＳ１１０の処理により演算した各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）に対して、Ｙ→Ｘ重みｗy→x（ｔ）が値１以上のときにはＹｄａｔａからＸｄａｔａへの因果関係があると判定し、Ｙ→Ｘ重みｗy→x（ｔ）が値１より小さいときにはＹｄａｔａからＸｄａｔａへの因果関係はないと判定するものとした。このように、ＹｄａｔａからＸｄａｔａへの各時刻ｔの因果関係の有無を判定することにより、ＹｄａｔａからＸｄａｔａに各時刻ｔで因果関係が存在しているか否かを判定できるようにすることができる。
【００４７】
次に、図４に例示するＸ→Ｙ重み演算処理により時系列データＸｄａｔａからＹｄａｔａへの各時刻ｔの因果関係の程度に対する各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）を演算する（ステップＳ１４０）。Ｘ→Ｙ重み演算処理は、実施例では、式（３１）および式（３２）により表わされる回帰モデルと各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）とを用いて行なうものとした。ここで、式（３１）および式（３２）中、「β」，「Ｃ」，「Ｄ」は係数ベクトル，「θ」，「ｑ」は定数，「ε(y)（ｔ）」はＹｄａｔａの時刻ｔの入力値ｙ（ｔ）とＹｄａｔａの時刻ｔ−１の状態ベクトルＹt-1を用いて予測される時刻ｔの予測値との誤差（以下、Ｙ→Ｙ誤差という），「ε(y,x)（ｔ）」はＹｄａｔａの時刻ｔの入力値ｙ（ｔ）とＸｄａｔａおよびＹｄａｔａの時刻ｔ−１の状態ベクトルＸt-1，Ｙt-1を用いて予測される時刻ｔの予測値との誤差（以下、ＹＸ→Ｙ誤差という）である。状態ベクトルＸt-1，Ｙt-1は、前述の式（２３）および式（２４）により表わされる。以下、図２の因果関係判定ルーチンの説明を一旦中断し、図４のＸ→Ｙ重み演算処理について説明する。
【００４８】
【数１７】
2008257595.tif 000019
【００４９】
Ｘ→Ｙ演算処理では、まず、任意の時刻に対して各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）に値１を設定する（ステップＳ３００）。そして、式（３１）および式（３２）により表わされる回帰モデルと各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）とを用いた式（３３）および式（３４）により係数ベクトルβ，Ｃ，Ｄおよび定数θ，ｑを演算する（ステップＳ３１０）。ここで、式（３３）はＹ→Ｘ重み演算処理で用いられた前述の式（２７）と同様の構造の式であり、式（３３）中、Ｘ→Ｙ重みｗx→y（ｔ）とＹｄａｔａの入力値ｙ（ｔ）と状態ベクトルＹt-1は定まっている。このため、式（３３）では、Ｙ→Ｙ誤差ε(y)（ｔ）の期待値が最小となるよう係数ベクトルβおよび定数θを演算することになる。また、式（３４）は前述の式（２８）と同様の構造の式であり、式（３４）中、Ｘ→Ｙ重みｗx→y（ｔ）とＸｄａｔａ，Ｙｄａｔａの入力値ｘ（ｔ），ｙ（ｔ）と状態ベクトルＸt-1，Ｙt-1は定まっている。このため、式（３４）では、ＹＸ→Ｙ誤差ε(y,x)（ｔ）の期待値が最小となるよう係数ベクトルＣ，Ｄおよび定数ｑを演算することになる。
【００５０】
【数１８】
2008257595.tif 000020
【００５１】
こうして係数ベクトルβ，Ｃ，Ｄおよび定数θ，ｑを演算すると、演算した係数ベクトルβおよび定数θとＹｄａｔａの入力値ｙ（ｔ）とを用いて前述の式（３１）により各時刻ｔのＹ→Ｙ誤差ε(y)（ｔ）を演算すると共に係数ベクトルＣ，Ｄおよび定数ｑとＹｄａｔａの入力値ｙ（ｔ）とを用いて式（３２）により各時刻ｔのＹＸ→Ｙ誤差ε(y,x)（ｔ）を演算する（ステップＳ３２０）。
【００５２】
そして、各時刻ｔのＹ→Ｙ誤差ε(y)（ｔ）およびＹＸ→Ｙ誤差ε(y,x)（ｔ）を用いて式（３５）の評価式ｃ2の値が最大となるよう各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）を演算する（ステップＳ３３０）。ここで、式（３５）中、Ｙ→Ｙ誤差ε(y)（ｔ）やＹＸ→Ｙ誤差ε(y,x)（ｔ）を含む項は、各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）と、Ｙ→Ｙ誤差ε(y)（ｔ）の２乗とＹＸ→Ｙ誤差ε(y,x)（ｔ）の２乗との差と、の積の期待値である。また、式（３５）中、「λ2×Ｅ（ｗx→y）」はＸ→Ｙ重みｗx→y（ｔ）が時刻の変化に対して激しく振動しないようにするための振動抑制項である。そして、「Ｅ（ｗx→y）」は、式（３６）により表わされ、式（３６）中、「ｅ」はパラメータである。また、「λ2」はパラメータであり、λ1と同様にＲＡＭ２５の所定アドレスに書き込まれたものを読み込んで用いることができる。式（３５）は、Ｙ→Ｘ重み演算処理で用いられた前述の式（２９）と同様の構造の式であるから、式（２９）と同様に演算することができる。この処理では、振動抑制項（λ2×Ｅ（ｗx→y））を無視すれば、ＸｄａｔａからＹｄａｔａへの因果関係の程度が大きい時刻ｔほど大きなＸ→Ｙ重みｗx→y（ｔ）を演算し、ＸｄａｔａからＹｄａｔａへの因果関係の程度が小さい時刻ｔほど小さなＸ→Ｙ重みｗx→y（ｔ）を演算する処理となる。そして、振動抑制項（λ2×Ｅ（ｗx→y））を加味して評価式ｃ6の値が最大となるよう各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）を演算することにより、各時刻ｔの重みｗx→y（ｔ）の時刻変化に対する変化を平滑化させることができる。
【００５３】
【数１９】
2008257595.tif 000021
【００５４】
【数２０】
2008257595.tif 000022
【００５５】
こうして各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）を演算すると、Ｘ→Ｙ重み演算処理を終了するための演算終了条件が成立したか否かを判定し（ステップＳ３４０）、演算終了条件が成立していないときには、ステップＳ３１０に戻る。こうしてステップＳ３１０〜Ｓ３４０の処理を繰り返し実行することにより各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）を更新していき、演算終了条件が成立したときに、各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）をＲＡＭ２５の所定アドレスに書き込んで（ステップＳ３５０）、Ｘ→Ｙ重み演算処理を終了する。ここで、演算終了条件は、Ｙ→Ｘ重み演算処理と同様の条件を用いることができる。
【００５６】
図２の因果関係判定ルーチンの説明に戻る。各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）を演算すると（ステップＳ１４０）、演算した各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）を用いてＸｄａｔａからＹｄａｔａへの各時刻ｔの因果関係の有無を判定し（ステップＳ１５０）、判定結果をユーザが指定した出力ファイルに出力したりディスプレイ４４に出力したりディスクドライブ装置４１にセットされたＣＤやＤＶＤに出力したりして（ステップＳ１６０）、因果関係判定プログラム３０を終了する。ステップＳ１５０の判定は、実施例では、ステップＳ１４０の処理により演算した各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）に対して、Ｘ→Ｙ重みｗx→y（ｔ）が値１以上のときにはＸｄａｔａからＹｄａｔａへの因果関係があると判定し、Ｘ→Ｙ重みｗx→y（ｔ）が値１より小さいときにはＸｄａｔａからＹｄａｔａへの因果関係はないと判定するものとした。このように、ＸｄａｔａからＹｄａｔａへの各時刻ｔの因果関係の有無を判定することにより、ＸｄａｔａからＹｄａｔａに各時刻ｔで因果関係が存在しているか否かを判定できるようにすることができる。
【００５７】
次に、二つの異なる時系列データの一例として、生体に関する時系列データ、特に、心拍数と呼吸の時系列データを用いて両者間の各時刻ｔの因果関係の有無を判定する場合について説明する。実施例では、心拍数をＸｄａｔａ，呼吸をＹｄａｔａとして、Ｘｄａｔａ，Ｙｄａｔａについて、引用文献２に収録されている睡眠時無呼吸症候群の患者の睡眠時の１１７５秒後から１７７５秒後のデータを用いるものとした。なお、引用文献２のデータは、引用文献３に収録されているデータの一部から作成されたものである。また、実施例では、式（２９），（３５）におけるパラメータλ1，λ2についてはλ1＝λ2として値１．０，値５．０，値１０．０のいずれかを用いるものとし、式（３０），（３６）におけるパラメータｅについては値２を用いるものとし、式（２３），（２４）における埋め込み次元ｍ，ｎについてはｍ＝ｎとして値３を用いるものとし、同じく式（２３），（２４）における項の次数の最大値ｏ1，ｏ2についてはｏ1＝ｏ2として値２を用いるものとした。
引用文献２：D.R. Rigney, A.L. Goldberger, W. Ocasio, Y. Ichimaru, G.B. Moody, and R.Mark. "Multi-channel physiological data: description and analysis," In A.S. Weigend and N.A. Gershenfeld, editors, Time Sereies Prediction: Forecasting the Future and Understanding the Past, Proceedings of the Santa Fe Institute Studies in the Science of Complexity, Vol. XV. Addison-Wesley, Reading, MA, 1993.
引用文献３：Y. Ichimaru and G.B. Moody. "Development of the polysomnographic database on cd-rom, "Psychiatry and Clinical Neurosciences, Vol. 53, No. 2, pp.175-177, 1999.
【００５８】
図５〜図７は、心拍数Ｘｄａｔａと呼吸Ｙｄａｔａとの間の各時刻ｔの因果関係の有無を判定した結果の一例を示す。図５はパラメータλ1，λ2について共に値１．０を用いた場合の判定結果であり、図６はパラメータλ1，λ2について共に値５．０を用いた場合の判定結果であり、図７はパラメータλ1，λ2について共に値１０．０を用いた場合の判定結果である。図５〜図７では、上段から順に、（ａ）心拍数Ｘｄａｔａ，（ｂ）呼吸Ｙｄａｔａ，（ｃ）呼吸Ｙｄａｔａから心拍数Ｘｄａｔａへの因果関係の有無，（ｄ）心拍数Ｘｄａｔａから呼吸Ｙｄａｔａへの因果関係の有無を示している。図５〜図７において、呼吸Ｙｄａｔａから心拍数Ｘｄａｔａへの因果関係については、呼吸が行なわれている時刻（１１７５，１２２５，１２７０，１３１０，１３５５，１４００秒後頃）に因果関係関係があると判定された。これは、引用文献４に記載されている呼吸性の心拍変動を表わしていると考えられる。なお、呼吸Ｙｄａｔａから心拍数Ｘｄａｔａへの因果関係の有無について、実施例では、パラメータλ1を大きくするのに伴って振動が抑制された。また、図５〜図７において、心拍数Ｘｄａｔａから呼吸Ｙｄａｔａへの因果関係については、心拍数が急降下する時刻や心拍数が極大値に達する時刻（１１９０，１２４０，１３５５秒後頃）に因果関係があると判定された。この時刻後頃には、呼吸の停止と再開が行なわれているため、この結果は、心拍数が呼吸の停止と再開に関連していることを示唆していると考えることができる。なお、心拍数Ｘｄａｔａから呼吸Ｙｄａｔａへの因果関係については、パラメータλ2に拘わらず大きな振動は生じなかった。
引用文献４：F. Yasuma and J. Hayano, "Respiratory Sinus Arrhythmia", Chest, Vol. 125, pp. 683-690, 2004.
【００５９】
以上説明した実施例の因果関係判定装置２０によれば、ＲＡＭ２５に記憶された二つの異なる時系列データＸｄａｔａ，Ｙｄａｔａに対して、ＹｄａｔａからＸｄａｔａへの各時刻ｔの因果関係の程度とその因果関係に対する各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）とを用いた評価式の値が最大となるよう各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）を演算してこれをＲＡＭ２５に記憶すると共に記憶したＹ→Ｘ重みｗy→x（ｔ）を用いてＹｄａｔａからＸｄａｔａへの各時刻ｔの因果関係の有無を判定してこれを出力ファイルなどに出力し、ＸｄａｔａからＹｄａｔａへの各時刻ｔの因果関係の程度とその因果関係に対する各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）とを用いた評価式の値が最大となるよう各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）を演算してこれをＲＡＭ２５に記憶すると共に記憶したＸ→Ｙ重みｗx→y（ｔ）を用いてＸｄａｔａからＹｄａｔａへの各時刻ｔの因果関係の有無を判定してこれを出力ファイルなどに出力するから、ＹｄａｔａからＸｄａｔａに及びＸｄａｔａからＹｄａｔａに各時刻ｔで因果関係が存在しているか否かを判定することができる。
【００６０】
実施例の因果関係判定装置２０では、式（２１），（２２），（３１），（３２）により表わされる回帰モデルと、式（２３）および式（２４）により表わされる状態ベクトルＸt-1，Ｙt-1と、を用いるものとしたが、式（２３）および式（２４）に代えて、式（２５）および式（２６）により表わされる状態ベクトルＸt-1，Ｙt-1を用いるものとしてもよいし、これ以外の状態ベクトルを用いるものとしてもよい。また、式（２１），（２２），（３１），（３２）により表わされる回帰モデルにおいて、状態ベクトルＸt-1，Ｙt-1に代えて、状態ベクトルに所定の処理を施したものを用いるものとしてもよい。
【００６１】
実施例の因果関係判定装置２０では、Ｙ→Ｘ重み演算処理において、式（２９）により表わされる評価式ｃ5が最大となるようＹ→Ｘ重みｗy→x（ｔ）を演算し、式（３１）により表わされる評価式ｃ6が最大となるようＸ→Ｙ重みｗx→y（ｔ）を演算するものとしたが、式（２９），（３０）における振動抑制項を無視するものとしてもよい。この場合、Ｙ→Ｘ重み演算処理では、評価式ｃ5に代えて式（３７）により表わされる評価式ｃ7の値が最大となるようＹ→Ｘ重みｗy→x（ｔ）を演算すればよく、Ｘ→Ｙ重み演算処理では、評価式ｃ6に代えて式（３８）により表わされる評価式ｃ8の値が最大となるようＸ→Ｙ重みｗx→y（ｔ）を演算すればよい。このように振動抑制項を無視するときには、評価式ｃ7を用いて得られるＹ→Ｘ重みｗy→x（ｔ）や、評価式ｃ8を用いて得られるＸ→Ｙ重みｗx→y（ｔ）が時刻の変化に対して激しく振動することがあり得るため、Ｙ→Ｘ重みｗy→x（ｔ），Ｘ→Ｙ重みｗx→y（ｔ）に対してなまし処理を施したりローパスフィルタを施したりしたものをそれぞれＹ→Ｘ重みｗy→x（ｔ），Ｘ→Ｙ重みｗx→y（ｔ）として再設定するものとしてもよい。
【００６２】
【数２１】
2008257595.tif 000023
【００６３】
実施例の因果関係判定装置２０では、式（２９）により表わされる評価式ｃ5の値が最大となるよう各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）を演算するものとしたが、評価式ｃ5のうちＸ→Ｘ誤差ε(x)（ｔ）およびＸＹ→Ｘ誤差ε(x,y)（ｔ）に関する項ｄ1について、式（２９）に代えて、式（３９）により表わされるものを用いるものとしてもよい。また、実施例では、式（３０）により表わされる評価式ｃ6の値が最大となるよう各時刻ｔのＸ→Ｙ重みｗx→y（ｔ）を演算するものとしたが、評価式ｃ6のうちＹ→Ｙ誤差ε(y)（ｔ）およびＹＸ→Ｙ誤差ε(y,x)（ｔ）に関する項ｄ2について、式（３０）に代えて、式（４０）により表わされるものを用いるものとしてもよい。
【００６４】
【数２２】
2008257595.tif 000024
【００６５】
実施例の因果関係判定装置２０では、各時刻ｔのＹ→Ｘ重みｗy→x（ｔ）およびＸ→Ｙ重みｗx→y（ｔ）に対して、値１以上であれば因果関係が存在すると判定し、値１未満であれば因果関係は存在しないと判定するものとしたが、自動閾値設定法（例えば、引用文献５参照）などを用いて因果関係が存在するか否かを判定するものとしてもよい。
引用文献５：大津展之・栗田多喜夫・関田巖．行動計量学シリーズ１２巻  パターン認識−理論と応用−．朝倉書店，１９９６．
【００６６】
実施例では、二つの異なる時系列データの一例として、生体に関する時系列データ、特に、心拍数と呼吸の時系列データを用いて両者間の各時刻ｔの因果関係の有無を判定するものとしたが、二つの異なる時系列データとしては、心拍数と呼吸の時系列データの組み合わせ以外の生体に関する時系列データを用いるものとしてもよい。また、二つの異なる時系列データとして、生体に関する時系列データ以外の如何なる自然界の時系列データを用いるものとしてもよい。例えば、気象に関する時系列データ（例えば、一地点の気温と気圧の時系列データや二地点の気温と気温に時系列データ，二地点の気温と気圧の時系列データなど）を用いるものとしてもよいし、医療解析に関する時系列データ（例えば、薬剤の投与量と体温の時系列データや、薬剤の投与量と血中のある成分の濃度との時系列データなど）を用いるものとしてもよいし、通信ネットワークに関する時系列データ（例えば、サーバの負荷と通信速度の時系列データなど）を用いるものとしてもよいし、環境アセスメントに関する時系列データ（例えば、工事箇所の数と環境指標の値の時系列データなど）を用いるものとしてもよいし、神経科学的分析に関する時系列データ（例えば、脳の二点の神経信号データ（電圧など）など）を用いるものとしてもよいし、遺伝子ネットワーク解析に関する時系列データ（例えば、二つの異なる遺伝子の発現状態を示す時系列データなど）を用いるものとしてもよいし、行動認識に関する時系列データ（例えば、画像データと人の位置の時系列データや、二人の位置の時系列データなど）を用いるものとしてもよいし、自動車運転行動分析に関する時系列データ（例えば、車間距離とアクセル操作量の時系列データや、歩行者の位置とハンドル操作量の時系列データ，他車の位置と視線方向の時系列データなど）を用いるものとしてもよい。さらに、二つの異なる時系列データとしては、こうした自然界における時系列データに限定されるものではなく、自然界以外の時系列データを用いるものとしてもよい。
【００６７】
実施例では、コンピュータ２２に因果関係判定プログラム３０がインストールされたものとしての因果関係判定装置２０として説明したが、因果関係判定プログラム３０の形態として用いるものとしてもよい。
【００６８】
以上、本発明を実施するための最良の形態について実施例を用いて説明したが、本発明はこうした実施例に何等限定されるものではなく、本発明の要旨を逸脱しない範囲内において、種々なる形態で実施し得ることは勿論である。
【産業上の利用可能性】
【００６９】
本発明は、因果関係判定装置の製造産業や、コンピュータを因果関係判定装置として機能させるプログラムの製造産業などに利用可能である。
【図面の簡単な説明】
【００７０】
【図１】本発明の一実施例としての因果関係判定装置２０の構成の概略を示す構成図である。
【図２】実施例の因果関係判定ルーチンの一例を示すフローチャートである。
【図３】Ｙ→Ｘ重み演算処理の一例を示すフローチャートである。
【図４】Ｘ→Ｙ重み演算処理の一例を示すフローチャートである。
【図５】パラメータλ1，λ2について共に値１．０を用いた場合の心拍数Ｘｄａｔａと呼吸Ｙｄａｔａとの間の因果関係の有無を判定した結果を示す説明図である。
【図６】パラメータλ1，λ2について共に値５．０を用いた場合の心拍数Ｘｄａｔａと呼吸Ｙｄａｔａとの間の因果関係の有無を判定した結果を示す説明図である。
【図７】パラメータλ1，λ2について共に値１０．０を用いた場合の心拍数Ｘｄａｔａと呼吸Ｙｄａｔａとの間の因果関係の有無を判定した結果を示す説明図である。
【符号の説明】
【００７１】
２０  因果関係判定装置、２２  コンピュータ、２３  ＣＰＵ、２４  ＲＯＭ、２５  ＲＡＭ、２６  クロック回路、２７  ハードディスク、２８  入出力処理回路、２９  バス、３０  因果関係判定プログラム、３２  データ入力モジュール、３４  因果関係判定モジュール、３６  結果出力モジュール、４１  ディスクドライブ装置、４２  キーボード、４３  マウス、４４  ディスプレイ。
【図１】
2008257595.tif 000034
【図２】
2008257595.tif 000035
【図３】
2008257595.tif 000036
【図４】
2008257595.tif 000037
【図５】
2008257595.tif 000038
【図６】
2008257595.tif 000039
【図７】
2008257595.tif 000040