JP7483180B2 - 因果関係推定装置、因果関係推定方法、及び因果関係推定プログラム - Google Patents

Description

本開示は、因果関係推定装置、因果関係推定方法、及び因果関係推定プログラムに関する。

例えば、特許文献１は、２変量時系列データ（「２次元時系列データ」とも呼ばれる。）のカーネル平均の差が最大になる値であるｍａｘｉｍｕｍ ｍｅａｎ ｄｉｓｃｒｅｐａｎｃｙ（ＭＭＤ）を計算し、教師あり学習によって変量間の因果関係を明らかにする装置を開示している。

特開２０１７－２２８２５６号公報

Ｊａｃｋ Ｍ． Ｗａｎｇ， Ｄａｖｉｄ Ｊ． Ｆｌｅｅｔ， ａｎｄ Ａａｒｏｎ Ｈｅｒｔｚｍａｎｎ， "Ｇａｕｓｓｉａｎ Ｐｒｏｃｅｓｓ Ｄｙｎａｍｉｃａｌ Ｍｏｄｅｌｓ"， ＮＩＰＳ´０５： Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ ｔｈｅ １８ｔｈ Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｃｏｎｆｅｒｅｎｃｅ ｏｎ Ｎｅｕｒａｌ Ｉｎｆｏｒｍａｔｉｏｎ Ｐｒｏｃｅｓｓｉｎｇ Ｓｙｓｔｅｍｓ， Ｄｅｃｅｍｂｅｒ ２００５， ｐｐ．１４４１－１４４８ Ｅｄｗｉｎ Ｖ． Ｂｏｎｉｌｌａ， Ｋｉａｎ Ｍｉｎｇ Ａ． Ｃｈａｉ， Ｃｈｒｉｓｔｏｐｈｅｒ Ｋ．Ｉ． Ｗｉｌｌｉａｍｓ， "Ｍｕｌｔｉ－ｔａｓｋ Ｇａｕｓｓｉａｎ Ｐｒｏｃｅｓｓ Ｐｒｅｄｉｃｔｉｏｎ"， Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ ｔｈｅ Ａｄｖａｎｃｅｓ ｉｎ Ｎｅｕｒａｌ Ｉｎｆｏｒｍａｔｉｏｎ Ｐｒｏｃｅｓｓｉｎｇ Ｓｙｓｔｅｍｓ ２０， （２００８）

しかしながら、特許文献１の装置では、任意数の変量（例えば、３変量以上）の時系列データ間の因果関係を推定することができないという課題がある。

本開示は、上記課題を解決するためになされたものであり、任意数の変量の時系列データ間の因果関係の推定を可能にする因果関係推定装置、推定方法、及び推定プログラムを提供することを目的とする。

本開示の因果関係推定装置は、複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む学習用データを取得するデータ取得部と、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の因果関係を示す因果関係パラメータを計算し、前記学習用データと前記因果関係パラメータとからガウス過程の分散共分散行列を計算し、前記因果関係パラメータをマルチタスクガウス過程モデルで表現する計算部と、前記分散共分散行列に基づいて最適化関数を計算し、前記最適化関数に基づいて前記因果関係パラメータを更新する最適化部と、を有することを特徴とする。

本開示の他の因果関係推定装置は、因果関係パラメータデータベースから、複数の状態変数の時系列データと複数の観測変数の時系列データとの間の階層ごとの因果関係を示す因果関係パラメータを読み出し、前記因果関係パラメータに基づいて、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとを前記階層ごとに配置換えする因果グラフ構築部と、複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む検証用データを取得するデータ取得部と、配置換えされた前記複数の状態変数の時系列データと前記複数の観測変数の時系列データに対し、グレンジャー因果の前記検証用データを用いた検証と疑似相関の前記検証用データを行いた検証との一方又は両方を行う因果グラフ検証部とを有することを特徴とする。

本開示によれば、任意数の変量の時系列データ間の因果関係の推定を可能にすることができる。

実施の形態１に係る因果関係推定装置の構成を示すブロック図である。 図１の因果関係推定装置のハードウェア構成の例を示す図である。 図１のデータ取得部の前処理部の構成を示すブロック図である。 図１のデータ取得部の前処理部によって次元圧縮される前の時刻情報と次元圧縮された後の時刻情報との例を示す図である。 図１のデータ取得部の前処理部によって次元拡張される前の角度情報と次元拡張された後の角度情報との例を示す図である。 図１の学習部の計算部によって行われる処理の例を示す図である。 図１の学習部の計算部の時間シフトオペレータによって行われる処理の例を示す図である。 図１の学習部の計算部で作成された因果関係パラメータの初期値と個数の例を表形式で示す図である。 図１の因果関係推定装置の動作を示すフローチャートである。 実施の形態２に係る因果関係推定装置の構成を示すブロック図である。 図１０の因果関係推定装置のハードウェア構成の例を示す図である。 図１０の構築部の因果グラフ構築部による因果グラフの構築処理を示す図である。 （Ａ）及び（Ｂ）は、図１０の構築部の因果グラフ構築部によって行われる状態情報及び観測情報の配置替え処理の例を示す図である。 （Ａ）及び（Ｂ）は、図１０の構築部の因果グラフ検証部によって行われる検証処理の例を示す図である。 （Ａ）及び（Ｂ）は、図１０の構築部の因果グラフ検証部によって行われる検証処理の他の例を示す図である。 図１０の因果関係推定装置の動作を示すフローチャートである。 （Ａ）及び（Ｂ）は、図１０の因果関係推定装置によって時系列データ間の関係性から観測されていない時系列データを予測する処理を示す図である。 図１０の因果関係推定装置が、Ｃｈａｎｇｅ Ｐｏｉｎｔ Ｋｅｒｎｅｌを導入することでセンサデータの変化点又は故障を学習したモデルを使用する場合の動作を示す図である。

以下に、実施の形態に係る因果関係推定装置、因果関係推定方法、及び因果関係推定プログラムを、図面を参照しながら説明する。以下の実施の形態は、例にすぎず、実施の形態を適宜組み合わせること及び各実施の形態を適宜変更することが可能である。

実施の形態に係る因果関係推定装置は、例えば、各種経済指標（例えば、円ドルの為替レート、石油価格、公共投資、など）に関する時系列データ（Ｘ）と企業の株価の時系列データ（Ｙ）とからなる多次元時系列データがサンプルとして与えられた場合に、時系列データ（Ｘ）と時系列データ（Ｙ）との間の因果関係を、「Ｘ→Ｙ」すなわち「時系列データ（Ｘ）が原因であり、時系列データ（Ｙ）が結果である。」というように推定するための装置である。時系列データ（Ｘ）は状態変数と呼ばれ、時系列データ（Ｙ）は観測変数と呼ばれる。また、状態変数は、説明変数又は潜在変数とも呼ばれる。観測変数は、目的変数、従属変数、又は被説明変数とも呼ばれる。また、時系列データ（Ｘ）及び時系列データ（Ｙ）のいずれも、１次元時系列データである必要はなく、２変量以上の多次元時系列データ（すなわち、多変量時系列データ）であってもよい。また、状態変数である時系列データ（Ｘ）は、別の状態変数の時系列データ（Ｘ´）の観測変数であってもよい。また、観測変数である時系列データ（Ｙ）は、別の観測変数の時系列データ（Ｙ´）の状態変数であってもよい。

