(57)【要約】【課題】  与えられた多次元の任意のデータのなかから、異常値といったような所定のデータを適切に可視化して検出する。【解決手段】  異常値可視化装置Ａを用いて対象データＹを構成する複数のデータｙのなかの所定のデータを、Ｋ個の自己回帰モデルを混合した自己回帰混合モデルを用いて可視化する。異常値可視化装置Ａは、モデルパラメータ推定部１１で、対象データＹを読み込んでメモリに記憶し、データｙとこのデータｙに基づいて作成される入力ベクトルｘとが対になった複数のデータＤを作成してメモリに記憶し、拡大確率ベクトル推定部１２で、このように記憶したデータＤとクラスｋとの同時確率を推定してメモリに記憶し、確率ベクトルに平均２シグマ値を結合することで拡大確率ベクトルを作成してメモリに記憶し、異常値可視化部１３で、この拡大確率ベクトルを用いて所定のデータを表示して可視化する。【選択図】    図１
2005258599.tif 000002
【特許請求の範囲】
【請求項１】
情報を記憶するメモリを作業領域として演算を行う演算手段を有するコンピュータを用いて、対象データＹを構成する複数のデータｙのなかの所定のデータを、Ｋ個の確率モデルの線形和で記述される混合モデルを用いて可視化するデータの可視化方法において、
前記演算手段が、
前記対象データＹを読み込んで前記メモリに記憶し、
前記データｙと各確率モデルとの同時確率をベクトルにしたものを確率ベクトルとして設定して前記メモリに記憶し、
前記確率ベクトルに、全確率モデルの所定の確率変数における確率密度の平均である平均確率密度を結合することで拡大確率ベクトルを作成して前記メモリに記憶し、
前記拡大確率ベクトルを用いて前記所定のデータを可視化する座標データを作成して前記メモリに記憶すること、
を特徴とするデータの可視化方法。
【請求項２】
情報を記憶するメモリを作業領域として演算を行う演算手段を有するコンピュータを用いて、対象データＹを構成する複数のデータｙのなかの所定のデータを、パラメータが異なるＫ個の自己回帰モデルの線形和で記述される自己回帰混合モデルを用いて可視化するデータの可視化方法において、
前記演算手段が、
前記対象データＹを読み込んで前記メモリに記憶する手順、
前記対象データＹに含まれる複数のデータｙと前記自己回帰モデルの次数τに基づいて、前記自己回帰モデルへの入力ベクトルｘを複数作成して前記メモリに記憶する手順、
前記複数のデータｙと複数の入力ベクトルｘとを用いて、前記自己回帰モデルのパラメータを算出して前記メモリに記憶する手順、
自己回帰混合モデルのパラメータを、前記複数のデータｙと前記複数の入力ベクトルｘとを用いて算出する手順を、前記自己回帰モデルの数だけ順次繰り返して前記メモリに記憶する手順、
前記入力ベクトルｘと当該入力ベクトルｘに対応する前記データｙとからなるデータＤと各自己回帰モデルとの同時確率をベクトルにしたものを確率ベクトルとして設定して前記メモリに記憶する手順、
全自己回帰モデルの所定の確率変数における確率密度の平均である平均確率密度を算出して前記メモリに記憶する手順、
前記確率ベクトルに前記平均確率密度を付加することにより拡大確率ベクトルを推定して前記メモリに記憶する手順、
前記拡大確率ベクトルを用いて視覚的に表示するための座標を算出して前記メモリに記憶する手順、
を備え、
このように算出した座標を用いて前記所定のデータを可視化することを特徴とするデータの可視化方法。
【請求項３】
前記パラメータの算出が、
ＥＭアルゴリズムのＥ−ｓｔｅｐとＭ−ｓｔｅｐを繰り返して実行し、各算出されるパラメータに基づいてデータの尤度の対数をとった判定値を算出し、この判定値を最大化する前記パラメータを検証されたパラメータとして採用する手順であること、
を特徴とする請求項２に記載のデータの可視化方法。
【請求項４】
前記拡大確率ベクトルを用いて可視化する手順として、
前記複数のデータＤのうちの１つのデータＤと、このデータＤとは異なる他の１つのデータＤとの類似度を、前記拡大確率ベクトルに基づいて算出して前記メモリに記憶する手順、
この類似度に基づいて前記複数のデータの座標を算出して前記メモリに記憶する手順、
を行う、ＣｏＰＥ法を用いること、
を特徴とする請求項２又は請求項３に記載のデータの可視化方法。
【請求項５】
前記対象データＹが、多次元の時系列のデータであることを特徴とする請求項１ないし請求項４のいずれか１項に記載のデータの可視化方法。
【請求項６】
前記所定の確率変数における確率密度が２シグマ値であり、平均確率密度が２シグマ値であること、
を特徴とする請求項１ないし請求項５のいずれか１項に記載のデータの可視化方法。
【請求項７】
情報を記憶するメモリを作業領域として演算を行う演算手段と前記演算手段での演算結果を表示する表示装置を有して、対象データＹを構成する複数のデータｙのなかの所定のデータを、パラメータが異なるＫ個の自己回帰モデルの線形和で記述される自己回帰混合モデルを用いて可視化するデータの可視化装置において、
前記演算手段は、
前記対象データＹに含まれる複数のデータｙと前記自己回帰モデルの次数τに基づいて、前記自己回帰モデルへの入力ベクトルｘを複数作成する機能、
前記複数のデータｙと複数の入力ベクトルｘを用いて、前記自己回帰モデルのパラメータを算出する機能、
自己回帰混合モデルのパラメータを、前記複数のデータｙと前記複数の入力ベクトルｘとを用いて算出する手順を、前記自己回帰モデルの数だけ順次繰り返して行う機能、
を備えたモデルパラメータ推定部と、
前記入力ベクトルｘを入力として対応する前記データｙを出力とする複数のデータＤと各自己回帰モデルとの同時確率をベクトルにしたものを確率ベクトルとして設定する機能、
全自己回帰モデルの所定の確率変数における確率密度の平均である平均確率密度を算出する機能、
前記確率ベクトルに前記平均確率密度を付加することにより拡大確率ベクトルを推定する機能、
を備えた拡大確率ベクトル推定部と、
前記複数のデータＤのうちの１つのデータＤと、このデータＤとは異なる他の１つのデータＤとの類似度を、前記拡大確率ベクトルに基づいて算出する機能、
この類似度に基づいて、前記複数のデータの座標を算出する機能、
この算出した座標を前記表示装置に表示するように指示する機能、
を備えた異常値可視化部と、
を含んで構成されることを特徴とするデータの可視化装置。
【請求項８】
請求項１ないし請求項６のいずれか１項に記載のデータの可視化方法を、情報を記憶するメモリを作業領域として演算を行う演算手段を有するコンピュータに実行させることを特徴とするデータの可視化プログラム。
【請求項９】
請求項１ないし請求項６のいずれか１項に記載のデータの可視化方法を、情報を記憶するメモリを作業領域として演算を行う演算手段を有するコンピュータに実行させる可視化プログラムを記憶したことを特徴とする記憶媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、経済時系列データ、テキストデータ、生物データなど、任意のデータが与えられたときに、例えば該データに含まれる所定のデータを可視化する方法や装置などに関する。
【背景技術】
【０００２】
近年、大量のデータが電子的に蓄積されつつある。このように電子的に蓄積されるデータから異常値を検出することによって、過去の大きな出来事を知ること、また将来の危機に備えることが、重要な研究課題となっている。
異常値検出のための主要な方法として、正常な状態のデータを利用してモデルを構築し、該モデルとの誤差の大きいデータを異常とするものがある（非特許文献１）。また、異常値を出力する方法としては、横軸にデータ、縦軸に異常の程度をプロットするものがある。
【非特許文献１】Early statistical detection of anthrax outbreaks by tracking over-the-counter medication sales, A. Goldenberg et.al. Proceedings of the National Academy of Sciences of the United States of America, 99(8):5237-5240,2002
【発明の開示】
【発明が解決しようとする課題】
【０００３】
しかしながら、正常な状態でのモデルを構築し、該構築したモデルとの誤差の大きいデータを異常とすると、モデルの表現能力の欠如のために、異常値ではないのに異常値とみなされてしまうことがある。正確に異常値を検出するために、モデルをより表現能力の高いものにする必要がある。
【０００４】
さらに、従来の異常値の出力法では、異常の程度を知ることはできるが、各異常の性質
を知ることは困難であるという問題点がある。
【０００５】
そこで、本発明は、与えられた任意のデータのなかから、例えば異常値といったような所定のデータを適切に可視化して検出することができるようにした、データの可視化方法、データの可視化装置、データの可視化プログラム、及び記憶媒体を提供することを主たる目的とする。
【課題を解決するための手段】
【０００６】
前記課題を解決した本発明は、与えられた任意のデータに基づき確率モデルに含まれるパラメータが異なる複数の確率モデルの線形和で記述される混合モデルを構築し、データと各確率モデルとの同時確率と平均２シグマ値を結合したものを表示することにより、異常値といったような所定のデータを可視化することとした。
【０００７】
即ち、本発明は、情報を記憶するメモリを作業領域として演算を行う演算手段を有するコンピュータを用いて、対象データＹを構成する複数のデータｙのなかの所定のデータを、Ｋ個の確率モデルの線形和で記述される混合モデルを用いて可視化するデータの可視化方法に関する。本発明のデータの可視化方法では、前記演算手段が、前記対象データＹを読み込んで前記メモリに記憶し、前記データｙと各確率モデルとの同時確率をベクトルにしたものを確率ベクトルとして設定して前記メモリに記憶し、前記確率ベクトルに、全確率モデルの所定の確率変数における確率密度である平均確率密度を結合することで拡大確率ベクトルを作成して前記メモリに記憶し、前記拡大確率ベクトルを用いて前記所定のデータを可視化する座標データを作成して前記メモリに記憶することを特徴とする。
【０００８】
ここで、確率モデルは、後記する実施形態では、自己回帰モデルであるが、正規分布モデルや多項分布モデルなどを適用することができる。また、混合モデルは、後記する実施形態では、自己回帰混合モデルであるが、確率モデルに対応した他の混合モデルでもよい。また、平均確率密度は、後記する実施形態では、平均２シグマ値であるが、他の平均確率密度でもよい。なお、詳しい解決手段については、後記する実施形態で詳しく説明する。
【発明の効果】
【０００９】
本発明によれば、与えられた任意のデータのなかから、例えば異常値といったような所定のデータを適切に可視化して検出することができる。
【発明を実施するための最良の形態】
【００１０】
以下、本発明の「データの可視化方法、データの可視化装置、データの可視化プログラム、及び記憶媒体」を実施するための最良の形態（以下「実施形態」という）を詳細に説明する。なお、以下説明する実施形態は、最初に可視化方法（データの可視化方法）の原理（「可視化方法の基本的な考え方」・「可視化方法の原理・数式の説明」）を説明し、その後、可視化方法などを具現化した「異常値可視化装置」を具体的に説明する。
【００１１】
≪可視化方法の基本的な考え方≫
最初に、可視化方法の基本的な考え方を、式（１）〜式（５）を参照して説明する。
【００１２】
（拡大確率ベクトルの推定）
本実施形態では、可視化するデータは、次の式（１）の構造を有するＤ*とする。
なお、この式（１）におけるＤ（ｉ）はｉ番目のデータ、Ｎは与えられたデータの総数である。このＤ（ｉ）は、後記するデータｙ（ｔ）と後記する入力ベクトルｘ（ｔ）が対になったものである。補足すると、Ｄ（ｉ）は、「入力ベクトルｘ（ｉ）とデータｙ（ｉ）とからなるデータ」である。
ちなみに、データＤ（ｉ）やデータｙ（ｔ）におけるカッコ内のｉやｔは、配列変数として定義されているデータを指定するインデックスである。インデックスの名称が異なっても、インデックスの値が同じならば同じデータを指定することができる。
【００１３】
【数１】
2005258599.tif 000003
【００１４】
データＤ*（つまりデータＤ（ｉ））は、Ｋ個の確率モデルの重み付き線形和によって記述される次の式（２）の混合モデルによって生成される。つまり、混合モデルは、式（２）に示すように、複数の確率モデルの線形和で定義される。なお、後記するように、Ｋは確率モデルの総数である。
【００１５】
【数２】
2005258599.tif 000004
【００１６】
なお、ΣＰ（ｋ）＝１、Ｐ（ｋ）≧０、Θは未知パラメータである。ここで、未知パラメータのΘは後記するＡとσの少なくとも１つに相当する。
未知パラメータは、未知パラメータのΘが与えられたもとでのデータＤ*が生成される確率を最大化することによって推定する。つまり、次の式（３）を最大化することで未知パラメータのΘが推定される。ここで、式（３）の左辺はΘが推定値であることを示す。また、式（３）の右辺は、右辺を最大とするΘを推定値とすることを意味する。ちなみに、式（３）は後記する式（１１）〜式（１７）に対応する。
【００１７】
【数３】
2005258599.tif 000005
【００１８】
データＤ（ｉ）と各確率モデル（ｋ＝１， ... ，Ｋ）との同時確率をベクトルにしたものを、次の式（４）で示す確率ベクトルｚ*（ｉ）とする。なお、ここでもカッコ内のｉは配列変数である確率ベクトルの任意の１つを指定するインデックスである。
ここで、式（４）は、「データｙと確率モデルとの同時確率をベクトルにしたものを確率ベクトルとして設定」したものである。
【００１９】
【数４】
2005258599.tif 000006
【００２０】
この確率ベクトルｚ*（ｉ）に平均２シグマ値（α〜）を結合（付加）したものを、次の式（５）の拡大確率ベクトルｚ（ｉ）とする。
【００２１】
【数５】
2005258599.tif 000007
【００２２】
平均２シグマ値とは、全確率モデルの２シグマ値（標準偏差の２倍）における確率密度の平均である。ちなみに、式（４）は後記する式（１９）に対応し、式（５）は後記する式（１８）に対応する。
【００２３】
（拡大確率ベクトルの可視化による異常値の可視化）
次に、拡大確率ベクトルを可視化することにより異常値を可視化する手法について説明する。拡大確率ベクトルを、例えば２次元や３次元のベクトルに次元圧縮して可視化する。次元圧縮法としては古典的多次元尺度法、ＣｏＰＥ法など数多く提案されており、どの手法を用いても構わない。後記するのはＣｏＰＥ法である。
【００２４】
拡大確率ベクトルのなかで、平均２シグマ値が相対的に高いデータは、異常値の程度が大きい（他のデータとの異なり方が大きい）ものであるということができる。そのため拡大確率ベクトルを可視化することにより、異常値の程度を表す軸ができ、異常値を視覚的に検出することが可能となる。なお、確率ベクトルは各確率モデルに属する確率を表しているため、拡大確率ベクトルを可視化することによって、異常値の程度の大きさだけでなく、異常値がどのような性質を持つかも知ることができる。
補足すると、表示されるデータが複数の確率モデルのうち、どの確率モデルにより生成されたかが検証可能であるので、例えば、異常値と判定されたデータがどの確率モデルにより生成されたかを知り、その確率モデルのパラメータを確認すれば、異常値がどのような性質を持つかを知ることができる。
ちなみに、表示するデータを、例えば確率モデルごとに色分けするなど、一の確率モデルにより処理されたデータと、パラメータが異なる他の確率モデルにより処理されたデータとを、識別可能に表示することで、異常値と判定されたデータがどの確率モデルにより生成されたかを視覚的に知ることができる。
【００２５】
≪可視化方法の原理・数式の説明≫
次に、与えられた任意のデータに含まれる異常値を、混合モデルを用いて可視化する可視化方法の原理と使用する数式を説明する。なお、数式の説明順序は、後で説明するフローチャート（図３など）での処理の説明順序とは必ずしも一致していない。
【００２６】
（データ）
まず、与えられた任意のデータ（つまり対象データＹ）を説明する。本実施形態の可視化方法の適用対象となる対象データＹは、次の式（６）の構造をしている。
例えば、ｙ（１）は１９８３年１月の経済時系列データ、ｙ（２）は翌月の経済時系列データ、ｙ（Ｔ）は１９８３年１月から数えてＴ番目の月の経済時系列データである。仮にＴが２４０ならば、ｙ（Ｔ）は２００２年１２月の経済時系列データである。つまり、データｙ（ｔ）は時刻ｔの経済時系列データ、Ｔは与えられた時系列の総数である。ちなみに、この経済時系列データは、１９８３年１月から２００２年１２までの、「日本のマネタリーベース」、「国際金利」、「卸売物価指数」、「機械受注」、「鉱工業生産指数」、「円ドル為替レート」の月ごとの、６次元の時系列データであるものとする（時刻ｔ＝１９８３年１月〜２００２年１２月、変数の次元数ｄ＝６、時系列の総数Ｔ＝２４０）。なお、式（６）の対象データＹは、「複数のデータｙからなる任意の対象データＹ」に相当する。ここでは、対象データＹとして６次元のデータを示しているが、次元は、１次元以上でよい。
【００２７】
【数６】
2005258599.tif 000008
【００２８】
（自己回帰モデルへの入力ベクトル）
対象データＹに含まれる異常値を可視化するに際して、対象データＹから式（７）で示す構造をした自己回帰モデルへの入力ベクトルｘ（ｔ）を作成する。なお、ｘ（ｔ）を適宜「入力ベクトル」と略称する。ちなみに、可視化の対象となる対象データＹに含まれる「異常値」は、「所定のデータ」に相当する。
【００２９】
【数７】
2005258599.tif 000009
【００３０】
この式（７）で、τは自己回帰モデルの次数を示す。例えば、次数τが２ならば時刻ｔから１月遡った時刻のデータｙ（ｔ−１）と２月遡った時刻のデータｙ（ｔ−２）に基づいて、時刻ｔにおける入力ベクトルｘ（ｔ）が作成されることになる。なお、式（７）における「´」は行列の転置である。ちなみに、入力ベクトルｘ（ｔ）は、ｄτのベクトルを有するデータになる。ここで、ｄは前記のとおり変数の次元数（ここでは６）である。例えば、次数τが２の場合は、ｄτは１２になる。
【００３１】
（自己回帰モデル）
自己回帰モデルは、次の式（８）で与えられる。すなわち、自己回帰モデルでは、時刻ｔのデータｙ（ｔ）を過去τ期間のデータの重み付き線形和と誤差との和で記述する。つまり、時刻ｔにおけるデータｙ（ｔ）は、次の式（８）で表現される。
なお、式（８）のＡは、ｄτ×ｄ行列の自己回帰モデルのパラメータ（適宜「パラメータ」という）、εは白色雑音の誤差項である。
【００３２】
【数８】
2005258599.tif 000010
【００３３】
（自己回帰混合モデル）
自己回帰混合モデルでは、次の式（９）のように、時刻ｔのデータｙ（ｔ）は、Ｋ個の自己回帰モデルの重み付き線形和によって生成されると考える（Ｋは１以上）。なお、ｋは、自己回帰モデルを特定（指定）するインデックスであり、１，２，…，Ｋの値をとる）。
【００３４】
【数９】
2005258599.tif 000011
【００３５】
ここで、式（９）の左辺は、入力ベクトルｘ（ｔ）が解っているときのデータｙ（ｔ）の確率密度である。また、式（９）の右辺のＰ（ｋ）は、自己回帰モデルがｋ番目のときの重みである。同じく右辺のｐ（ｙ（ｔ）｜ｘ（ｔ），ｋ）は、入力ベクトルｘ（ｔ）と自己回帰モデルがｋ番と解っているときの入力データｙ（ｔ）の確率密度である。つまり、この式（９）は、Ｋ個の確率モデルの重み付き線形和を示すものでもあるといえる。
【００３６】
パラメータを推定するに際して、式（９）の右辺のｐ（ｙ（ｔ）｜ｘ（ｔ），ｋ）は、次の式（１０）で求まる。すなわち、自己回帰モデルでは、誤差を正規分布と仮定しているため、Ｋ個の自己回帰モデルのうちのｋ番目の自己回帰モデルに属すとしたときの時刻ｔのデータｙ（ｔ）の確率分布（つまり、入力ベクトルｘ（ｔ）と自己回帰モデルが解っているときの入力データｙ（ｔ）の確率密度）は次の式（１０）で求まる。
なお、Ａkはｋ番目の自己回帰モデルのパラメータ、σk2はｋ番目の自己回帰モデルのの分散である。ここで、すべての変数の分散は一定かつ独立という仮定をおいている。
【００３７】
【数１０】
2005258599.tif 000012
【００３８】
（モデルパラメータの推定）
自己回帰混合モデルの未知パラメータの推定は次のように行う。
【００３９】
Ｓｔｅｐ１；  まず、Ｋ*＝１として、次の式（１１）で表される自己回帰モデルの２乗誤差Ｅ1を最小にすることでパラメータＡ1を推定する。
この計算は、最小２乗法で行う。なお、時刻ｔは、初期値がτ＋１であり、終値がＴである。前記したとおり、τは自己回帰モデルの次数である。Ｋ*はループカウンタである。
【００４０】
【数１１】
2005258599.tif 000013
【００４１】
Ｓｔｅｐ２；  １≦ｓ≦Ｋ*を満たすｓを選び、パラメータを次の式（１２）、式（１３）のように設定する。ｓは、パラメータＡを選択するインデックスである。ΔＡは例えば乱数である。ちなみに、最初にＳｔｅｐ２を実行するときは、Ｓｔｅｐ１でＫ*＝１としていることからｓは１になる。２度目以降にＳｔｅｐ２を実行するときは、再びｓが１になることもあり得る。
【００４２】
【数１２】
2005258599.tif 000014
【００４３】
Ｓｔｅｐ３；  ループカウンタのＫ*を１つインクリメントする。
【００４４】
Ｓｔｅｐ４；  データの尤度の対数をとった次の式（１３）で、ＬK*を最大化することによりパラメータ（Ａ、σ、Ｐ）を推定する。なお、右辺のΣのなかの第２項の確率密度は、式（１０）に相当する。
【００４５】
【数１３】
2005258599.tif 000015
【００４６】
Ｓｔｅｐ５；  Ｋ*＜Ｋならば、つまりループカウンタのＫ*が、自己回帰モデルの総数であるＫに満たない場合は、Ｓｔｅｐ２に戻り、Ｓｔｅｐ２以降の処理を行う。
【００４７】
ところで、Ｓｔｅｐ４において、解析的にＬK*を最大にするパラメータを算出することはできない。このため、ＥＭアルゴリズムを反復することによって、つまり次のＥ−Ｓｔｅｐ、Ｍ−Ｓｔｅｐを反復することによって、近似的にＬK*を最大にするパラメータを算出する。
【００４８】
Ｅ−Ｓｔｅｐ；  次の式（１４）を用いて、ｘ（ｔ）とｙ（ｔ）が解っているときのｋ番目の自己回帰モデルの重みＰ（ｋ｜ｘ（ｔ），ｙ（ｔ））を算出する。
【００４９】
【数１４】
2005258599.tif 000016
【００５０】
Ｍ−Ｓｔｅｐ；  次のパラメータ、Ｐ（ｋ）、Ａk、σk2を次の式（１５）〜式（１７）により推定する。なお、Ｐ（ｋ）、Ａk、σk2は、自己回帰モデルがｋ番目のときの、それぞれ、確率密度、パラメータ、分散である。
【００５１】
【数１５】
2005258599.tif 000017
【００５２】
（拡大確率ベクトルの推定）
パラメータの推定が終わると自己回帰混合モデルが構築される。この後は、拡大確率ベクトルのｚ（ｔ）を推定する。ｚ（ｔ）は、次の式（１８）により推定（定義）される。
なお、式（１８）の拡大確率ベクトルｚ（ｔ）は、次の式（１９）の確率ベクトルｚ*（ｔ）に次の式（２０）の自己回帰混合モデルの平均２シグマ値（α〜）を付加したものである。ちなみに、式（１８）は式（５）に対応し、式（１９）は式（４）に対応する。
【００５３】
【数１６】
2005258599.tif 000018
【００５４】
なお、式（１８）と式（１９）において、重みＰについても確率密度ｐについても、ｋは１，２，...，Ｋまでの値をとる。また、確率ベクトルｚ*（ｔ）と拡大確率ベクトルｚ（ｔ）は、各時刻ｔ（ｔ＝１〜Ｔ）ごとにそれぞれ推定される。ちなみに、式（２０）において、分散σk2は式（１７）のものを使用することができる。
ここで、式（１９）は、「入力ベクトルｘと当該入力ベクトルｘに対応する前記データｙとからなるデータＤと各自己回帰モデルとの同時確率をベクトルにしたものを確率ベクトルとして設定」したものである。
【００５５】
（可視化による異常値検出）
与えられた対象データＹからパラメータを推定、自己回帰混合モデルの構築、拡大確率ベクトルの推定が完了すると可視化による異常値検出を行う。本実施形態では、異常値は、各時刻（ｔ＝１〜Ｔ）の拡大確率ベクトルｚ（ｔ）をもとに、ＣｏＰＥ法（Connectivity Preserving Embedding法）によって２又は３次元で可視化するものとする。
なお、ＣｏＰＥ法は、「クロスエントロピ最小化に基づくネットワークデータの埋め込み」、情報処理学会論文誌、44(9):1234-1231,2003に詳しく記載されている。
【００５６】
ＣｏＰＥ法では、まず、時刻ｔがｉのときのデータＤ（ｉ）と時刻ｔがｊのときのデータＤ（ｊ）の類似度ｓi,jを、時刻ｉと時刻ｊの拡大確率ベクトルｚ（ｉ），ｚ（ｊ）のコサイン類似度によって算出する。この際使用する数式は次の式（２１）である。なお、データＤ（ｉ）は、入力ベクトルｘ（ｉ）を入力としてデータｙ（ｉ）を出力とするデータである（ｙとｘとが対になったデータ）。同様に、データＤ（ｊ）は、入力ベクトルｘ（ｊ）を入力としてデータｙ（ｊ）を出力とするデータである。
【００５７】
【数１７】
2005258599.tif 000019
【００５８】
この式（２１）において、分子のｚ（ｉ）´は、時刻ｉにおける式（１８）の拡大確率ベクトルｚ（ｉ）の転置を示す。また、分母は拡大確率ベクトルｚ（ｉ）の絶対値と同ｚ（ｊ）の絶対値の積を示す。
【００５９】
ここで、ｒiをデータＤ（ｉ）の座標とする。座標（ｉがτ＋１〜Ｔまでの全座標）を、ｒiとｒjの確率密度ｐ（ｒi，ｒj）と類似度ｓi,jのクロスエントロピの和に正則化項を加えた次の式（２２）を最小化することにより算出する。
【００６０】
【数１８】
2005258599.tif 000020
【００６１】
なお、時刻ｉは、初期値がτ＋１であり、終値がＴ−１である。また、時刻ｊは、初期値がτ＋２であり、終値がＴである。μは正則化項の重みで、例えば可視化結果に基づいて適宜設定される。‖ｒi‖2は、座標ｒiの原点からの距離である。
【００６２】
ここで、ｐ（ｒi，ｒj）のクロスエントロピは次の式（２３）により算出され、この式（２３）で必要とするｐ（ｒi，ｒj）は次の式（２４）により算出される。
【００６３】
【数１９】
2005258599.tif 000021
【００６４】
この座標ｒiを２次元又は３次元でグラフ化して、可視化する。これにより、異常値が視覚的に判別可能になる。
【００６５】
≪異常値可視化装置の構成≫
次に、前記した原理を踏まえて、与えられたデータ（対象データＹ）に含まれる異常値を、コンピュータにより可視化する異常値可視化装置を具体的に説明する。
図１は、本発明の一実施形態の異常値可視化装置の構成図である。
【００６６】
図１に示すように、異常値可視化装置Ａは、演算装置１、入力装置２、記憶装置３、表示装置４がバスに接続された構成をしている。演算装置１は、ハードウェア的には、演算手段としてのＣＰＵ（Central Processing Unit）や情報を記憶するメモリとしてのＲＡＭ（Random Access Memory）などから構成されている。この演算装置１は、ソフトウェア構成として、モデルパラメータ推定部１１、拡大確率ベクトル推定部１２、異常値可視化部１３を備えている。各部１１，１２，１３の機能は、後で詳しく説明する。なお、モデルパラメータ推定部１１、拡大確率ベクトル推定部１２、異常値可視化部１３は、その基になるプログラムが記憶装置３に記憶されており、異常値可視化装置Ａを起動した際、ＣＰＵにより記憶装置３からＲＡＭ上に読み出されて各部１１，１２，１３として機能するものとする。ＣＰＵは、各種情報を記憶するメモリであるＲＡＭを作業領域として、演算結果や演算経過を随時ＲＡＭに記憶しつつ各種演算を行う。もちろん、ＣＰＵは、ＣＰＵが内蔵するキャッシュメモリも作業領域としつつ各種演算を行う。
【００６７】
入力装置２は、キーボード、ネットワークインタフェイスカード、ＦＤＤのようなディスクドライブ装置などから構成される入力手段であり、前記した対象データＹ（経済時系列のデータ）や自己回帰モデルの次数τ、自己回帰モデルの総数Ｋなどの各種データを異常値可視化装置Ａに入力するために用いられる。記憶装置３は、ハードディスク装置などから構成される記憶手段である。本実施形態では、入力装置２から入力された各種データは、記憶装置３に記憶されるものとする。また、記憶装置３は、前記したように、各部１１，１２，１３の基になるプログラムを記憶しているものとする。表示装置４は、例えばグラフィックボード及びそれに接続された液晶モニタであり、異常値可視化部１３により可視化されたデータを視覚的に表示する。
【００６８】
次に、演算装置１にソフトウェア的に（プログラムモジュールとして）構成されるモデルパラメータ推定部１１、拡大確率ベクトル推定部１２、異常値可視化部１３の機能を説明する。
【００６９】
（モデルパラメータ推定部）
演算装置１のモデルパラメータ推定部１１は、式（９）及び式（１０）で示す自己回帰混合モデルのパラメータを推定して、自己回帰混合モデルを構築する機能を有する。このため、モデルパラメータ推定部１１は、（ａ）記憶装置３に記憶された式（６）で示す対象データＹと次数τ、その他を演算装置１に読み込む処理、（ｂ）演算装置１に読み込んだ対象データＹと次数τとから、式（７）の自己回帰モデルへの入力ベクトルを算出する処理、（ｃ）対象データＹと入力ベクトルｘ（ｔ）から、式（１１）〜式（１７）を用い、自己回帰混合モデルのパラメータなど（Ａ、σ、Ｐ（ｋ））を推定する処理を行う機能を有する。なお、読み込んだ対象データＹなどや演算の結果・途中経過は、適宜図示しないＲＡＭに記憶して、後段の拡大確率ベクトル推定部１２での処理に活用するものとする。本実施形態では、重みのＰ（ｋ）もパラメータの１つとして説明する。
【００７０】
（拡大確率ベクトル推定部）
拡大確率ベクトル推定部１２は、前記した式（１８）〜式（２０）を用いて、また、適宜ＲＡＭに記憶されているモデルパラメータ推定部１１の演算の結果や途中経過を用いて、拡大確率ベクトルｚ（ｔ）を各時刻ｔ（ｔ＝１〜Ｔ）ごとにそれぞれ推定する処理を行う機能を有する。
【００７１】
（異常値可視化部）
異常値可視化部１３は、前記した式（２１）〜式（２４）を用いて対象データＹに含まれる異常値をＣｏＰＥ法により可視化する処理を行う機能を有する。
【００７２】
≪異常値可視化装置の動作≫
次に、フローチャートを用いて、前記説明した構成を有する異常値可視化装置Ａの動作を説明する。
図２は、読み込んだデータから異常値を可視化するまでの処理の概要を示したフローチャートである。この図２に示すように、異常値可視化装置Ａは、ステップＳ１００で、記憶装置３から対象データＹを読み込み、「自己回帰混合モデルの構築処理」を行う。なお、この処理を行うのは、図１に示すモデルパラメータ推定部１１である。次に、ステップＳ２００で、「拡大確率ベクトルの推定処理」を行う。なお、この処理を行うのは、図１に示す拡大確率ベクトル推定部１２である。拡大確率ベクトルの推定処理後は、Ｓ３００で、「異常値の可視化処理」を行う。なお、この処理を行うのは、図１に示す異常値可視化部１３である。
【００７３】
（自己回帰混合モデルの構築）
次に、図３のフローチャートを参照して、図２のステップＳ１００の「自己回帰混合モデルの構築処理」を詳細に説明する（適宜図１など参照）。ここで行うのは、式（９）及び式（１０）で示す自己回帰混合モデルのパラメータを推定して、自己回帰混合モデルを構築する処理である。以下の説明における動作主体は、モデルパラメータ推定部１１である。なお、モデルパラメータ推定部１１は、ＣＰＵやＲＡＭなどから構成される演算装置１で機能するプログラムモジュールである。
ちなみに、この図３のフローチャートは、コンピュータでの演算処理の順序を考慮して表現しているので、「可視化方法の原理・数式の説明」のところで説明した数式の登場順序とは必ずしも一致しない。この点は、以降の各フローチャート（図４、図５）においても同じである。
【００７４】
まず、図３に示すように、モデルパラメータ推定部１１は、記憶装置３から演算装置１に、式（６）で示す構造をした対象データＹを読み込む（Ｓ１０１）。対象データＹは、前記したとおり６次元の経済時系列データである。読み込んだ対象データＹについて、式（７）で示す構造をした自己回帰モデルへの入力データｘ（ｔ）の作成を行う（Ｓ１０２）。なお、式（７）で使われるτは、前記したとおり自己回帰モデルの次数であり、このτは、別途演算装置１に読み込まれているものとする。
【００７５】
次に、ループカウンタＫ*を初期化して１にする（Ｓ１０３）。なお、ループカウンタＫ*は、１〜Ｋまでの整数値をとる。このステップＳ１０３は、すべての自己回帰モデルのパラメータを推定するための準備ステップである。
【００７６】
続いて、前記ステップＳ１０１で読み込んだ対象データＹに含まれるデータｙ（ｔ）と前記ステップＳ１０２で作成した入力ベクトルｘ（ｔ）を用いて、式（１１）の２乗最小誤差Ｅ1を最小にすることで、パラメータＡ1を推定する（Ｓ１０４）。ここまでは、自己回帰モデルのパラメータの推定で、次のステップＳ１０５以降は、自己回帰混合モデルのパラメータの推定である。
【００７７】
パラメータＡ1を推定した後は、１≦ｓ≦Ｋ*を満たす任意のインデックスｓを選び、パラメータを式（１２）のように設定する（Ｓ１０５）。なお、ループカウンタのＫ*が１のときは、必ずｓは１になる。式（１２）のＡK*+1＝Ａｓ−ΔＡは、添字の部分でＫ*＋１としていることから、次のＡK*の値を設定していることになる。
【００７８】
パラメータのＡを式（１２）のように設定すると、ループカウンタＫ*を１つインクリメントする（Ｓ１０６）。そして、前記したＥＭアルゴリズムのＥ−ｓｔｅｐを実行する（Ｓ１０７）。つまり、式（１４）により、データｙ（ｔ）と入力ベクトルｙ（ｔ）が解っているときのｋ番目の自己回帰モデルの重みを算出する（Ｓ１０７）。
ちなみに、式（１１）と式（１３）では、各時刻ｔ＝τ＋１〜Ｔまでの値を使って、１つのＥ1、１つのＬK*が算出される。一方、式（１４）については、複数個のＰ（ｋ｜ｘ（ｔ）、ｙ（ｔ））が算出される。
【００７９】
次に、前記したＥＭアルゴリズムのＭ−ｓｔｅｐを実行する（Ｓ１０８）。つまり、式（１５）によりＰ（ｋ）を、式（１６）によりσk2を、式（１７）によりＡkをそれぞれ推定する。ちなみに、ステップＳ１０７で算出した重みＰ（ｋ｜ｘ（ｔ），ｙ（ｔ））を用いて式（１５）〜式（１７）を実行することから、実行するごとに異なるＰ（ｋ）、σk2、Ａkが推定される。
【００８０】
ステップＳ１０９とステップＳ１１０では、式（１３）を用い、推定したパラメータを検証する。パラメータの検証は、式（１３）によりデータの尤度の対数をとったＬK*を算出し（Ｓ１０９）、これが最大化したときのパラメータを検証済みのパラメータとして採用するものである。なお、本実施形態では、ＬK*の最大化をＬK*が収束したか否かにより判断する（Ｓ１１０）。従って、ＬK*が収束したならば（Ｓ１１０においてＹｅｓ）、インクリメントされたＫ*におけるパラメータの推定が完了する。一方、ＬK*が収束しないならば（Ｓ１１０においてＮｏ）、ステップＳ１０７に戻ってステップＳ１０７以降を再度実行し、新たなパラメータを推定して検証を行う。
【００８１】
ちなみに、本実施形態では、ＬK*が収束したか否かの判断は、ステップＳ１０９において算出したＬK*の前回値と今回値を比較し、その差の絶対値が所定値未満（以下）ならば収束したと判断するものである。もちろん、ＬK*が収束したか否かを判断するのはＬK*が最大化したか否かを判断するための一例であり、ＬK*が最大化したか否かの判断の仕方については、本実施形態以外の判断の仕方でもよい。
【００８２】
ステップＳ１１１では、ループカウンタのＫ*が自己回帰モデルの総数であるＫよりも小さいか（Ｋ*＜Ｋ）否かを判断する。Ｋ*＜Ｋならば（ステップＳ１１１においてＹｅｓ）、ステップＳ１０５に戻り、ステップＳ１０５以降を実行する。つまり、次のＫ*（次の自己回帰モデル）についてのパラメータの推定を行う。一方、Ｋ*＜Ｋでないならば、つまりＫ*＝Ｋならば（ステップＳ１１１においてＮｏ）、自己回帰混合モデルのすべてのパラメータ（Ｐ（ｋ），σk2、Ａk）が推定されたので、この図３のフローチャートの処理を終了する。これにより、自己回帰混合モデルが構築される。
【００８３】
（拡大確率ベクトルの推定処理）
次に、図４のフローチャートを参照して、図２のステップＳ２００の「拡大確率ベクトルの推定処理」を詳細に説明する（適宜図など１参照）。以下の説明における動作主体は、拡大確率ベクトル推定部１２である。なお、拡大確率ベクトル推定部１２は、演算装置１で機能するプログラムモジュールである。
【００８４】
まず、図４に示すように、拡大確率ベクトル推定部１２は、ステップＳ２０１でｔ＝０として、時刻ｔを初期化する。これは、ｔ＝１〜Ｔまでのすべての時刻ｔについて、拡大確率ベクトルを推定するための準備である。ステップＳ２０２で、時刻ｔをインクリメントする（データｙ（０）が存在しないため）。ステップＳ２０３で、「データｙ（ｔ）と入力ベクトルｘ（ｔ）」と「各自己回帰モデル（ｋ＝１，... ，Ｋ）」との同時確率（つまりデータＤ（ｔ）と各確率モデルとの同時確率）をベクトルにしたものを、式（１８）の確率ベクトルｚ*（ｔ）として設定する。なお、式（１８）は、前記したステップＳ１００での演算結果に基づいて設定される。
【００８５】
次に、自己回帰混合モデルの平均２シグマ値を、ステップＳ１００での演算結果を利用して式（２０）により算出する（Ｓ２０４）。そして、この平均２シグマ値を確率ベクトルｚ*（ｔ）に付加して、式（１９）の拡大確率ベクトルｚ（ｔ）を推定する（Ｓ２０５）。次に、時刻ｔが時系列の総数であるＴより小さいか（ｔ＜Ｔか）を判断する（Ｓ２０６）。ｔがＴに満たない場合（Ｓ２０６でＹｅｓ）は、ステップＳ２０２に戻ってステップＳ２０２以降を実行する。つまり、次の時刻ｔの拡大確率ベクトルｚ（ｔ）を推定する。一方、ステップＳ２０６で、ｔ＝Ｔとなった場合は（Ｙｅｓ）、すべての時刻ｔについての拡大確率ベクトルが推定されたので、この図４のフローチャートの処理を終了する。
【００８６】
（異常値の可視化処理）
次に、図５のフローチャートを参照して、図２のステップ３００の「異常値の可視化処理」を詳細に説明する（適宜図１など参照）。以下の説明における動作主体は、異常値可視化部１３である。なお、異常可視化部１３は、ＣＰＵやＲＡＭなどから構成される演算装置１で機能するプログラムモジュールである。
【００８７】
まず、図５に示すように、異常値可視化部１３は、ステップＳ３０１により、時刻ｉのデータＤ（ｉ）と時刻ｉとは異なる時刻ｊのデータＤ（ｊ）の類似度ｓi,jを、前記したステップＳ２００で算出した時刻ｉと時刻ｊの拡大確率ベクトルｚ（ｉ），ｚ（ｊ）を用いて算出する（式（２１）参照）。
【００８８】
次に、ステップＳ３０２により、式（２２）〜式（２４）を用いてｒiをデータＤ（ｉ）の座標とする。座標ｒｉを、式（２３）・式（２４）を用いつつ、確率密度ｐ（ｒｉ，ｒｊ）とのクロスエントロピの和に正則化項を加えた式（２２）を最小化することにより求める。
【００８９】
以上の結果、他のデータから離れたところに位置するものが異常値であり、異常値が可視化される。
【００９０】
≪可視化結果≫
本発明の可視化方法などの有効性を、図６を参照しつつ、経済時系列データを用いた実施例で示す（適宜図１など参照）。なお、図６の（ａ）から（ｄ）は、本実施形態の異常値可視化装置により異常値を３次元で可視化した結果を示す図である（図面では平面的に表示されている）。
【００９１】
用いるデータ（対象データＹ）は、前記したとおり、１９８３年１月から２００２年１２月までの「日本のマネタリーベース」、「国際金利」、「卸売物価指数」、「機械受注」、「鉱工業生産指数」、「円ドル為替レート」の月ごとの時系列データで、変数の数ｄ＝６、データ数Ｔ＝２４０である。この対象データＹに対して、前記した自己回帰混合モデルを用いた時系列データの異常値の可視化方法を採用し、自己回帰モデルの次数をτ＝１、自己回帰モデルの総数Ｋ＝１〜４として可視化を行った。ちなみに、Ｋが１の場合が図６（ａ）であり、Ｋが２の場合が図６（ｂ）であり、Ｋが３の場合が図６（ｃ）であり、Ｋが４の場合が図６（ｄ）である。各図中の１点が１つの月である。なお、異常値可視化装置Ａは、自己回帰モデルごとに点の種類が異なるように表示する機能を有しており、自己回帰モデルが複数ある場合でも、表示されている点が何番目の自己回帰モデルにより生成されたものであるのかが、視覚的に解るようになっている。例えば、図６（ｂ）は、１番目の自己回帰モデルが○印の点で示してあり、２番目の自己回帰モデルが×印の点で示してある。図６（ｃ）は、３番目の自己回帰モデルが△印の点で示してある。また、図６（ｄ）は、４番目の自己回帰モデルが□の点で示してある。
【００９２】
図６に示すように、この期間に起こった大きな経済事件として、１９９８年８月のロシア危機がある。異常値の可視化結果を見ると、自己回帰モデルの総数Ｋ＝１〜４のすべての場合において、ロシア危機の翌月である１９９８年９月が他の多くの月から大きく離れており、異常値であることが解る。しかしながら、自己回帰モデルの総数Ｋが１の場合は、他の月も１９９８年９月近辺に位置されており、異常値が複数存在する。自己回帰モデルの総数Ｋを２、３としていっても異常値とみなされる月は複数あるが、自己回帰モデルの総数Ｋが４の場合は、１９９８年９月のみ異常値となっている。
【００９３】
これは、自己回帰モデルの総数Ｋが１の場合、自己回帰混合モデルは線形モデルであるため、データが非線形性を有していた場合、適切にモデル化することができない。そのため、本来は異常値ではないにもかかわらず、モデルの表現能力が低いために異常値とされることがある。しかし、自己回帰モデルの総数Ｋを増やしていくことによって、モデルの表現能力が高まり、本来の異常値を検出することができる。
このように、本発明によれば、表現力のあるモデルで可視化するので、従来手法よりも良好な異常値の検出が実現できる。
【００９４】
また、従来の異常値検出手法は、データ空間のなかでの平均からのデータの乖離をもとにしている。しかし、本実施形態では、各データをモデルのパラメータ空間に写像することにより、換言すると、式（２）や式（９）のように、データをある確率モデルから生成されるようにすることにより、対象データＹといった生データが有する特徴をより明確にすることができる。
【００９５】
なお、異常値可視化装置は、プログラムにより可視化方法を実行するものであり、この可視化方法をコンピュータに実行させるプログラムは、ＣＤＲ−ＲＯＭやデジタル多目的ディスクなどの記憶媒体に記憶されて流通され、また、通信回線を介して流通され、コンピュータにインストールされて機能する。この可視化方法の実施形態としては、可視化方法を実行するプログラムがインストールされたコンピュータ（つまり異常値可視化装置）が通信回線を介して対象データを取得し、これを可視化し、その結果を、通信回線を介して回答するようにしてもよい。また、拡大確率ベクトルを用いての異常値の可視化処理（図２のステップＳ３００）は、前記ＣｏＰＥ法以外に種々の手法を適用することができる。
また、確率モデルについて、前記実施形態では、自己回帰モデルを用いた例を説明したが、この自己回帰モデルに限定されるものではなく、正規分布モデルや多項分布モデルといった、他の確率モデルを用いてもよい。
また、前記した実施形態では、所定のデータの一例として異常値を可視化したが、所定のデータとして、例えば好ましい値を可視化するようにしてもよい。
つまり、本発明は、前記した実施形態に限定されることなく、その技術思想の及ぶ範囲で幅広く変形実施することができる。
【図面の簡単な説明】
【００９６】
【図１】本発明の一実施形態の異常値可視化装置の構成図である。
【図２】読み込んだデータから異常値を可視化するまでの処理の概要を示したフローチャートである。
【図３】図２のステップＳ１００の「自己回帰混合モデルの構築処理」を示すフローチャートである。
【図４】図２のステップＳ２００の「拡大確率ベクトルの推定処理」を示すフローチャートである。
【図５】図２のステップＳ３００の「異常値の可視化処理」を示すフローチャートである。
【図６】図１の異常値可視化装置により異常値を可視化した結果を示す図であり、（ａ）はＫが１の場合を示し、（ｂ）はＫが２の場合を示し、（ｃ）はＫが３の場合を示し、（ｄ）はＫが４の場合を示す。
【符号の説明】
【００９７】
Ａ      異常値可視化装置
１      演算装置（演算手段）
１１    モデルパラメータ推定部
１２    拡大確率ベクトル推定部
１３    異常値可視化部
２      入力装置
３      記憶装置
４      表示装置
【図１】
2005258599.tif 000022
【図２】
2005258599.tif 000023
【図３】
2005258599.tif 000024
【図４】
2005258599.tif 000025
【図５】
2005258599.tif 000026
【図６】
2005258599.tif 000027