実施の形態１に係る因果関係推定装置は、学習用時系列データから学習済みモデル（因果関係パラメータを含む）を生成する機械学習装置である。実施の形態２に係る因果関係推定装置は、生成された学習済みモデルの因果関係パラメータ（例えば、実施の形態１に係る因果関係推定装置で生成された因果関係パラメータ）と検証用時系列データとから因果グラフを構築して出力する装置である。実施の形態２に係る因果関係推定装置は、実施の形態１に係る因果関係推定装置と別個の装置であるが、実施の形態１に係る因果関係推定装置の構成を有していてもよい。

実施の形態１に係る因果関係推定装置の学習部は、例えば、ガウス過程モデルを用いて、データ系列間の相関とデータ系列間のラグ情報とを含む多次元時系列データをモデリングする。ガウス過程モデルとしては、ガウス過程動的力学モデル（Ｇａｕｓｓ Ｐｒｏｃｅｓｓ Ｄｙｎａｍｉｃａｌ Ｍｏｄｅｌ：ＧＰＤＭ）及びマルチタスクガウス過程（Ｍｕｌｔｉ Ｔａｓｋ Ｇａｕｓｓ Ｐｒｏｃｅｓｓ：ＭＴＧＰ）モデル、などの公知の方法を用いることができる。ＧＰＤＭは、例えば、非特許文献１に記載されている。ＭＴＧＰモデルは、例えば、非特許文献２に記載されている。

ＧＰＤＭを用いることで、多次元時系列データの非線形な状態の時間発展を表現することができる。また、ＭＴＧＰモデルを用いることで、複数の状態関数と複数の観測関数との間の相関及び複数の状態関数と複数の観測関数との間のラグを、カーネル関数の因果関係パラメータとして表すことができる。実施の形態１に係る因果関係推定装置は、学習によって得られた因果関係を、記憶装置の因果関係パラメータデータベース（因果関係パラメータＤＢ）に因果関係パラメータとして保存する。実施の形態２に係る因果関係推定装置は、因果関係パラメータＤＢに保存されている因果関係パラメータから因果グラフを構築し、時系列データＤＢに保存されている検証用の時系列データを用いて因果グラフを検証して、検証済みの因果グラフを出力する。

《１》実施の形態１
《１－１》構成
《１－１－１》因果関係推定装置１００
図１は、実施の形態１に係る因果関係推定装置１００の構成を示すブロック図である。因果関係推定装置１００は、実施の形態１に係る因果関係推定方法（すなわち、学習方法）を実施することができる装置である。図１に示されるように、因果関係推定装置１００は、データ取得部１０と、学習部２０とを有してる。因果関係推定装置１００は、例えば、コンピュータである。学習によって得られた因果関係パラメータは、記憶部の因果関係パラメータＤＢ２６に記憶される。因果関係パラメータＤＢ２６は、因果関係推定装置１００の一部としての記憶部（後述の図２の記憶部１０３）に格納されてもよいが、因果関係推定装置１００とは別の外部装置（例えば、因果関係推定装置１００と通信可能なネットワーク上のサーバ）の記憶部に格納されてもよい。

因果関係推定装置１００は、複数の状態変数の時系列データの集合Ｘ（０：ｔ）と複数の観測変数の時系列データの集合Ｙ（０：ｔ）とを含む学習用データを取得し、複数の状態変数の時系列データと複数の観測変数の時系列データとの間の因果関係を示す因果関係パラメータθを計算し、学習用データと因果関係パラメータθとからガウス過程の分散共分散行列Ｋ（Ｘ，Ｘ´）を計算し、因果関係パラメータθをマルチタスクガウス過程モデルで表現し、分散共分散行列に基づいて最適化関数を計算し、最適化関数に基づいて因果関係パラメータθを更新する。因果関係パラメータθは、複数の状態変数の時系列データと複数の観測変数の時系列データとの間の相関を示す相関係数と、複数の状態変数の時系列データと複数の観測変数の時系列データとの間のラグを示すラグ係数とを含み、因果関係推定装置１００は、相関を、前記マルチタスクガウス過程モデルのＬＭＣ（Ｌｉｎｅａｒ Ｍｏｄｅｌ ｏｆ Ｃｏｒｅｇｉｏｎ）カーネルの線形相関で表現する。

データ取得部１０は、時系列データＤＢ１２から時系列データが入力され、前処理部１６に時系列データを出力する入力部１４と、入力部１４から出力された時系列データに前処理を施して学習部２０に前処理済みの時系列データを出力する前処理部１６とを有している。時系列データＤＢ１２は、因果関係推定装置１００の一部としての記憶部（後述の図２の記憶部１０３）に格納されてもよいが、因果関係推定装置１００とは別の外部装置（例えば、因果関係推定装置１００と通信可能なネットワーク上のサーバ）の記憶部に格納されてもよい。

学習部２０は、計算部２２と、最適化部２４とを有している。計算部２２は、因果関係パラメータを初期化し、分散共分散行列を計算する。最適化部２４は、分散共分散行列に基づいて最適化関数を計算し、最適化関数に基づいて因果関係パラメータを更新する。各部の詳細については後述する。

図２は、因果関係推定装置１００のハードウェア構成の例を示す図である。因果関係推定装置１００は、プロセッサ１０１と、メモリ１０２と、不揮発性の記憶装置である記憶部１０３とを有している。因果関係推定装置１００は、ネットワークを介して他の装置と通信を行う通信部を備えてもよい。プロセッサ１０１は、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）などである。メモリ１０２は、例えば、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）などの、揮発性の半導体メモリである。記憶部１０３は、ハードディスクドライブ（ＨＤＤ）又はソリッドステートドライブ（ＳＳＤ）などの記憶装置である。記憶部１０３は、情報（例えば、各種のデータベース）及びプログラムを記憶する。

因果関係推定装置１００の各機能は、処理回路により実現される。処理回路は、専用のハードウェアであってもよいし、メモリ１０２に格納されるプログラムを実行するプロセッサ１０１であってもよい。プロセッサ１０１は、処理装置、演算装置、マイクロプロセッサ、マイクロコンピュータ、及びＤＳＰ（Ｄｉｇｉｔａｌ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｏｒ）のいずれであってもよい。

処理回路が専用のハードウェアである場合、処理回路は、例えば、単一回路、複合回路、プログラム化したプロセッサ、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ－Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）、又はこれらのうちのいずれかを組み合わせたものである。

処理回路がプロセッサ１０１である場合、実施の形態１に係る因果関係推定プログラム（すなわち、学習プログラム）は、ソフトウェア、ファームウェア、又はソフトウェアとファームウェアとの組み合わせにより実現される。ソフトウェア及びファームウェアは、プログラムとして記述され、メモリ１０２に格納される。プロセッサ１０１は、メモリ１０２に記憶されたプログラムを読み出して実行することにより、図１に示される各部の機能を実現することができる。因果関係推定プログラムは、ネットワークを介してのダウンロードにより、又は、光ディスクなどのような情報を記録する記録媒体（すなわち、コンピュータで読取可能な記憶媒体）によって提供され、因果関係推定装置１００にインストールされる。なお、因果関係推定装置１００は、一部を専用のハードウェアで実現し、一部をソフトウェア又はファームウェアで実現するようにしてもよい。このように、処理回路は、ハードウェア、ソフトウェア、ファームウェア、又はこれらのうちのいずれかの組み合わせによって、図１に示される各機能ブロックの機能を実現することができる。

《１－１－２》入力部１４
データ取得部１０の入力部１４は、時系列データＤＢ１２から状態変数の時系列データの集合Ｘ（０：ｔ）（すなわち、複数個の時系列データ）及び観測変数の時系列データの集合Ｙ（０：ｔ）（すなわち、複数個の時系列データ）を取得し、これらを前処理部１６に出力する。

具体的に言えば、入力部１４は、時系列データＤＢ１２から、原因と結果の関係にあると考えられる、Ｑ個（Ｑは正の整数である。）の状態変数の時系列データｘ_１（０：ｔ），ｘ_２（０：ｔ），…，ｘ_Ｑ（０：ｔ）と、Ｄ個（Ｄは正の整数である。）の観測変数の時系列データｙ_１（０：ｔ），ｙ_２（０：ｔ），…，ｙ_Ｄ（０：ｔ）とを、それぞれ選択する。このことは、Ｑ個の状態変数の時系列データが原因となって、Ｄ個の観測変数の時系列データが観測されたことを意味する。なお、観測変数の時系列データの各々ｙ_ｄ（０：ｔ）（ｄ＝１，２，…，Ｄ）は、１次元時系列データである必要があるが、状態変数の時系列データの各々ｘ_ｑ（０：ｔ）（ｑ＝１，２，…，Ｑ）は、１次元時系列データである必要はなく、多次元時系列データであってもよい。例えば、状態変数の時系列データｘ_ｑ（０：ｔ）として状態情報である角度αの時系列データを選択し、後述するように状態情報である角度αを次元拡張してｓｉｎαの時系列データとｃｏｓαの時系列データとした場合、ｓｉｎαの時系列データとｃｏｓαの時系列データとのそれぞれを１次元時系列データの状態変数とせずに、これらのセットである（ｓｉｎα，ｃｏｓα）を１つの２次元時系列データに変換してもよい。

続いて、入力部１４では、時系列データＤＢ１２から選択したＱ個の状態変数の時系列データと、Ｄ個の観測変数の時系列データとから、長さＴ（Ｔは正の整数である。）の状態変数の時系列データの集合Ｘ（０：ｔ）と観測変数の時系列データの集合Ｙ（０：ｔ）を取得し、これらを前処理部１６に渡す。なお、Ｘ（０：ｔ）及びＹ（０：ｔ）において、括弧内の０：ｔは、時刻０から時刻ｔまでの時系列データを表す。ただし、時系列データは、０で始まる必要はなく、時系列データＤＢ１２に格納されている時系列データから任意の個所の長さＴの時系列データを選択してもよい。

Ｘ（０：ｔ）は、状態変数の時系列データの集合を表す。すなわち、
Ｘ（０：ｔ）＝｛ｘ_１（０：ｔ），ｘ_２（０：ｔ），…，ｘ_Ｑ（０：ｔ）｝である。
Ｙ（０：ｔ）は、観測変数の時系列データの集合を表す。すなわち、
Ｙ（０：ｔ）＝｛ｙ_１（０：ｔ），ｙ_２（０：ｔ），…，ｙ_Ｄ（０：ｔ）｝である。
なお、ｘ_ｑ（０：ｔ）は、状態ｑの状態変数の時系列データを表す。また、ｙ_ｄ（０：ｔ）は、ｄ次元目の観測変数の時系列データを表す。

《１－１－３》前処理部１６
データ取得部１０の前処理部１６は、入力部１４から出力された状態変数の時系列データの集合Ｘ（０：ｔ）と観測変数の時系列データの集合Ｙ（０：ｔ）とを取得し、これらに前処理を施し、前処理済みの状態変数の時系列データの集合Ｘ（０：ｔ）と前処理済みの観測変数の時系列データの集合Ｙ（０：ｔ）とを、学習部２０の計算部２２に出力する。

図３は、データ取得部１０の前処理部１６の構成を示すブロック図である。図３に示されるように、前処理部１６は、次元変更部１７と、正規化部１８とを有している。

次元変更部１７は、次元の変更が必要な状態情報の時系列データｘ_ｑ（０：ｔ）に対して、各々の状態情報に対して適切な次元削減処理（すなわち、次元圧縮処理）又は次元拡張処理を行う。状態情報の時系列データｘ_ｑ（０：ｔ）の次元の変更が必要な場合は、例えば、状態情報が時刻情報又は角度情報などの周期性のある情報である場合、及び長方形の各辺の長さよりも対角線の長さのほうが有効である場合、などである。以下に、次元圧縮処理と次元拡張処理の例を示す。なお、以下の例のそれぞれの逆の処理を行うことも可能である。

図４は、前処理部１６の次元変更部１７によって次元圧縮される前の時刻情報（左側の表）と次元圧縮された後の時刻情報（右側の表）との例を示す図である。図４に示されるように、時刻情報は、観測データが１日周期の場合には、２４時間ごとの周期性を表す。このため、２つの情報である「分」と「時」とを統合して１つの情報である「時分」に圧縮することで、時刻情報を３次元時系列データから２次元時系列データに圧縮することができる。

図５は、前処理部１６の次元変更部１７によって次元拡張される前の角度情報と次元拡張された後の角度情報との例を示す図である。図５に示されるように、角度情報は、１つの情報である「角度α°」から２つの情報の組合せである「ｓｉｎα」と「ｃｏｓα」に拡張することで、角度情報を１次元時系列データから２次元時系列データに拡張することができる。なお、次元変更部１７は、別の次元圧縮手法（例えば、主成分分析等）又は別の次元拡張方法を用いてもよい。

正規化部１８は、次元変更された状態変数の時系列データの集合Ｘ（０：ｔ）と次元変更された観測変数の時系列データの集合Ｙ（０：ｔ）とを、平均が０であり分散が１であるように正規化する。

《１－１－４》計算部２２
学習部２０の計算部２２は、前処理部１６から前処理済の状態変数の時系列データの集合Ｘ（０：ｔ）と前処理済の観測変数の時系列データの集合Ｙ（０：ｔ）を受け取り、ガウス過程の分散共分散行列Ｋ（Ｘ，Ｘ´）の計算を行い、最適化部２４に分散共分散行列Ｋ（Ｘ，Ｘ´）を出力する。

図６は、学習部２０の計算部２２によって行われる処理の例を示す図である。図６の例では、ＧＰＤＭで表される観測変数の時系列データ間の相関を、ＭＴＧＰモデルのＬＭＣカーネルの線形相関で表現する。

図６においては、ｔｉｍｅ－ｓｈｉｆｔ ｏｐｅｒａｔｏｒ（時間シフトオペレータ）（「Ｌａｇ Ｏｐｅｒａｔｏｒ（ラグオペレータ）」とも呼ばれる。）Ｆ^Ｌを導入することで、状態関数ｇ_１（ｘ_１（ｔ））～ｇ_Ｑ（ｘ_Ｑ（ｔ））が遅れて（又は、進んで）観測関数ｆ_１（ｘ_１（ｔ））～ｆ_Ｄ（ｘ_Ｄ（ｔ））に影響を及ぼすラグ効果を表現することができる。ラグ効果とは、原因が結果に同時刻に影響を及ぼすのではなく、「原因が遅れて結果に影響を及ぼす」又は「原因が速く結果に影響を及ぼす」ことを表す。ラグ効果の具体例は、原油価格が高くなって、ガソリン価格が連動して高くなるわけではなく、ある期間（例えば、次の仕入れが行われるまでの期間）を経てガソリン価格が高くなることである。あるいは、ラグ効果の他の具体例は、将来、原油価格が高くなることを見越して、先にガソリン価格を値上げすることである。

状態情報の一つに時刻情報を加えることで、計算部２２は、時系列性又は周期性を考慮することができる。ただし、計算部２２は、時刻情報に関してのラグ効果を考慮しないことも可能である。

以下、実施の形態１における分散共分散行列Ｋ（Ｘ，Ｘ´）の計算方法について説明する。ＧＰＤＭでは、状態方程式は式（１）のように定義され、観測方程式は式（２）のように定義される。

ここで、Ｘ（ｔ）は、状態変数の集合｛ｘ_１（ｔ），ｘ_２（ｔ），…，ｘ_Ｑ（ｔ）｝を表す。ｕ_ｑ（ｔ）及びｖ_ｄ（ｔ）は、ホワイトガウスノイズを表す。
ガウス過程状態空間モデルでは、モデルの状態遷移関数ｆ（ｘ）及び状態関数ｇ（ｘ）のガウス過程によるモデル化が行われ、ガウス過程の表記を用いて、一般に、以下の式（３）、（４）のように表される。
ｇ（ｘ）～ｇｐ（０，Ｋ_ｇ） （３）
ｆ（ｘ）～ｇｐ（０，Ｋ_ｆ） （４）
式（３）は、状態の非線形な時間発展を表し、式（４）は、状態関数から観測関数への変換を表す。また、Ｋｇ、Ｋｆは、グラム行列を表す。

式（１）の状態関数ｇ_ｑ（ｘ）は、状態情報ｘ_ｑから生成されるガウス過程のモデルｇｐ（０，ｋ（ｘ_ｑ，ｘ_ｑ´））を表す。

ｋ_ｑは、正定値カーネルを表す。使用する正定値カーネルとしては、ＲＢＦカーネル（Ｒａｄｉａｌ ｂａｓｉｓ ｆｕｎｃｔｉｏｎ ｋｅｒｎｅｌ）等のデータに適したものを選択する。例えば、ＲＢＦカーネルは、以下の式（５）で与えられる。

式（５）において、

は、ＲＢＦカーネルの因果関係パラメータである。

従来手法のＬＭＣの観測関数ｆ_ｄ（Ｘ（ｔ））は、式（６）で表される。

これに対し、実施の形態１で提案するＬＭＣの観測関数ｆ_ｄ（Ｘ（ｔ））は、以下の式（７）で表される。

従来手法の式（６）において、ａ_ｄ，ｑは、従来手法のＬＭＣの線形相関を表し、状態変数ｑから観測変数ｄへの相関係数を表す。

実施の形態１の式（７）において、

は、実施の形態１で新たに導入したｔｉｍｅ－ｓｈｉｆｔ ｏｐｅｒａｔｏｒであり、Ｌ_ｄ，ｑは、状態変数ｑから観測変数ｄへのラグ係数を表す。

Ｆ^Ｌは、時間軸に沿って状態関数の状態情報を未来（あるいは、過去）にずらすｔｉｍｅ－ｓｈｉｆｔ ｏｐｅｒａｔｏｒ（ラグオペレータ）を表し、以下の式（８）のように定義される。

図７は、学習部２０の計算部２２のｔｉｍｅ－ｓｈｉｆｔ ｏｐｅｒａｔｏｒによって行われる処理の例を示す図である。図７は、ｔｉｍｅ－ｓｈｉｆｔ ｏｐｅｒａｔｏｒ Ｆ^Ｌの動作を示す。図７の例は、情報を時間軸ｔに沿って未来にずらす場合（Ｌ＞０）を示している。情報を時間軸ｔに沿って過去にずらす場合（Ｌ＜０）は、ずらす方向は図７の場合の方向の逆方向である。

図７の例のように、状態情報のうち、時刻情報についてはｔｉｍｅ－ｓｈｉｆｔ ｏｐｅｒａｔｏｒを適用しないことも可能であり、適用しない場合にはＬ＝０を適用する。

図７の例では、計算部２２は、例えば、３次元時系列データである状態情報＃２にＬ＝１を適用しており、この場合、３次元時系列データのいずれも時間軸ｔに沿ってデータをシフトさせる。計算部２２は、シフト後、未来にはみ出した部分（すなわち、ｔ＝Ｔ＋１の部分）を削除し、空いた部分（ｔ＝０の部分）は、直前値（すなわち、ｔ＝１の部分の値）で穴埋めする。

図７の例では、計算部２２は、例えば、２次元時系列データである状態情報＃３にＬ＝２を適用しており、この場合、２次元時系列データのいずれも時間軸ｔに沿ってデータをシフトさせる。計算部２２は、シフト後、未来にはみ出した部分（すなわち、ｔ＝Ｔ＋１、ｔ＝Ｔ＋２の部分）を削除し、空いた部分（すなわち、ｔ＝０、ｔ＝１の部分）は直前値（ｔ＝２）で穴埋めする。

図７の例では、計算部２２は、例えば、３次元時系列データである状態情報＃ＱにＬ＝０を適用しており、この場合、３次元時系列データのシフトはなく、３次元時系列データは変更されない。

〈提案モデルの分散共分散行列の計算の仕方〉
既存手法により、式（７）の

の観測データがガウス過程に従うと、ｆは、以下の式（９）の多次元ガウス分布で表される。

式（９）において、

である。

Ｋ（Ｘ，Ｘ´）は、分散共分散行列（又は、グラム行列）と呼ばれ、状態変数間（すなわち、Ｘ⇔Ｘ´）の類似度を表す行列である。

式（９）において、μ（Ｘ）は、状態変数の平均行列を表す。正規化部１８で状態変数を正規化しているので、μ（Ｘ）＝０である。

分散共分散行列Ｋ（Ｘ，Ｘ´）の各成分（Ｋ（Ｘ，Ｘ´））_ｄ，ｄは、式（７）のｆ_ｄ（Ｘ）を用いて以下の式（１０）のように計算できる。

なお、式（１０）における２行目から３行目への式変換では、異なる状態間（ｑ≠ｑ´）は、独立であるため、

となる。

分散共分散行列Ｋ（Ｘ，Ｘ´）は、以下の式（１１）及び式（１２）のように表される。

ここで、Ｂ_ｑは、ｃｏｒｅｇｉｏｎａｌｉｚａｔｉｏｎ ｍａｔｒｉｘと呼ばれ、状態関数から観測関数への線形変換を表し、以下のように表記される。

分散共分散行列Ｋ（Ｘ，Ｘ´）の計算は、状態変数の時系列データの集合Ｘ（０：ｔ）と、以下に示される因果関係パラメータθを用いて行うことができる。

《１－１－５》最適化部２４
図８は、因果関係パラメータθの初期値と個数の例を表形式で示す図である。これらの因果関係パラメータは、学習部２０の最適化部２４で最適化される。

最適化部２４は、計算部２２によって計算された分散共分散行列Ｋ（Ｘ，Ｘ´）を受け取り、周辺尤度の計算と周辺尤度が最小化するように因果関係パラメータθを最適化する処理を行い、最適化した因果関係パラメータθを因果関係パラメータＤＢに保存する。

周辺尤度は、以下のように、計算することで得られる。観測情報が観測される確率は、以下のように式（１３）で得ることができる。

式（１３）の両辺の対数をとると、以下の周辺尤度を計算することができる。ただし、Ｋ_θ（Ｘ，Ｘ´）は、因果関係パラメータを用いて計算した分散共分散行列Ｋ（Ｘ，Ｘ´）を表す。

ただし、Ｎは、特徴量ベクトルＸの長さ、Ｄは、ｙの出力次元数である。

因果関係パラメータθの最適化には、周辺尤度ｌｏｇｐ（ｙ｜Ｘ，θ）が最大になるようにすればよい。一般的な最適化問題と形を合わせるため、式（１４）の両辺にマイナスをかけた以下の式（１５）で最適化関数Ｅを最小化するようにする。

最適化部２４は、式（１５）の最適化関数Ｅが最小化するようにθを更新する。最適化の際に、因果関係パラメータθを更新すると、Ｋ_θ（Ｘ，Ｘ´）の更新も必要になるため計算部２２でＫ_θ（Ｘ，Ｘ´）を計算する。

最適化部２４は、最適化に際し、既知技術である確率的勾配降下法などを用いることができる。例えば、Ｌ_ｄ，ｑは、グリッドサーチ用いて最適化することができ、残りの因果関係パラメータは、確率的勾配降下法を用いて最適化することができる。

以上の動作を実行することにより、最適化された多変量因果関係の因果関係パラメータθを求めることができる。従来は、状態変数（Ｑ＝１）、観測変数（Ｄ＝１）にした場合の２変量の時系列データから「状態情報→観測情報」の因果関係を推定していた。これに対し、実施の形態１では、状態変数及び観測変数の各々を１次元以上とすることが可能であるため、任意数の変量の時系列データ間の因果関係の推定が可能であり、多変量（Ｑ＋Ｄ）の時系列データの因果関係を推定することができる。つまり、実施の形態１では、複数の状態変数と複数の観測変数との間の因果関係がわかるようになり、２変量の時系列データの因果関係の推定だけでなく、３変量以上の時系列データの因果関係の推定を行うことができる。

《１－２》動作
図９は、実施の形態１に係る因果関係推定装置１００の動作（すなわち、学習方法）を示すフローチャートである。まず、ステップＳ１０１で、入力部１４は、時系列データＤＢ１２から状態変数の時系列データの集合Ｘ（０：ｔ）と観測変数の時系列データの集合Ｙ（０：ｔ）を取得する。

次に、ステップＳ１０２で、前処理部１６は、状態変数の時系列データの集合Ｘ（０：ｔ）において次元変更が必要な状態変数ｘ_ｑ（０：ｔ）の次元変更（すなわち、次元圧縮又は次元拡張）を行う。

次に、ステップＳ１０３で、前処理部１６は、状態変数の時系列データの集合Ｘ（０：ｔ）と観測変数の時系列データの集合Ｙ（０：ｔ）の正規化を行う。

次に、ステップＳ１０４で、計算部２２は、以下の因果関係パラメータθを初期値に設定する。

次に、ステップＳ１０５で、計算部２２は、状態変数の時系列データの集合Ｘ（０：ｔ）と因果関係パラメータθを用いて、式（１１）の分散共分散行列Ｋ_θ（Ｘ，Ｘ´）を計算する。

次に、ステップＳ１０６で、最適化部２４は、式（１５）の最適化関数Ｅを計算する。

次に、ステップＳ１０７で、最適化部２４は、最適化関数Ｅが最小になるように因果関係パラメータθを最適化（すなわち、更新）する。

分散共分散行列Ｋ_θ（Ｘ，Ｘ´）の更新の際には、ステップＳ１０５で、更新された因果関係パラメータθを用いて、式（１１）のＫ（Ｘ，Ｘ´）が計算される。

最適化関数Ｅが収束するまでステップＳ１０５～Ｓ１０７が繰り返される。最適化関数Ｅが収束したと判定されたときに、因果関係パラメータである相関係数ａ_ｄ，ｑ及びラグ係数Ｌ_ｄ，ｑが因果関係パラメータＤＢ２６に保存される。

《１－３》効果
実施の形態１によれば、複数の状態変数と複数の観測変数との間の因果関係がわかるようになり、任意数の変量の時系列データの因果関係を推定することができる。

《２》実施の形態２
《２－１》構成
《２－１－１》因果関係推定装置２００
因果関係推定装置２００は、因果関係パラメータＤＢから、複数の状態変数の時系列データと複数の観測変数の時系列データとの間の階層ごとの因果関係を示す因果関係パラメータθを読み出し、因果関係パラメータに基づいて、複数の状態変数の時系列データと複数の観測変数の時系列データとを階層ごとに配置換えし、複数の状態変数の時系列データの集合Ｘ（０：ｔ）と複数の観測変数の時系列データの集合Ｙ（０：ｔ）とを含む検証用データを取得し、配置換えされた複数の状態変数の時系列データと複数の観測変数の時系列データに対し、グレンジャー因果の検証用データを用いた検証と疑似相関の前記検証用データを行いた検証との一方又は両方を行う。

具体的に言えば、実施の形態２に係る因果関係推定装置２００は、因果グラフを構築する因果グラフ構築機能を備えた因果グラフ構築装置である。因果グラフとは、実施の形態１で得られた因果関係の情報に基づいて、「原因→結果」の関係でデータ項目一覧をグラフ構造にしたものである。構築された因果グラフの要素である状態関数及び観測関数は、因果関係パラメータＤＢに記憶された相関係数ａ_ｄ，ｑ及びラグ係数Ｌ_ｄ，ｑに基づいて、配置換えされる。構築された因果グラフは、グレンジャー因果と疑似相関とを用いて検証される。検証済みの因果グラフは、出力される。

実施の形態２に係る因果関係推定装置２００は、因果関係パラメータである相関係数ａ_ｄ，ｑ及びラグ係数Ｌ_ｄ，ｑを用いて因果グラフを構築する。因果関係パラメータである相関係数ａ_ｄ，ｑ及びラグ係数Ｌ_ｄ，ｑは、例えば、実施の形態１に係る因果関係推定装置１００で作成された学習済みモデルの因果関係パラメータである。因果関係推定装置２００は、例えば、コンピュータである。因果関係推定装置２００は、実施の形態１に係る因果関係推定装置１００を構成するコンピュータと同じコンピュータであってもよく、又は、異なるコンピュータであってもよい。

図１０は、実施の形態２に係る因果関係推定装置２００の構成を示すブロック図である。因果関係推定装置２００は、実施の形態２に係る因果関係推定（すなわち、因果グラフ構築方法）を実施することができる装置である。図１０に示されるように、因果関係推定装置２００は、構築部３０と、データ取得部４０と、出力部９０とを有している。構築部３０は、因果グラフ構築部３２と、因果グラフ検証部３４とを有している。構築部３０は、因果グラフ構築部３２に因果関係パラメータを提供する因果関係パラメータＤＢ８０を有してもよい。データ取得部４０は、入力部４４と、前処理部４６とを有している。データ取得部４０は、入力部４４に時系列データを提供する時系列データＤＢ４２を有してもよい。出力部９０は、構築部３０によって構築された因果グラフを出力する。時系列データＤＢ４２及び因果関係パラメータＤＢ８０は、因果関係推定装置２００の一部としての記憶部（後述の図１１の記憶部２０３）に格納されてもよいが、因果関係推定装置１００とは別の外部装置（例えば、因果関係推定装置２００と通信可能なネットワーク上のサーバ）の記憶部に格納されてもよい。

図１１は、因果関係推定装置２００のハードウェア構成の例を示す図である。因果関係推定装置２００は、プロセッサ２０１と、メモリ２０２と、不揮発性の記憶装置である記憶部２０３とを有している。因果関係推定装置２００は、外部の装置とのインタフェース、ネットワークを介して他の装置と通信を行う通信部、などを備えてもよい。プロセッサ２０１は、ＣＰＵなどである。メモリ２０２は、例えば、ＲＡＭなどの、揮発性の半導体メモリである。記憶部２０３は、ＨＤＤ又はＳＳＤなどの記憶装置である。記憶部２０３は、情報（例えば、各種のデータベース）及びプログラムを記憶する。

因果関係推定装置２００の各機能は、処理回路により実現される。処理回路は、専用のハードウェアであってもよいし、メモリ２０２に格納されるプログラムを実行するプロセッサ２０１であってもよい。

処理回路がプロセッサ２０１である場合、実施の形態２に係る推定プログラム（すなわち、因果グラフ構築プログラム）は、ソフトウェア、ファームウェア、又はソフトウェアとファームウェアとの組み合わせにより実現される。ソフトウェア及びファームウェアは、プログラムとして記述され、メモリ２０２に格納される。プロセッサ２０１は、メモリ２０２に記憶されたプログラムを読み出して実行することにより、図１０に示される各部の機能を実現することができる。プログラムは、ネットワークを介してのダウンロードにより、又は、光ディスクなどのような情報を記録する記録媒体から、因果関係推定装置２００にインストールされる。

《２－１－２》データ取得部４０
データ取得部４０は、実施の形態１のデータ取得部１０と同様の機能を有している。ただし、入力部４４は、後述するグレンジャー因果と疑似相関を用いた因果グラフの検証用のデータを取得する。

《２－１－３》因果グラフ構築部３２
構築部３０の因果グラフ構築部３２は、因果関係パラメータＤＢ８０に保存されている因果関係パラメータ

を取得し、因果関係パラメータを用いて因果グラフの構築を行い、構築された因果グラフを因果グラフ検証部３４に出力する。

図１２は、構築部３０の因果グラフ構築部３２による因果グラフの構築処理を示す図である。図１２において、ｙ（観測情報）よりもｘ（状態情報）がラグ（遅れ）があって、ｘ（状態情報）とｙ（観測情報）との間の相関係数が高い場合には、「ｘ´を原因としてｙ´が変動する」という因果関係（すなわち、ｘ´→ｙ´）がある可能性が高いと考えられる。

図１２おいて、因果関係パラメータは、以下のとおりである。

これらの因果関係パラメータは、階層ごとに「状態情報→観測情報」の順番に並べられる。なお、ｈは、階層番号を示す正の整数である。

図１３（Ａ）及び（Ｂ）は、構築部３０の因果グラフ構築部３２によって行われる状態情報及び観測情報の配置替え処理の例を示す図である。図１３（Ａ）及び（Ｂ）に示されるように、状態情報と観測情報との間は、状態情報から観測情報に向かう矢印で結ばれる。

図１３（Ａ）においてより、状態関数よりも観測情報のほうがラグ（進み）がある場合には、矢印の向きを逆にする。図１３（Ｂ）より、相関係数が低い場合には、交絡因子又は中間因子である可能性を考え、因果グラフの上位などに持っていくような配置換えを行う。

図１３（Ａ）及び（Ｂ）で配置換えを行った際には、配置換えを行った後の状態変数と観測変数を用いて実施の形態１に係る因果関係推定装置１００で因果関係パラメータを再度求める。

《２－１－４》因果グラフ検証部３４
図１４（Ａ）及び（Ｂ）は、構築部３０の因果グラフ検証部３４によって行われる検証処理の例を示す図である。グレンジャー因果は、例えば、図１４（Ａ）に示されるように、状態変数ｘ_２を含めて状態変数ｘ_１，ｘ_２，ｘ_３によって観測変数ｙ_１，ｙ_２を予測した場合（ケース１）に比べて、図１４（Ｂ）に示されるように、状態変数ｘ_２を削除して状態変数ｘ_１，ｘ_３によって観測変数ｙ_１，ｙ_２を予測した場合（ケース２）に予測精度が下がると、状態変数ｘ_２は観測変数ｙ_１，ｙ_２に対してグレンジャー因果があると考える。また、状態変数ｘ_１，ｘ_２，ｘ_３によって観測変数ｙ_１，ｙ_２を予測した場合（ケース１）に比べて、状態変数ｘ_１，ｘ_３によって観測変数ｙ_１，ｙ_２を予測した場合（ケース２）に予測精度が上がると、状態変数ｘ_２は観測変数ｙ_１，ｙ_２に対してグレンジャー因果がないと考える。グレンジャー因果がない場合には、状態変数ｘ_２の配置換え又は状態変数ｘ_２の因果グラフからの削除を行う。

予測には、因果関係推定装置１００で因果関係の推定に使用した状態変数の時系列データの集合Ｘ（０：ｔ）及び観測変数の時系列データの集合Ｙ（０：ｔ）から、未来の状態変数の時系列データの集合Ｘ（ｔ＋１：ｔ＋Δｔ）及び未来の観測変数の時系列データの集合Ｙ（ｔ＋１：ｔ＋Δｔ）を、予測の検証用データ（すなわち、テストデータ）として用いてもよい。

また、訓練データで因果関係パラメータ

を求めて、既知技術のガウス過程回帰で、以下の予測値

を求めてもよい。

予測誤差は、例えば、以下のＲＭＳＥ（Ｒｏｏｔ Ｍｅａｎ Ｓｑｕａｒｅｄ Ｅｒｒｏｒ）で評価することができる。

図１５（Ａ）及び（Ｂ）は、構築部３０の因果グラフ検証部３４によって行われる検証処理の他の例を示す図である。図１５（Ａ）及び（Ｂ）より、非定常状態（ある因子が外的要因などで急激な変化が起こる（インパルス応答））で、どのくらい因子に影響するか（又は、どれくらい変化がなく一定値をとるか）、を見積もることで検証することも可能である。

例えば、図１５（Ａ）に示される例では、状態変数ｘ_２を介して疑似相関がある観測変数ｙ_１とｙ_２を用いる方法があると考えられる。図１５（Ｂ）に示される例では、状態変数ｘ_２が観測変数ｙ_１とｙ_２の両方に影響を及ぼすが、状態変数ｘ_１とｘ_３は、観測変数ｙ_１とｙ_２のそれぞれ以外に影響を及ぼさない場合を考える。ｘ_２が固定値をとるなど変化がない場合には、ｘ_２が変化しないため観測変数ｙ_１とｙ_２に影響を及ばさなくなるため、観測変数ｙ_１とｙ_２の間には、相関がなくなると考えられる。

疑似相関を用いた因果関係の検証において、ｘ_２の影響がなくなっても観測変数ｙ_１とｙ_２の相関が残っている、あるいは、ｘ_１と観測変数ｙ_１、ｘ_２とｙ_２が高い相関を示されるようになる場合には、状態変数ｘ_２の配置換え又は因果グラフからの削除を行う。

因果グラフを構築することで、直接測定できない場所のセンサ値の予測又は異常原因の探索等を可能となる。適用事例として交通量又は気象情報からの路面損傷の予測又は、付近の雨量と水位データからの河川の水位予測、天候又は周辺地域の経済状況データを活用した電力需要予測などが考えられる。

《２－２》動作
図１６は、因果関係推定装置２００の動作（すなわち、推論動作）を示すフローチャートである。まず、ステップＳ２０１で、構築部３０は、因果関係パラメータＤＢから各階層の

を取得する。

次に、ステップＳ２０２で、因果グラフ構築部３２は、状態情報及び観測情報を、階層ごとに「状態情報→観測情報」の順に並べる。

次に、ステップＳ２０３で、因果グラフ構築部３２は、因果方向（すなわち、矢印の向き）を修正する。

次に、ステップＳ２０４で、因果グラフ構築部３２は、因果グラフの状態情報及び観測情報の配置換えを行う。

次に、ステップＳ２０５で、入力部４４は、時系列データＤＢ４２からグレンジャー因果の検証用データを取得する。

次に、ステップＳ２０６で、前処理部４６は、検証用データの次元変更を行う。

次に、ステップＳ２０７で、前処理部４６は、検証用データの正規化を行う。

次に、ステップＳ２０８で、因果グラフ検証部３４は、グレンジャー因果を用いて因果グラフを検証する。

次に、ステップＳ２０９で、入力部４４は、時系列データＤＢ４２から疑似相関による検証用のデータを取得する。

次に、ステップＳ２１０で、前処理部４６は、検証用データの次元変更を行う。

次に、ステップＳ２１１で、前処理部４６は、検証用データの正規化を行う。

次に、ステップＳ２１２で、因果グラフ検証部３４は、疑似相関を用いた因果関係の検証をする。

以上に示した処理に基づいて、因果グラフ検証部３４が検証した因果グラフは出力部９０に出力される。

《２－３》効果
実施の形態２によれば、因果グラフを構築することで、直接測定できない場所のセンサ値の予測又は異常原因の探索等が可能となる。

また、実施の形態２は、例えば、交通量又は気象情報からの路面損傷の予測、付近の雨量と水位データからの河川の水位予測、天候又は周辺地域の経済状況データを活用した電力需要予測、などに適用可能である。

《２－４》適用例
［観測値の予測］
図１７（Ａ）及び（Ｂ）は、因果関係推定装置２００によって時系列データ間の関係性から観測されていない時系列データを予測する処理を行う学習部と推論部とを示す。この場合、因果関係推定装置１００は、図１７（Ａ）に示される過去のＡ１駅、Ａ２駅、Ａ３駅の混雑情報をマルチタスク学習することで、学習済みモデルを生成する。因果関係推定装置１００は、Ａ１駅、Ａ２駅、Ａ３駅のうちの、Ａ１駅とＡ２駅の混雑情報のみがわかり、Ａ３駅の混雑情報がわからない場合（図１７（Ｂ）の場合）に、マルチタスク学習で生成された学習済みモデルを用いて、お互いの混雑度の相関又はラグ情報からＡ３駅の混雑度を推論部で予測（推論）することが可能である。

［異常検知又は故障予知］
図１８は、因果関係推定装置１００にＣｈａｎｇｅ Ｐｏｉｎｔ Ｋｅｒｎｅｌを導入することでセンサデータの変化点又は故障を学習した場合の動作を示す図である。センサデータの変化点又は故障の学習によって生成された学習済みモデルを用いることで、故障予知を行うことができる。

［Ｃａｕｓａｌ Ｉｍｐａｃｔへの適用］
公知の時系列因果の推論フレームワーク（例えば、Ｃａｕｓａｌ Ｉｍｐａｃｔ）において、標準の線形状態空間モデルを、実施の形態の方法（ＧＰＤＭ＋ＭＴＧＰモデル）で置き換えることも可能である。この場合、非線形又は非定常の相関も考慮できるようになり、予測精度を向上させることができる。

［計算の高速化］
上述したＭＴＧＰモデルの計算は、計算量Ｏ（Ｄ^３Ｎ^３）、メモリＯ（Ｄ^２Ｎ^２）と計算負荷又はメモリコストが高い。このため、例えば、乱択化フーリエ特徴（Ｒａｎｄｏｍ Ｆｏｕｒｉｅｒ ｆｅａｔｕｒｅｓ）又は変分化フーリエ特徴（Ｖａｒｉａｔｉｏｎａｌ Ｆｏｕｒｉｅｒ ｆｅａｔｕｒｅｓ）を用いることで、マルチタスク学習を高速化し、メモリコストを削減することが可能である。

１００ 因果関係推定装置、 ２００ 因果関係推定装置（因果グラフ構築装置）、 １０ データ取得部、 １２、４２ 時系列データＤＢ、 １４ 入力部、 １６ 前処理部、 １７ 次元変更部、 １８ 正規化部、 ２０ 学習部、 ２２ 計算部、 ２４ 最適化部、 ２６、８０ 因果関係パラメータＤＢ、 ３０ 構築部、 ３２ 因果グラフ構築部、 ３４ 因果グラフ検証部、 ４０ データ取得部、 ４４ 入力部、 ４６ 前処理部、 ９０ 出力部。

Claims (12)

1. 複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む学習用データを取得するデータ取得部と、
   前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の因果関係を示す因果関係パラメータを計算し、前記学習用データと前記因果関係パラメータとからガウス過程の分散共分散行列を計算し、前記因果関係パラメータをマルチタスクガウス過程モデルで表現する計算部と、
   前記分散共分散行列に基づいて最適化関数を計算し、前記最適化関数に基づいて前記因果関係パラメータを更新する最適化部と、
   を有することを特徴とする因果関係推定装置。
2. 前記因果関係パラメータは、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の相関を示す相関係数と、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間のラグを示すラグ係数とを含み、
   前記計算部は、前記相関を、前記マルチタスクガウス過程モデルのＬＭＣカーネルの線形相関で表現する
   ことを特徴とする請求項１に記載の因果関係推定装置。
3. 前記データ取得部は、前記複数の状態変数の時系列データの次元変更を行い、次元変更された前記状態変数の時系列データを前記計算部に提供する
   ことを特徴とする請求項１又は２に記載の因果関係推定装置。
4. 前記複数の状態変数の時系列データは、時刻情報の時系列データを含む
   ことを特徴とする請求項１又は２に記載の因果関係推定装置。
5. 前記複数の状態変数の時系列データは、角度情報の時系列データを含む
   ことを特徴とする請求項１又は２に記載の因果関係推定装置。
6. 前記最適化部は、更新された前記因果関係パラメータを因果関係パラメータデータベースに保存する
   ことを特徴とする請求項１又は２に記載の因果関係推定装置。
7. 前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の階層ごとの因果関係を示す前記因果関係パラメータに基づいて、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとを前記階層ごとに配置換えする因果グラフ構築部と、
   複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む検証用データを取得する他のデータ取得部と、
   配置換えされた前記複数の状態変数の時系列データと前記複数の観測変数の時系列データに対し、グレンジャー因果の前記検証用データを用いた検証と疑似相関の前記検証用データを行いた検証との一方又は両方を行う因果グラフ検証部と、
   を有することを特徴とする請求項１又は２に記載の因果関係推定装置。
8. 因果関係パラメータデータベースから、複数の状態変数の時系列データと複数の観測変数の時系列データとの間の階層ごとの因果関係を示す因果関係パラメータを読み出し、前記因果関係パラメータに基づいて、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとを前記階層ごとに配置換えする因果グラフ構築部と、
   複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む検証用データを取得するデータ取得部と、
   配置換えされた前記複数の状態変数の時系列データと前記複数の観測変数の時系列データに対し、グレンジャー因果の前記検証用データを用いた検証と疑似相関の前記検証用データを行いた検証との一方又は両方を行う因果グラフ検証部と、
   を有することを特徴とする因果関係推定装置。
9. 前記因果関係パラメータは、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の相関を示す相関係数と、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間のラグを示すラグ係数とを含む
   ことを特徴とする請求項８に記載の因果関係推定装置。
10. 予測対象に関する複数の状態変数の時系列データと複数の観測変数の時系列データとを取得し、前記予測対象に関する前記複数の状態変数の時系列データと前記複数の観測変数の時系列データから未観測の次元の観測情報を予測するための、前記因果関係パラメータに基づく学習済みモデルを用いて、前記未観測の次元の観測情報を予測する推論部をさらに有する
    ことを特徴とする請求項１、２、８、９のいずれか１項に記載の因果関係推定装置。
11. 因果関係推定装置によって実施される因果関係推定方法であって、
    複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む学習用データを取得するステップと、
    前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の因果関係を示す因果関係パラメータを計算し、前記学習用データと前記因果関係パラメータとからガウス過程の分散共分散行列を計算し、前記因果関係パラメータをマルチタスクガウス過程モデルで表現するステップと、
    前記分散共分散行列に基づいて最適化関数を計算し、前記最適化関数に基づいて前記因果関係パラメータを更新するステップと、
    を有することを特徴とする因果関係推定方法。
12. 複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む学習用データを取得するステップと、
    前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の因果関係を示す因果関係パラメータを計算し、前記学習用データと前記因果関係パラメータとからガウス過程の分散共分散行列を計算し、前記因果関係パラメータをマルチタスクガウス過程モデルで表現するステップと、
    前記分散共分散行列に基づいて最適化関数を計算し、前記最適化関数に基づいて前記因果関係パラメータを更新するステップと、
    をコンピュータに実行させることを特徴とする因果関係推定プログラム。
    

Images