(57)【要約】【課題】複数のオブジェクトの間の関係の有無を示す関係データの時系列から、オブジェクトのクラスタの存在を推定し、クラスタの時間変化を推定する。【解決手段】公知のＩＲＭ（Infinite Relational Model）の数学的モデルが所属クラスタｚ、クラスタ間関連度η、混合比β、および関連データｘによって表される。そのＩＲＭを拡張して、クラスタの時間変化を表すパラメータπを導入し、所属クラスタｚの推定方法も時間に従って考慮可能なように改造した。そして、図１に示す所属クラスタｚ推定演算（ステップＳ１０４）、時間変化情報π推定演算（ステップＳ１０５）、クラスタ間関連度η推定演算（ステップＳ１０６）、および混合比β推定演算（ステップＳ１０７）を、終了条件を満足するまで、繰り返し実行することによって、各パラメータの推定値が算出される。【選択図】図１
2011048583.tif 000002
【特許請求の範囲】
【請求項１】
複数のオブジェクトの間の関係の有無を表す関係データを用いて、オブジェクトのクラスタリングを実行するクラスタリング装置であって、
前記クラスタリング装置は、
前記関係データをクラスタリングする関数の一つである無限関係モデル（ＩＲＭ、Infinite Relational Model）において算出される混合比βおよびクラスタｋとクラスタｌとの間の関係の強さを示すクラスタ間関連度ηｋ，ｌと、クラスタｋとクラスタｌとの間の関連の有無を示す関係データを所定の時間間隔で観測した関係データｘｔ，ｋ，ｌ、時刻ｔ−１においてクラスタｋに所属していたオブジェクトが次の時刻ｔにどのクラスタに所属しやすいかを示す時間変化情報πｔ，ｋ、および時刻ｔにおいてオブジェクトｉが所属するクラスタを示す所属クラスタｚｔ，ｉと、ハイパーパラメータα０，κと、クラスタ数Ｋとを記憶する記憶部と、
前記記憶部から前記所属クラスタｚｔ，ｉ、前記時間変化情報πｔ，ｋ、前記混合比β、前記ハイパーパラメータα０，κ、および前記クラスタ数Ｋとを取得し、所定期間ｔ＝１〜Ｔにおいて、ｚｔ−１，ｉ＝ｋかつｚｔ，ｉ＝ｌとなるオブジェクトの数をｍｔ，ｋ，ｌとしたとき、ディリクレ分布Dirichlet（α０β１＋ｍｔ，ｋ，ｌ，・，α０βｋ＋ｍｔ，ｋ，ｌ＋κ，・，α０βＫ＋ｍｔ，ｋ，Ｋ，α０（１−ΣＫｋ＝１βｋ））からサンプリングして前記時間変化情報πｔ，ｋを算出し、当該算出した時間変化情報πｔ，ｋによって前記記憶部に記憶してある時間変化情報πｔ，ｋを更新し記憶する時間変化情報推定部と、
前記記憶部から前記所属クラスタｚｔ，ｉ、前記時間変化情報πｔ，ｋ、前記クラスタ間関連度ηｋ，ｌ、前記混合比β、所定期間ｔ＝１〜Ｔの前記関係データｘｔ，ｋ，ｌ、および前記クラスタ数Ｋを取得し、取得した前記混合比β、前記時間変化情報πｔ，ｋ、前記クラスタ間関連度ηｋ，ｌ、前記関係データｘｔ，ｋ，ｌ、およびクラスタ数Ｋを、式（１)、式（２）、式（３）、式（４）に適用して、式（１）からｕｔ，ｊをサンプリングし、メッセージ変数ｐｔ，ｉ，ｋを式（３）で定義したとき、ｔ＝１からｔ＝Ｔまで順番に式（２）を用いてメッセージ変数ｐｔ，ｉ，ｋを算出し、当該算出したメッセージ変数ｐｔ，ｉ，ｋに対して、ｔ＝Ｔからｔ＝１まで順番に式（４）を用いてメッセージ変数ｐｔ，ｉ，ｋを算出し、ｐ（ｚｔ，ｉ＝ｋ｜ｚｔ−１，ｉ＝ｌ）が０とならない場合の所属クラスタｚｔ，ｉを算出し、当該算出した所属クラスタｚｔ，ｉによって前記記憶部に記憶してある所属クラスタｚｔ，ｉを更新し記憶する所属クラスタ推定部と、
前記クラスタ間関連度ηｋ，ｌを前記無限関係モデルによって算出し、前記記憶部に記憶するクラスタ間関連度推定部と、
前記混合比βを前記無限関係モデルによって算出し、前記記憶部に記憶する混合比推定部と、
前記混合比推定部、前記クラスタ間関連度推定部、前記時間変化情報推定部、および前記所属クラスタ推定部における演算を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す終了判定部と、
を備える
ことを特徴とするクラスタリング装置。
【数１】
2011048583.tif 000020
【請求項２】
前記クラスタリング装置は、
さらに、ハイパーパラメータγ，ξ，Ψを記憶する前記記憶部を備え、
前記クラスタ間関連度推定部は、前記記憶部から前記所属クラスタｚｔ，ｉ、前記クラスタ間関連度ηｋ，ｌ、前記ハイパーパラメータξ，Ψ、前記クラスタ数Ｋ、および所定期間ｔ＝１〜Ｔの前記関係データｘｔ，ｋ，ｌを取得して、ｚｔ，ｉ＝ｋかつｚｔ，ｊ＝ｌとなる（ｔ，ｉ，ｊ）の組の数をＮｋ，ｌとし、当該Ｎｋ，ｌの中の前記関係データｘｔ，ｋ，ｌが関係有りを示す数をｎｋ，ｌとしたとき、ベータ分布Beta（ξ＋ηｋ，ｌ，Ψ＋Ｎｋ，ｌ−ｎｋ，ｌ）からサンプリングして前記クラスタ間関連度ηｋ，ｌを算出し、当該算出したクラスタ間関連度ηｋ，ｌによって前記記憶部に記憶してあるクラスタ間関連度ηｋ，ｌを更新し記憶し、
前記混合比推定部は、前記記憶部から前記所属クラスタｚｔ，ｉ、前記混合比β、および前記ハイパーパラメータγ，α０，κを取得し、取得した前記所属クラスタｚｔ，ｉ、前記混合比β、および前記ハイパーパラメータα０，κを、式（５）および式（６）に適用して、それぞれ補助変数Ｒｔ，ｋ，ｌおよび補助変数Ｏｔ，ｋを算出し、前記算出した補助変数Ｒｔ，ｋ，ｌおよび補助変数Ｏｔ，ｋを式（７）に適用して補助変数＾Ｒｔ，ｋ，ｌを算出し、ディリクレ分布Dirichlet（Σｔ，ｋ＾Ｒｔ，ｋ，１，Σｔ，ｋ＾Ｒｔ，ｋ，２，・，Σｔ，ｋ＾Ｒｔ，ｋ，γ）からサンプリングして前記混合比βを算出し、当該算出した混合比βによって前記記憶部に記憶してある混合比βを更新し記憶する
ことを特徴とする請求項１に記載のクラスタリング装置。
【数２】
2011048583.tif 000021
【請求項３】
前記クラスタリング装置は、演算結果を表示する表示装置と接続され、
前記終了判定部は、前記終了条件を満足した場合、前記混合比β、前記クラスタ間関連度ηｋ，ｌ、前記時間変化情報πｔ，ｋ、および前記所属クラスタｚｔ，ｉのいずれか一つまたはいずれかの組み合わせを前記表示装置に出力する
ことを特徴とする請求項１または請求項２に記載のクラスタリング装置。
【請求項４】
複数のオブジェクトの間の関係の有無を表す関係データを用いて、オブジェクトのクラスタリングを実行するクラスタリング装置において用いられるクラスタリング方法であって、
前記クラスタリング装置は、
前記関係データをクラスタリングする関数の一つである無限関係モデル（ＩＲＭ、Infinite Relational Model）において算出される混合比βおよびクラスタｋとクラスタｌとの間の関係の強さを示すクラスタ間関連度ηｋ，ｌと、クラスタｋとクラスタｌとの間の関連の有無を示す関係データを所定の時間間隔で観測した関係データｘｔ，ｋ，ｌ、時刻ｔ−１においてクラスタｋに所属していたオブジェクトが次の時刻ｔにどのクラスタに所属しやすいかを示す時間変化情報πｔ，ｋ、および時刻ｔにおいてオブジェクトｉが所属するクラスタを示す所属クラスタｚｔ，ｉと、ハイパーパラメータα０，κと、クラスタ数Ｋとを記憶する記憶部と処理部とを備え、
前記処理部は、
前記記憶部から前記所属クラスタｚｔ，ｉ、前記時間変化情報πｔ，ｋ、前記混合比β、前記ハイパーパラメータα０，κ、および前記クラスタ数Ｋとを取得し、所定期間ｔ＝１〜Ｔにおいて、ｚｔ−１，ｉ＝ｋかつｚｔ，ｉ＝ｌとなるオブジェクトの数をｍｔ，ｋ，ｌとしたとき、ディリクレ分布Dirichlet（α０β１＋ｍｔ，ｋ，ｌ，・，α０βｋ＋ｍｔ，ｋ，ｌ＋κ，・，α０βＫ＋ｍｔ，ｋ，Ｋ，α０（１−ΣＫｋ＝１βｋ））からサンプリングして前記時間変化情報πｔ，ｋを算出し、当該算出した時間変化情報πｔ，ｋによって前記記憶部に記憶してある時間変化情報πｔ，ｋを更新し記憶する時間変化情報推定ステップと、
前記記憶部から前記所属クラスタｚｔ，ｉ、前記時間変化情報πｔ，ｋ、前記クラスタ間関連度ηｋ，ｌ、前記混合比β、所定期間ｔ＝１〜Ｔの前記関係データｘｔ，ｋ，ｌ、および前記クラスタ数Ｋを取得し、取得した前記混合比β、前記時間変化情報πｔ，ｋ、前記クラスタ間関連度ηｋ，ｌ、前記関係データｘｔ，ｋ，ｌ、およびクラスタ数Ｋを、式（１)、式（２）、式（３）、式（４）に適用して、式（１）からｕｔ，ｊをサンプリングし、メッセージ変数ｐｔ，ｉ，ｋを式（３）で定義したとき、ｔ＝１からｔ＝Ｔまで順番に式（２）を用いてメッセージ変数ｐｔ，ｉ，ｋを算出し、当該算出したメッセージ変数ｐｔ，ｉ，ｋに対して、ｔ＝Ｔからｔ＝１まで順番に式（４）を用いてメッセージ変数ｐｔ，ｉ，ｋを算出し、ｐ（ｚｔ，ｉ＝ｋ｜ｚｔ−１，ｉ＝ｌ）が０とならない場合の所属クラスタｚｔ，ｉを算出し、当該算出した所属クラスタｚｔ，ｉによって前記記憶部に記憶してある所属クラスタｚｔ，ｉを更新し記憶する所属クラスタ推定ステップと、
前記クラスタ間関連度ηｋ，ｌを前記無限関係モデルによって算出し、前記記憶部に記憶するクラスタ間関連度推定ステップと、
前記混合比βを前記無限関係モデルによって算出し、前記記憶部に記憶する混合比推定ステップと、
前記混合比推定ステップ、前記クラスタ間関連度推定ステップ、前記時間変化情報推定ステップ、および前記所属クラスタ推定ステップにおける演算を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返し演算させる終了判定ステップと、
を実行することを特徴とするクラスタリング方法。
【数１】
2011048583.tif 000022
【請求項５】
前記クラスタリング装置は、
さらに、ハイパーパラメータγ，ξ，Ψを記憶する前記記憶部を備え、
前記処理部は、
前記クラスタ間関連度推定ステップにおいて、前記記憶部から前記所属クラスタｚｔ，ｉ、前記クラスタ間関連度ηｋ，ｌ、前記ハイパーパラメータξ，Ψ、前記クラスタ数Ｋ、および所定期間ｔ＝１〜Ｔの前記関係データｘｔ，ｋ，ｌを取得して、ｚｔ，ｉ＝ｋかつｚｔ，ｊ＝ｌとなる（ｔ，ｉ，ｊ）の組の数をＮｋ，ｌとし、当該Ｎｋ，ｌの中の前記関係データｘｔ，ｋ，ｌが関係有りを示す数をｎｋ，ｌとしたとき、ベータ分布Beta（ξ＋ηｋ，ｌ，Ψ＋Ｎｋ，ｌ−ｎｋ，ｌ）からサンプリングして前記クラスタ間関連度ηｋ，ｌを算出し、当該算出したクラスタ間関連度ηｋ，ｌによって前記記憶部に記憶してあるクラスタ間関連度ηｋ，ｌを更新し記憶し、
前記混合比推定ステップにおいて、前記記憶部から前記所属クラスタｚｔ，ｉ、前記混合比β、および前記ハイパーパラメータγ，α０，κを取得し、取得した前記所属クラスタｚｔ，ｉ、前記混合比β、および前記ハイパーパラメータα０，κを、式（５）および式（６）に適用して、それぞれ補助変数Ｒｔ，ｋ，ｌおよび補助変数Ｏｔ，ｋを算出し、前記算出した補助変数Ｒｔ，ｋ，ｌおよび補助変数Ｏｔ，ｋを式（７）に適用して補助変数＾Ｒｔ，ｋ，ｌを算出し、ディリクレ分布Dirichlet（Σｔ，ｋ＾Ｒｔ，ｋ，１，Σｔ，ｋ＾Ｒｔ，ｋ，２，・，Σｔ，ｋ＾Ｒｔ，ｋ，γ）からサンプリングして前記混合比βを算出し、当該算出した混合比βによって前記記憶部に記憶してある混合比βを更新し記憶する
ことを特徴とする請求項４に記載のクラスタリング方法。
【数２】
2011048583.tif 000023
【請求項６】
前記クラスタリング装置は、演算結果を表示する表示装置と接続され、
前記処理部は、
前記終了条件を満足した場合、前記混合比β、前記クラスタ間関連度ηｋ，ｌ、前記時間変化情報πｔ，ｋ、および前記所属クラスタｚｔ，ｉのいずれか一つまたはいずれかの組み合わせを前記表示装置に出力する
ことを特徴とする請求項４または請求項５に記載のクラスタリング方法。
【請求項７】
請求項４ないし請求項６のいずれか一項に記載のクラスタリング方法を、コンピュータとしてのクラスタリング装置に実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、複数のオブジェクトの間の関係の有無を示す関係データから、オブジェクトのクラスタの存在を推定し、クラスタの生成、隆盛、消滅等のクラスタの時間変化を推定する技術に関する。
【背景技術】
【０００２】
複数のオブジェクトの間の関係（例えば、人と人との関係やインターネットにおけるリンク関係等）の有無を示す関係データを用いて、オブジェクトのクラスタ（グループ）間の関係を表現し、そのクラスタ間の関係を最も良く表すオブジェクトのクラスタを求める技術として、非特許文献１に記載のStochastic Block Model（ＳＢＭ）や、そのＳＢＭを拡張した非特許文献２に記載のInfinite Relational Model（ＩＲＭ）とが知られている。ＳＢＭがクラスタ数を予め設定してからクラスタリングを実行するのに対して、ＩＲＭは、クラスタ数を事前に決めることなく、関係データに基づいて最適なクラスタ数を自動的に推定することができる。
【０００３】
ＩＲＭは、オブジェクトｉとオブジェクトｊとの関係データＸ＝｛ｘｉ，ｊ｝に基づいて、各オブジェクトｉ，ｊ間の関係を複数のクラスタに分類する。このクラスタの分類では、ノンパラメトリックベイズモデルの一つである非特許文献３に記載のDirichlet Process Mixture（ＤＰＭ）を適用することによって、クラスタ数の推定が可能になっている。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】K.Nowicki and T.A.B.Snijders,“Estimation and Prediction for Stochastic Blockstructures”, Journal of the American Statistical Association, Vol.96, No.455, p.1077-1087, 2001
【非特許文献２】C.Kemp, J.B.Tenenbaum, T.L.Griffiths, T.Yamada and N.Ueda, “Learning Systems Of Concepts With An Infinite Relational Model”, Proceedings of the 21st National Conference on Artificial Intelligence, 2006
【非特許文献３】T.S.Ferguson,“A Bayesian Analysis of Some Nonparametric Problems”, The Annals of Statistics, Vol.1, No.2, p.353-355, 1973.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、ＳＢＭやＩＲＭでは、それらのモデルの構成上、クラスタの時間変化を扱うことができない。すなわち、ＳＢＭやＩＲＭは、ある一時刻のスナップショットの関係データだけを用いたり、または一定時間の平均的な関係を示すデータを用いたりすることしかできない。そのため、人と人との関係やインターネットにおけるリンク関係等のようにオブジェクト間の関係が時間的に変化するデータに対して、時間方向の情報を解析することができなかった。
【０００６】
そこで、本発明の課題は、複数のオブジェクトの間の関係の有無を示す関係データの時系列から、オブジェクトのクラスタの存在を推定し、クラスタの時間変化を推定する技術を提供することを目的とする。
【課題を解決するための手段】
【０００７】
本発明は、複数のオブジェクトの間の関係の有無を表す関係データを用いて、オブジェクトのクラスタリングを実行するクラスタリング装置であって、前記クラスタリング装置が、前記関係データをクラスタリングする関数の一つである無限関係モデル（ＩＲＭ、Infinite Relational Model）において算出される混合比βおよびクラスタｋとクラスタｌとの間の関係の強さを示すクラスタ間関連度ηｋ，ｌと、クラスタｋとクラスタｌとの間の関連の有無を示す関係データを所定の時間間隔で観測した関係データｘｔ，ｋ，ｌ、時刻ｔ−１においてクラスタｋに所属していたオブジェクトが次の時刻ｔにどのクラスタに所属しやすいかを示す時間変化情報πｔ，ｋ、および時刻ｔにおいてオブジェクトｉが所属するクラスタを示す所属クラスタｚｔ，ｉと、ハイパーパラメータα０，κと、クラスタ数Ｋとを記憶する記憶部と、前記記憶部から前記所属クラスタｚｔ，ｉ、前記時間変化情報πｔ，ｋ、前記混合比β、前記ハイパーパラメータα０，κ、および前記クラスタ数Ｋとを取得し、所定期間ｔ＝１〜Ｔにおいて、ｚｔ−１，ｉ＝ｋかつｚｔ，ｉ＝ｌとなるオブジェクトの数をｍｔ，ｋ，ｌとしたとき、ディリクレ分布Dirichlet（α０β１＋ｍｔ，ｋ，ｌ，・，α０βｋ＋ｍｔ，ｋ，ｌ＋κ，・，α０βＫ＋ｍｔ，ｋ，Ｋ，α０（１−ΣＫｋ＝１βｋ））からサンプリングして前記時間変化情報πｔ，ｋを算出し、当該算出した時間変化情報πｔ，ｋによって前記記憶部に記憶してある時間変化情報πｔ，ｋを更新し記憶する時間変化情報推定部と、前記記憶部から前記所属クラスタｚｔ，ｉ、前記時間変化情報πｔ，ｋ、前記クラスタ間関連度ηｋ，ｌ、前記混合比β、所定期間ｔ＝１〜Ｔの前記関係データｘｔ，ｋ，ｌ、および前記クラスタ数Ｋを取得し、取得した前記混合比β、前記時間変化情報πｔ，ｋ、前記クラスタ間関連度ηｋ，ｌ、前記関係データｘｔ，ｋ，ｌ、およびクラスタ数Ｋを、式（１)、式（２）、式（３）、式（４）に適用して、式（１）からｕｔ，ｊをサンプリングし、メッセージ変数ｐｔ，ｉ，ｋを式（３）で定義したとき、ｔ＝１からｔ＝Ｔまで順番に式（２）を用いてメッセージ変数ｐｔ，ｉ，ｋを算出し、当該算出したメッセージ変数ｐｔ，ｉ，ｋに対して、ｔ＝Ｔからｔ＝１まで順番に式（４）を用いてメッセージ変数ｐｔ，ｉ，ｋを算出し、ｐ（ｚｔ，ｉ＝ｋ｜ｚｔ−１，ｉ＝ｌ）が０とならない場合の所属クラスタｚｔ，ｉを算出し、当該算出した所属クラスタｚｔ，ｉによって前記記憶部に記憶してある所属クラスタｚｔ，ｉを更新し記憶する所属クラスタ推定部と、前記クラスタ間関連度ηｋ，ｌを前記無限関係モデルによって算出し、前記記憶部に記憶するクラスタ間関連度推定部と、前記混合比βを前記無限関係モデルによって算出し、前記記憶部に記憶する混合比推定部と、前記混合比推定部、前記クラスタ間関連度推定部、前記時間変化情報推定部、および前記所属クラスタ推定部における演算を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す終了判定部と、を備えることを特徴とする。
【数１】
2011048583.tif 000003
【０００８】
また、本発明は、複数のオブジェクトの間の関係の有無を表す関係データを用いて、オブジェクトのクラスタリングを実行するクラスタリング装置において用いられるクラスタリング方法であって、前記クラスタリング装置が、前記関係データをクラスタリングする関数の一つである無限関係モデル（ＩＲＭ、Infinite Relational Model）において算出される混合比βおよびクラスタｋとクラスタｌとの間の関係の強さを示すクラスタ間関連度ηｋ，ｌと、クラスタｋとクラスタｌとの間の関連の有無を示す関係データを所定の時間間隔で観測した関係データｘｔ，ｋ，ｌ、時刻ｔ−１においてクラスタｋに所属していたオブジェクトが次の時刻ｔにどのクラスタに所属しやすいかを示す時間変化情報πｔ，ｋ、および時刻ｔにおいてオブジェクトｉが所属するクラスタを示す所属クラスタｚｔ，ｉと、ハイパーパラメータα０，κと、クラスタ数Ｋとを記憶する記憶部と処理部とを備え、前記処理部が、前記記憶部から前記所属クラスタｚｔ，ｉ、前記時間変化情報πｔ，ｋ、前記混合比β、前記ハイパーパラメータα０，κ、および前記クラスタ数Ｋとを取得し、所定期間ｔ＝１〜Ｔにおいて、ｚｔ−１，ｉ＝ｋかつｚｔ，ｉ＝ｌとなるオブジェクトの数をｍｔ，ｋ，ｌとしたとき、ディリクレ分布Dirichlet（α０β１＋ｍｔ，ｋ，ｌ，・，α０βｋ＋ｍｔ，ｋ，ｌ＋κ，・，α０βＫ＋ｍｔ，ｋ，Ｋ，α０（１−ΣＫｋ＝１βｋ））からサンプリングして前記時間変化情報πｔ，ｋを算出し、当該算出した時間変化情報πｔ，ｋによって前記記憶部に記憶してある時間変化情報πｔ，ｋを更新し記憶する時間変化情報推定ステップと、前記記憶部から前記所属クラスタｚｔ，ｉ、前記時間変化情報πｔ，ｋ、前記クラスタ間関連度ηｋ，ｌ、前記混合比β、所定期間ｔ＝１〜Ｔの前記関係データｘｔ，ｋ，ｌ、および前記クラスタ数Ｋを取得し、取得した前記混合比β、前記時間変化情報πｔ，ｋ、前記クラスタ間関連度ηｋ，ｌ、前記関係データｘｔ，ｋ，ｌ、およびクラスタ数Ｋを、前記式（１)、前記式（２）、前記式（３）、前記式（４）に適用して、前記式（１）からｕｔ，ｊをサンプリングし、メッセージ変数ｐｔ，ｉ，ｋを前記式（３）で定義したとき、ｔ＝１からｔ＝Ｔまで順番に前記式（２）を用いてメッセージ変数ｐｔ，ｉ，ｋを算出し、当該算出したメッセージ変数ｐｔ，ｉ，ｋに対して、ｔ＝Ｔからｔ＝１まで順番に前記式（４）を用いてメッセージ変数ｐｔ，ｉ，ｋを算出し、ｐ（ｚｔ，ｉ＝ｋ｜ｚｔ−１，ｉ＝ｌ）が０とならない場合の所属クラスタｚｔ，ｉを算出し、当該算出した所属クラスタｚｔ，ｉによって前記記憶部に記憶してある所属クラスタｚｔ，ｉを更新し記憶する所属クラスタ推定ステップと、前記クラスタ間関連度ηｋ，ｌを前記無限関係モデルによって算出し、前記記憶部に記憶するクラスタ間関連度推定ステップと、前記混合比βを前記無限関係モデルによって算出し、前記記憶部に記憶する混合比推定ステップと、前記混合比推定ステップ、前記クラスタ間関連度推定ステップ、前記時間変化情報推定ステップ、および前記所属クラスタ推定ステップにおける演算を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返し演算させる終了判定ステップと、を実行することを特徴とする。
【０００９】
このような構成によれば、公知のＩＲＭを拡張して、クラスタの時間変化を表す時間変化情報πｔ，ｋをパラメータとして導入し、所属クラスタｚｔ，ｉの推定方法も時間に従って考慮可能なようにモデルを構成したことによって、クラスタの時間変化を推定することが可能となる。
【００１０】
前記クラスタリング装置が、さらに、ハイパーパラメータγ，ξ，Ψを記憶する前記記憶部を備え、前記クラスタ間関連度推定部は、前記記憶部から前記所属クラスタｚｔ，ｉ、前記クラスタ間関連度ηｋ，ｌ、前記ハイパーパラメータξ，Ψ、前記クラスタ数Ｋ、および所定期間ｔ＝１〜Ｔの前記関係データｘｔ，ｋ，ｌを取得して、ｚｔ，ｉ＝ｋかつｚｔ，ｊ＝ｌとなる（ｔ，ｉ，ｊ）の組の数をＮｋ，ｌとし、当該Ｎｋ，ｌの中の前記関係データｘｔ，ｋ，ｌが関係有りを示す数をｎｋ，ｌとしたとき、ベータ分布Beta（ξ＋ηｋ，ｌ，Ψ＋Ｎｋ，ｌ−ｎｋ，ｌ）からサンプリングして前記クラスタ間関連度ηｋ，ｌを算出し、当該算出したクラスタ間関連度ηｋ，ｌによって前記記憶部に記憶してあるクラスタ間関連度ηｋ，ｌを更新し記憶し、前記混合比推定部は、前記記憶部から前記所属クラスタｚｔ，ｉ、前記混合比β、および前記ハイパーパラメータγ，α０，κを取得し、取得した前記所属クラスタｚｔ，ｉ、前記混合比β、および前記ハイパーパラメータα０，κを、式（５）および式（６）に適用して、それぞれ補助変数Ｒｔ，ｋ，ｌおよび補助変数Ｏｔ，ｋを算出し、前記算出した補助変数Ｒｔ，ｋ，ｌおよび補助変数Ｏｔ，ｋを式（７）に適用して補助変数＾Ｒｔ，ｋ，ｌを算出し、ディリクレ分布Dirichlet（Σｔ，ｋ＾Ｒｔ，ｋ，１，Σｔ，ｋ＾Ｒｔ，ｋ，２，・，Σｔ，ｋ＾Ｒｔ，ｋ，γ）からサンプリングして前記混合比βを算出し、当該算出した混合比βによって前記記憶部に記憶してある混合比βを更新し記憶することを特徴とする。
【数２】
2011048583.tif 000004
【００１１】
また、本発明は、前記クラスタリング装置が、さらに、ハイパーパラメータγ，ξ，Ψを記憶する前記記憶部を備え、前記処理部が、前記クラスタ間関連度推定ステップにおいて、前記記憶部から前記所属クラスタｚｔ，ｉ、前記クラスタ間関連度ηｋ，ｌ、前記ハイパーパラメータξ，Ψ、前記クラスタ数Ｋ、および所定期間ｔ＝１〜Ｔの前記関係データｘｔ，ｋ，ｌを取得して、ｚｔ，ｉ＝ｋかつｚｔ，ｊ＝ｌとなる（ｔ，ｉ，ｊ）の組の数をＮｋ，ｌとし、当該Ｎｋ，ｌの中の前記関係データｘｔ，ｋ，ｌが関係有りを示す数をｎｋ，ｌとしたとき、ベータ分布Beta（ξ＋ηｋ，ｌ，Ψ＋Ｎｋ，ｌ−ｎｋ，ｌ）からサンプリングして前記クラスタ間関連度ηｋ，ｌを算出し、当該算出したクラスタ間関連度ηｋ，ｌによって前記記憶部に記憶してあるクラスタ間関連度ηｋ，ｌを更新し記憶し、前記混合比推定ステップにおいて、前記記憶部から前記所属クラスタｚｔ，ｉ、前記混合比β、および前記ハイパーパラメータγ，α０，κを取得し、取得した前記所属クラスタｚｔ，ｉ、前記混合比β、および前記ハイパーパラメータα０，κを、前記式（５）および前記式（６）に適用して、それぞれ補助変数Ｒｔ，ｋ，ｌおよび補助変数Ｏｔ，ｋを算出し、前記算出した補助変数Ｒｔ，ｋ，ｌおよび補助変数Ｏｔ，ｋを前記式（７）に適用して補助変数＾Ｒｔ，ｋ，ｌを算出し、ディリクレ分布Dirichlet（Σｔ，ｋ＾Ｒｔ，ｋ，１，Σｔ，ｋ＾Ｒｔ，ｋ，２，・，Σｔ，ｋ＾Ｒｔ，ｋ，γ）からサンプリングして前記混合比βを算出し、当該算出した混合比βによって前記記憶部に記憶してある混合比βを更新し記憶することを特徴とする。
【００１２】
このような構成によれば、時間期間に対応した混合比βおよびクラスタ間関連度ηｋ，ｌを用いて、時間変化情報πｔ，ｋおよび所属クラスタｚｔ，ｉの推定方法を行うことが可能なようにモデルを構成したことによって、クラスタの時間変化を推定することが可能となる。
【００１３】
本発明は、前記クラスタリング装置が、演算結果を表示する表示装置と接続され、前記終了判定部が、前記終了条件を満足した場合、前記混合比β、前記クラスタ間関連度ηｋ，ｌ、前記時間変化情報πｔ，ｋ、および前記所属クラスタｚｔ，ｉのいずれか一つまたはいずれかの組み合わせを前記表示装置に出力することを特徴とする。
【００１４】
また、本発明は、前記クラスタリング装置は、演算結果を表示する表示装置と接続され、前記処理部が、前記終了条件を満足した場合、前記混合比β、前記クラスタ間関連度ηｋ，ｌ、前記時間変化情報πｔ，ｋ、および前記所属クラスタｚｔ，ｉのいずれか一つまたはいずれかの組み合わせを前記表示装置に出力することを特徴とする
【００１５】
このような構成によれば、時間変化情報πｔ，ｋと他のパラメータとを合わせて表示装置に出力することができる。
【００１６】
本発明は、前記クラスタリング方法を、コンピュータとしてのクラスタリング装置に実行させるためのプログラムとした。
【００１７】
このようなプログラムをインストールされたコンピュータは、このプログラムに基づいた機能を実現することができる。
【発明の効果】
【００１８】
本発明によれば、複数のオブジェクトの間の関係の有無を示す関係データの時系列から、オブジェクトのクラスタの存在を推定し、クラスタの時間変化を推定する技術を提供することができる。
【図面の簡単な説明】
【００１９】
【図１】本実施形態におけるｄＩＲＭの処理の流れを示す図である。
【図２】クラスタリング装置の構成の一例を示す図である。
【図３】所属クラスタ推定部の処理フローを示す図である。
【図４】時間変化情報推定部の処理フローを示す図である。
【図５】クラスタ間関連度推定部の処理フローを示す図である。
【図６】混合比推定部の処理フローを示す図である。
【図７】（ａ）は、ｄＩＲＭに対するｒａｎｄ  ｉｎｄｅｘの結果を示し、（ｂ）は、ＩＲＭに対するｒａｎｄ  ｉｎｄｅｘの結果を示す図である。
【図８】（ａ）オブジェクト間の関係の一例を示す図であり、（ｂ）ＩＲＭ適用後のオブジェクト間の関係の一例を示す図である。
【発明を実施するための形態】
【００２０】
次に、本発明を実施するための形態（以降「本実施形態」と称す）について、適宜図面を参照しながら詳細に説明する。
【００２１】
≪公知のＩＲＭの概要≫
初めに、ＩＲＭの概要について、図８を用いて説明する。図８（ａ）は、Ｎ種類のオブジェクトからなるドメインＤ＝｛１，２, ・・・, Ｎ｝上の二値の二項関係（２つのオブジェクトの間の関係）を示す図である。縦方向には一方のオブジェクトのインデックスｉ（Ｎ＝９）を表し、横方向には他方のオブジェクトのインデックスｊ（Ｎ＝９）を表している。図８（ａ）では、２つのオブジェクト間に関連が有る場合には、その升目を黒く塗りつぶして示している。一般的に、コンピュータによってこの関係の有無を数値で表す場合、演算を行いやすくするために、オブジェクトのインデックスｉ，ｊ 間に関係が無いときにはｘｉ，ｊ＝０と表し、オブジェクトのインデックスｉ，ｊ 間に関係が有る場合にはｘｉ，ｊ＝１と表すものとする。このデータｘｉ，ｊは、関係データと呼ばれる。
【００２２】
そして、図８（ａ）に示す関係データをＩＲＭに適用すると、図８（ｂ）に示すように、オブジェクトのインデックスｉ，ｊ間に関係の有るインデックスとそうでないインデックスとがそれぞれグループ化され、太い実線で区切られたクラスタが生成される。図８（ｂ）では、クラスタ数が自動的に３つに決定され、クラスタを識別する変数（以降、クラスタ番号とも称す。）ｋが１〜３と設定される。このとき、オブジェクトのインデックスｉが所属するクラスタをｚｉ＝ｋと表現する。
【００２３】
ＩＲＭの数学的モデルを、式（８）〜式（11）に示す。
【数３】
2011048583.tif 000005
【００２４】
ここで、β，ｚｉ，ηｋ，ｌ，ｘｉ，ｊについては後記する。また、「〜」は確率分布からサンプリングすることを表す。Stick()は、ＤＰＭにおいて用いられる分布である。Multinomial()、Beta()、およびBernoulli()は、それぞれ、多項分布、ベータ分布、ベルヌイ分布を表す。また、γ，ξ，Ψは事前に設定するハイパーパラメータである。
【００２５】
まず、式（８）では、無限次元のクラスタ混合比ベクトルβ（以降、単に混合比と称する。）を生成する。式（８）は具体的には式（12）のように計算される。
【数４】
2011048583.tif 000006
【００２６】
ここで、式（12）における矢印⇔は、左辺のｋ番目のクラスタが右辺で表されることを示す記号である。そして、式（12）の右辺の混合比βｋは、ｋ番目のクラスタにデータが所属している確率を表している。定義によって、Σｋβｋ＝１であれば、βは無限個のクラスタの混合比として用いることができる。
【００２７】
式（９）では、混合比βを用いて、オブジェクトのインデックスｉが所属するクラスタｚｉ＝ｋを多項分布（Multinomial()）からサンプリングする。式（10）および式（11）は、クラスタリング処理後のクラスタｚｉ＝ｋが与えられたのちに、実際に観測された関係データｘｉ，ｊを生成する過程を表す。式（10）では、クラスタｋ，ｌ間の関係の強さを示すクラスタ間関連度ηｋ，ｌをベータ分布（Beta()）からサンプリングする。このサンプリングされた値は、図８（ｂ）において、（ｋ，ｌ）で表されるクラスタ内の升目に黒（＝１）を含む確率を表す。具体的には、クラスタ（ｋ＝１，ｌ＝３）のクラスタ間関連度η１，３は０．１７（＝１／６）となる。また、個別の関係データｘｉ，ｊの値は、各オブジェクトｉ，ｊの所属する各クラスタｚｉ，ｚｊによって規定されるブロックの関係の強さを表すクラスタ間関連度ηｚｉ，ｚｊに基づいて、ベルヌイ分布（Bernoulli()）からサンプリングされる。
【００２８】
≪時間を含む時系列関係データをＩＲＭに適用する際の問題≫
次に、前記したＩＲＭの数学的モデルに、時間を含む時系列関係データを適用する際の問題の一例を説明する。今、時刻データを含んだ関係データ時系列Ｘ＝｛ｘｔ，ｉ，ｊ∈｛０，１｝,１≦ｔ≦Ｔ｝が与えられたとする。ここで、ｘｔ，ｉ，ｊ＝１は時刻ｔにおいてオブジェクトのインデックスｉ，ｊ間に関係が有ることを示す。なお、ｘｔ，ｉ，ｊは、時刻ｔと異なる時刻におけるオブジェクト間の関係を定義されないものとする。
【００２９】
まず、ＩＲＭの数学的モデルの式（８），式（９），式（10），式（11）には時刻インデックスｔが含まれていないため、時刻を含んだ関係データを適用することはできない。この問題を解決するためには、いくつかの方法が考えられる。最も単純には、時刻データを含んだ関係データ時系列Ｘから時刻データを含まない関係データ~Ｘを生成して、前記したＩＲＭの数学的モデルに適用する。例えば、関係データ~ｘｉ，ｊを、式（13）のように表す。
【００３０】
【数５】
2011048583.tif 000007
【００３１】
ここで、σは例えばσ＝０．５等、予め設定しておく閾値である。関係データ~ｘｉ，ｊにおいて、時間情報が欠落する以上当然ではあるが、クラスタリングの時間変化等は完全に無視されることになる。この場合、時刻のインデックスを含んだ関係データ~ｘｉ，ｊをＩＲＭに適用しようとすると、例えば、式（14）〜式（17）の様に表される。
【００３２】
【数６】
2011048583.tif 000008
【００３３】
ここで、ｚｔ，ｉ＝ｋは時刻ｔにおけるオブジェクトのインデックスｉの所属クラスタを表す。このｚｔ，ｉ＝ｋは、ＩＲＭにおいて、オブジェクトの所属クラスタに時刻ｔのインデックスを追加しただけであるので、時間の情報が活かされていない。そして、前記式（14）〜式（17）に示す数学的モデルでは、どの時刻のどのオブジェクトの所属クラスタインデックスも独立に同じ分布からサンプリングされることになる。したがって、時刻ｔ−１と時刻ｔではクラスタリングの結果に相関関係がある、といったクラスタリングの時間ダイナミクスを一切表現することができない。
【００３４】
≪本実施形態におけるクラスタリングモデル（ｄＩＲＭ）≫
次に、ＩＲＭを拡張して、時間変化する関係データを考慮可能とする、新しいクラスタリングモデルdynamic Infinite Relational Model（ｄＩＲＭ）について説明する。まず、時間変化を表すために、オブジェクトのクラスタリングに隠れマルコフモデル（ＨＭＭ：Hidden Markov Model）に似た構造を導入し、時刻ｔ−１でのクラスタリングと時刻ｔでのクラスタリングとの間に依存関係を持たせる。これによって、例えば、時刻ｔ−１においてクラスタ１に所属したオブジェクトは次の時刻ｔにおいてクラスタ１とクラスタ２に所属しやすい、等のように、関係データ時系列の様々な側面をモデル化することが可能となる。ＨＭＭの詳細は、例えば、「L.R.Rabiner,“A Tutorial on Hidden Markov Models and Selected applications in Speech Recognition”, Proceedings of the IEEE, Vol.77, No.2, p. 257-286, 1989」に記載されている。
【００３５】
ＨＭＭは、非常に幅広い分野で利用されている時系列モデルである。ＨＭＭでは、各時刻ｔにおいて、観測不能な隠れ状態ｓｔと観測データｙｔとを定義する。ＨＭＭは、隠れ状態間の時間遷移確率ｐ（ｓｔ｜ｓｔ−１）と隠れ状態からの観測モデルｐ（ｙｔ｜ｓｔ）とで特徴付けられる。時間遷移確率ｐ（ｓｔ｜ｓｔ−１） は、時刻ｔにおける隠れ状態ｓｔの値が時刻ｔ−１における隠れ状態の値ｓｔ−１に依存して確率的に決定されることを表す。観測モデルｐ（ｙｔ｜ｓｔ）は、時刻ｔにおける観測データが同時刻の隠れ状態量に依存して決定されることを表す。
【００３６】
このように、ｄＩＲＭでは、ＨＭＭを応用し、オブジェクトのインデックスｉが時刻ｔにおいて所属するクラスタｚｔ，ｉを隠れ状態量と定義することで、クラスタリングの時間ダイナミクスをモデル化している。
【００３７】
まず、公知のＨＭＭと同様に、時刻ｔ−１において、あるクラスタに所属していたオブジェクトが、次の時刻ｔにどのクラスタに所属しやすいか、という遷移確率を表すパラメータを用意する。また、時間ごとにクラスタ間の遷移確率自体が変化することも考えられることから、遷移確率パラメータは時刻ごとに分けて考える。なお、ｄＩＲＭではクラスタ数、すなわちＨＭＭの隠れ状態の取りうる値が不定（無限）であるため、単純なＨＭＭは利用できない。そこで、ｄＩＲＭでは、公知のinfinite  ＨＭＭで示される隠れ状態数が各時刻で１つ(ｓｔのみ)であることを拡張して、隠れ状態数が各時刻でＮ個のｚｔ，ｉとなるように設定する。具体的なモデルは次の様になる。
【００３８】
【数７】
2011048583.tif 000009
【００３９】
ここで、γ，α０，κ，ξ，Ψは事前に設定されるハイパーパラメータである。式（18）、式（21）、および 式（22）はＩＲＭの場合の式（14）、式（16）、および式（17）と同じである。式（19）は、新たに導入されたパラメータであり、式（20）は、ＩＲＭの場合の式（15）と異なる。式（19）に示すπｔ，ｋは、時刻ｔ−１においてクラスタｋに所属していたオブジェクトが、時刻ｔにどのクラスタに所属しやすいかを表すパラメータであり時間変化情報と称する。この時間変化情報πｔ，ｋは、公知のＨＭＭにおける状態間の遷移確率を時間依存で変化するよう拡張したことに対応する。また、式（20）はオブジェクトｉの時刻ｔにおける状態ｚｔ，ｉが時刻ｔ−１における状態に依存することを表している点で、公知のＨＭＭモデルのアイデアを利用している。式（19）中のＤＰ()はDirichlet Processを表している。このＤＰ()は、直感的には無限次元のDirichlet分布と考えることができる。したがって、無限次元のパラメータα０βから無限次元のベクトルπを生成することになる。この演算を簡略化するために、有限次元Ｌで近似した場合は式（23）の様なDirichlet分布となる。
【００４０】
【数８】
2011048583.tif 000010
【００４１】
ここで、β（Ｌ）はＬ次元ベクトルである。また、式（23）中の矢印→は、近似を表す。より詳しい数学的解説は「上田修功，山田武士,“ノンパラメトリックベイズモデル”,応用数理, Vol.8, No.3, p.16-214, 2007」に記載されている。
【００４２】
また、式（19）のκ（ただし、κ＞０）は、グローバルな混合比α０×βのｋ番目の要素にκを加えて、サンプリングされるπｔ，ｋのｋ番目の要素の値が大きくなりやすくするものである。なお、κについて詳細は、「E.B.Fox, E.B.Sudderth, M.I.Jordan and A.S.Willsky, “An HDPHMM for Systems with State Persistence”, Proceedings of the International Conference on Machine Learning (ICML), 2008」に記載されている。また、式（19）のδｋは、クラスタ番号ｋのときに「１」の値となり、それ以外では「０」となるデルタ関数である。つまり、式（19）でのκの作用は、ＨＭＭによる隠れ状態ｓｔ（ここではｚｔ，ｉ）の決定は確率的なので、例えばｓ１＝ｓ２＝・・・＝ｓt＝ｋという値が正しい場合にも、確率的に一部の値がｋ以外の値をとる可能性が高くなることを避けることである。
【００４３】
≪ｄＩＲＭの処理の流れ≫
図１を用いて、まず、ｄＩＲＭの処理の大まかな流れを説明し、個々の処理の詳細については後記する。なお、図１に示す処理は、クラスタリング装置８０（図２参照）が実行する。
【００４４】
ステップＳ１０１では、処理に用いるハイパーパラメータγ，α０，κ，ξ，Ψ、および観測値（関係データ時系列）Ｘ＝｛ｘｔ，ｉ，ｊ｝を取得する。ステップＳ１０２では、所属クラスタｚｔ，ｉ、時間変化情報πｔ，ｋ、クラスタ間関連度ηｋ，ｌ、混合比β、クラスタ数Ｋを初期化する。クラスタ数Ｋは、クラスタを表す変数ｋの取りうる数である。ステップＳ１０３では、繰り返し変数ｉｔｒ＝１と設定する。そして、ステップＳ１０４〜Ｓ１０８まで繰り返し演算を実行する。
【００４５】
ステップＳ１０４では、所属クラスタ（ｚ）推定部２３の処理（図３参照）を用いて、、所属クラスタｚｔ，ｉをサンプリングして、そのサンプリングした所属クラスタｚｔ，ｉによって記憶部３０（図６参照）に記憶されている所属クラスタｚｔ，ｉを更新し保存（記憶）する。ステップＳ１０５では、時間変化情報（π）推定部２４の処理（図４参照）を用いて、時間変化情報πｔ，ｋをサンプリングして、そのサンプリングした時間変化情報πｔ，ｋによって記憶部３０に記憶されている時間変化情報πｔ，ｋを更新し保存（記憶）する。ステップＳ１０６では、クラスタ間関連度（η）推定部２５の処理（図５参照）を用いて、クラスタ間関連度ηｋ，ｌをサンプリングして、そのサンプリングしたクラスタ間関連度ηｋ，ｌによって記憶部３０に記憶されているクラスタ間関連度ηｋ，ｌを更新し保存（記憶）する。ステップＳ１０７では、混合比（β）推定部２６の処理（図６参照）を用いて、混合比βをサンプリングして、そのサンプリングした混合比βによって記憶部３０に記憶されている混合比βを更新し保存（記憶）する。ステップＳ１０８では、繰り返し変数ｉｔｒを「１」加算してカウント数を増加する。
【００４６】
ステップＳ１０９では、終了条件を満足したか否かを判定する。終了条件は、例えば、予め決めておいた所定の繰り返し回数になったことである。そして、終了条件を満足していない場合（ステップＳ１０９でＮｏ）、処理はステップＳ１０４へ戻る。また、終了条件を満足する場合（ステップＳ１０９でＹｅｓ）、ステップＳ１１０では、演算結果を出力し、処理を終了する。
【００４７】
なお、前記した終了条件の代わりとして、終了条件を、１つ前のサンプリング処理におけるパラメータの値との差分の絶対値が、予め設定しておいた所定の閾値以下になったこと、としても構わない。また、演算結果の出力においては、所属クラスタｚｔ，ｉ，時間変化情報πｔ，ｋ，クラスタ間関連度ηｋ，ｌ，混合比βのパラメータ中から、必要なパラメータのデータのみを出力する。また、ステップＳ１０１においてハイパーパラメータγ，κ，α０，ξ，Ψを事前に与える代わりに、繰り返し演算のループの中で、他のパラメータと同時に推定しても良い（例えば、「Y.W.Teh, M.I.Jordan, M.J.Beal and D.M.Blei,“Hierarchical Dirichlet Process”, Journal of the American Statistical Association, Vol.101, No.476, p.1566-1581, 2006」参照）。
【００４８】
次に、図１の各ステップの処理を実行する、クラスタリング装置８０（後記する図２参照）の構成および処理の詳細について説明する。まず、処理の前提について、先に説明する。例えば、所属クラスタｚｔ，ｉ、時間変化情報πｔ，ｋ、クラスタ間関連度ηｋ，ｌ、混合比βを推定するために、ｂｅａｍ  ｓａｍｐｌｉｎｇを用いる場合を示す。この理由は、ｂｅａｍ  ｓａｍｐｌｉｎｇが、無限個存在するクラスタを有限個に打ち切って推定するので、無限次元ベクトルとして扱わなければならなかったβやπｔ，ｋを有限次元ベクトルで扱えるためである。そして、ｂｅａｍ  ｓａｍｐｌｉｎｇでは、有限個に打ち切って推定したクラスタの数は、サンプリングを繰り返すことで、理論上無限次元ベクトルを正当に扱った場合と等価になる。なお、ｂｅａｍ  ｓａｍｐｌｉｎｇの代わりにＧｉｂｂｓサンプリングによる方法を用いることもできるが、ｂｅａｍ  ｓａｍｐｌｉｎｇは、前記したように有限個に打ち切って演算を実行するので、Ｇｉｂｂｓサンプリングよりも高速に処理を行うことができるという利点がある。
【００４９】
また、以下の説明では、２つのドメインのオブジェクトの数は同じ（インデックスｉ，ｊがどちらも１，２，・・・，Ｎ）ものとするが、これらは同じでなくても良い。また、本実施形態では、クラスタを表す変数ｋを正整数と仮定して説明を行うが、変数ｋはクラスタリングの各グループの名前を識別できることが目的であるため、変数ｋの値そのものには意味はなく、正整数のかわりにａ，ｂ，ｃ等の記号（シンボル）や文字列で表現しても良い。
【００５０】
≪クラスタリング装置の構成≫
次に、クラスタリング装置８０の構成について、図２を用いて説明する。図２に示すように、クラスタリング装置８０は、処理部２０と記憶部３０とを備え、処理部２０がコンピュータにおけるＣＰＵ（Central Processing Unit）に相当し、記憶部３０が主記憶装置やＨＤＤ（Hard Disc Drive）やＵＳＢ（Universal Serial Bus）メモリ等に相当する。また、クラスタリング装置８０には、クラスタリング装置８０へデータを入力する入力装置１０、およびクラスタリング装置８０の演算結果を表示する表示装置４０が接続可能である。入力装置１０は、例えば、キーボードやマウス等であり、ユーザによるデータ入力や処理操作の指示に用いられる。表示装置４０は、例えば、ディスプレイ等であって、ユーザによる演算結果の確認に用いられる。なお、入力装置１０および表示装置４０は、必須ではなく、クラスタリング装置８０に外部接続されるＵＳＢメモリ等の記憶部（図示せず）を介して、該記憶部に記憶されているデータを入力とし、演算結果を該記憶部に出力するようにしても構わない。ただし、本実施形態では、クラスタリング装置８０に入力装置１０および表示装置４０が接続されているケースについて説明する。
【００５１】
クラスタリング装置８０の処理部２０は、入力データの初期化を実行する初期設定部２１と、繰り返し演算を実行する推定演算部２２とを備える。推定演算部２２は、機能として、所属クラスタ（ｚ）推定部２３、時間変化情報（π）推定部２４、クラスタ間関連度（η）推定部２５、混合比（β）推定部２６、および終了判定部２８を備える。初期設定部２１は、図１に示すステップＳ１０１およびＳ１０２を実行する。所属クラスタ（ｚ）推定部２３、時間変化情報（π）推定部２４、クラスタ間関連度（η）推定部２５、および混合比（β）推定部２６は、それぞれ、図１に示すステップＳ１０４、Ｓ１０５、Ｓ１０６、およびＳ１０７を実行する。そして、終了判定部２８は、図１に示すステップＳ１０９およびＳ１１０を実行する。なお、図２中に破線で示したハイパーパラメータ算出部２７は、ハイパーパラメータを事前に与える代わりに、繰り返し演算のループの中で推定する場合に用いられる。また、各部２３〜２７の演算する順番は、図示している順に限られない。
【００５２】
記憶部３０は、処理部２０の演算に用いられる各変数として、推定値３１、補助変数３２、メッセージ変数３３、カウント変数３４、クラスタ数３６、および観測値３７を記憶する。また、記憶部３０には、処理部２０によって実行されるアプリケーションプログラムが記憶されている。
【００５３】
次に、図２に示すクラスタリング装置８０の処理部２０の各部２１〜２８の処理の詳細について説明する。
【００５４】
（初期設定部）
初期設定部２１は、入力装置１０から、関係データ時系列Ｘ＝｛ｘｔ，ｉ，ｊ｝（ｔ＝１，２，・，Ｔ，ｉ＝１，２，・，Ｎ，ｊ＝１，２，・，Ｎ）とハイパーパラメータγ，α０，κ，ξ，Ψとを取得し（図１のステップＳ１０１）、それぞれ記憶部３０の観測値３７とハイパーパラメータ３５として記憶する。また、初期設定部２１は、式（18）〜式（22）において定義されたパラメータβ，ｚｔ，ｉ，πｔ，ｋ，ηｋ，ｌおよびクラスタ数Ｋの初期値を設定し、記憶部３０の推定値３１とクラスタ数３６として記憶する（図１のＳ１０２)。初期値の設定は次の様に行う。クラスタ数Ｋの初期値は、ランダムな正整数を設定する。すべてのｚｔ，ｉについては、１からクラスタ数Ｋのうちのいずれかの整数値を設定する。これは、完全ランダムに設定すれば良い。Ｌ次元のベクトルであるπｔ，ｋについては、ベクトルの要素の和が１であり、かつ、すべての要素が非負であるという条件を満たすように、ランダムな値を設定する。ηｋ，ｌについては、０≦ηｋ，ｌ≦１なる実数値をランダムに割り当てる。Ｌ次元のベクトルであるβについては、ベクトルの要素の和が１であり、かつ、すべての要素が非負であるという条件を満たすように、ランダムな値を設定する。
【００５５】
（所属クラスタ（ｚ）推定部）
所属クラスタ（ｚ）推定部２３は、ｚｔ，ｉのサンプリングを行い、そのサンプリングした値によって記憶部３０に保存してあるｚｔ，ｉを更新し記憶する。所属クラスタ（ｚ）推定部３１の処理フローについて図３を用いて説明する。まず、ステップＳ３０１では、所属クラスタ（ｚ）推定部３１は、現在記憶部３０に保存してあるｚｔ，ｉ，πｔ，ｋ，ηｋ，ｌ，β，Ｋ，Ｘ＝｛ｘ１，１，１，ｘ１，１，２，・・・，ｘＴ，Ｎ，Ｎ｝を取得する。次に、ステップＳ３０２では、ｔ＝１〜Ｔ，ｉ＝１〜Ｎについて、式（１）を用いて、補助変数ｕｔ，ｉを次のように算出する。
【００５６】
【数９】
2011048583.tif 000011
【００５７】
ここで、Uniform()は、一様分布を表す。つまり、式（１）は、ｕｔ，ｉ（t＝１〜Ｔ，ｉ＝１〜Ｎ）が、一様分布からサンプリングされることを表す。
【００５８】
続いて、ステップＳ３０３〜Ｓ３０９では、ｔ＝１からｔ＝Ｔまですべてのｉ＝1〜Ｎおよび必要なｋについてメッセージ変数ｐｔ，ｉ，ｋを算出する。このメッセージ変数の算出は、ｔ＝１からｔ＝Ｔまで順番に行うので、ｆｏｒｗａｒｄ  ｆｉｌｔｅｒｉｎｇと呼ぶ。具体的には、ステップＳ３０３ではｔ＝１と設定し、ステップＳ３０４ではｉ＝１と設定し、ステップＳ３０５では式（２）を用いてメッセージ変数ｐｔ，ｉ，ｋを算出し、ステップＳ３０６ではｉを「１」増加し、ステップＳ３０７でｉ＞Ｎの場合はステップＳ３０８でｔを「１」増加し、ステップＳ３０９でｔ＞Ｔの場合にステップＳ３１０へ処理を移行する。なお、ステップＳ３０７でｉ≦Ｎの場合はステップＳ３０５へ戻り、ステップＳ３０９でｔ≦Ｔの場合はステップＳ３０４へ戻る。
【００５９】
なお、ステップＳ３０５では、メッセージ変数ｐｔ，ｉ，ｋは、式（２）によって算出される。ただし、メッセージ変数は式（３）を満たすような変数である。
【数１０】
2011048583.tif 000012
【００６０】
ここで、式（２）は、ｚｔ，ｉに関するすべての時系列関係データの尤度を既存のメッセージ変数に乗算することを表す。なお、ｋに関してはｕｔ＋１，ｉ＜πｔ＋１，ｋ，ｚｔ＋１，ｉとなるｋ、すなわち有限個のｋに関してのみ計算し、その他のｋに関してはｐｔ，ｉ，ｋ＝０とする。つまり、ｚｔ，ｉ＝ｋのサンプリングにおいて、ｕｔ，ｉ＞πｔ，ｚｔ―１，ｉ，ｋとなるクラスタを表すクラスタ番号ｋにはオブジェクトのインデックスｉをアサイン（関連付け）できなくすることによって、ｚｔ，ｉの取りうるクラスタ番号の数を無限個から有限個に削減することができる。また、ｔ＝１の場合には、右辺最終項を無視して演算する。
【００６１】
また、式（３）において、ｘ１：ｔ，・，・は、ｘｔ´，ｉ，ｊの全集合のうち、１≦ｔ´≦ｔ, １≦ｉ，ｊ≦Ｎを満たすすべての値の集合を表す。つまり、ｘ１：ｔ，・，・は、時刻ｔ＝１〜ｔのすべての関係データである。ｕ１：ｔ，・は、ｕ１：ｔ´，ｉのうち、１≦ｔ´≦ｔ, １≦ｉ≦Ｎを満たす集合を表す。なお、「・」はすべてのインデックスを示す。
【００６２】
次に、ステップＳ３１０〜Ｓ３１６では、ｔ＝Ｔ，Ｔ−１，・，１，ｉ＝１〜Ｎについて、ｂａｃｋｗａｒｄ  ｓａｍｐｌｉｎｇ によってｚｔ，ｉ＝ｋの値を算出する。具体的には、ステップＳ３１０ではｔ＝Ｔと設定し、ステップＳ３１１ではｉ＝１と設定し、ステップＳ３１２では式（４）を用いてｚｔ，ｉをサンプリングして更新保存し、ステップＳ３１３ではｉを「１」増加し、ステップＳ３１４でｉ＞Ｎの場合はステップＳ３１５でｔから「１」減少し、ステップＳ３１６でｔ＜１の場合はステップＳ３１７へ処理を移行する。なお、ステップＳ３１４でｉ≦Ｎの場合はステップＳ３１２へ戻り、ステップＳ３１６でｔ≧１の場合はステップＳ３１１へ戻る。
【００６３】
ステップＳ３１２では、メッセージ変数ｐｔ，ｉ，ｋは、式（４）を用いて算出される。
【数１１】
2011048583.tif 000013
【００６４】
ここで、右辺のＩ（）は、カッコ内の条件式が満たされれば１、そうでなければ０の値をとる関数である。したがって、ここでもｋに関してはｕｔ＋１，ｉ＜πｔ＋１，ｋ，ｚｔ＋１，ｉとなるｋ、すなわち有限個のｋに関してのみ計算すれば良い。式（４）を、ｔ＝Ｔからｔ＝１ まで、すべてのｉ，ｋについて演算してｚｔ，ｉ＝ｋをサンプリングする。このとき、式（４）の右辺が０となる場合は、そのときのｚｔ，ｉ＝ｋは無視される（サンプリングされない）。そして、そのサンプリングしたｚｔ，ｉによって記憶部３０に記憶されているｚｔ，ｉを更新し保存する。その結果として、ｚｔ，ｉのクラスタ番号の取りうるクラスタ数Ｋ個の値だけが選択される。このように、ｄＩＲＭでは、サンプリングによって選択されたｚｔ，ｉのクラスタ番号がＫ個に限定されるが、一般性を失うことはない。
【００６５】
ステップＳ３１７では、以上の処理によってサンプリングして保存されたすべてのｚｔ，ｉのクラスタ番号の取りうるクラスタ数をＫとして記憶部３０のＫを更新し保存する。次に、ステップＳ３１８では、記憶部３０に保存されているすべての変数について、取りうるクラスタ番号を１〜Ｋに更新し、保存する。
【００６６】
このステップＳ３１８の処理の具体例を、以下に示す。いま、時間ステップ数ｔ＝１でオブジェクトの数Ｎ＝５の場合に、記憶されているｚ１，ｉが次のようになっていたとする。
ｚ１，１＝１
ｚ１，２＝３
ｚ１，３＝４
ｚ１，４＝６
ｚ１，５＝６
【００６７】
この場合、使用されているクラスタ番号は１，３，４，６の４種類であるためＫ＝４となる。上記のクラスタ番号の範囲１〜６のうち、クラスタ番号２と５は使用されていないので、クラスタ番号のインデックス「３」を「２」に、「４」を「３」に、「６」を「４」に置き換えることにより、無駄な（使用されていない）クラスタの番号がないようにする。
ｚ１，１＝１  （更新しない）
ｚ１，２＝３ → ２  
ｚ１，３＝４ → ３
ｚ１，４＝６ → ４
ｚ１，５＝６ → ４
【００６８】
これにより、すべての変数ｚｔ，ｉについて、とりうるクラスタの番号が１〜４のいずれかになるように更新される。そして、このようにｚｔ，ｉを書き換えると、他のパラメータπ、η、βや補助変数のインデックスｋ，ｌにも影響を及ぼすので、これらのすべての変数について上記ｚｔ，ｉと同様にインデックスの書き換え処理（ステップＳ３１８）を行う。
【００６９】
（時間変化情報（π）推定部）
図２に示す時間変化情報（π）推定部２４は、πｔ，ｋのサンプリングを行い、そのサンプリングしたπｔ，ｋによって記憶部３０に保存してあるπｔ，ｋを更新し記憶する。時間変化情報（π）推定部２４の処理フローについて図４を用いて説明する。
【００７０】
まず、ステップＳ４０１では、記憶部３０に現在保存してあるｚｔ，ｉ，πｔ，ｋ，β，Ｋとハイパーパラメータα０，κとを取得する。ここで、ｔを固定したときにｚｔ−１，ｉ＝ｋかつｚｔ，ｉ＝ｌとなるオブジェクトの数をｍｔ，ｋ，ｌと表し、そのｍｔ，ｋ，ｌをカウント変数とする。ステップＳ４０２では、すべてのｔ＝１〜Ｔ、ｋ＝１〜Ｋ、ｌ＝１〜Ｋについて、ｍｔ，ｋ，ｌ＝０に初期化（設定）する。
【００７１】
次に、ステップＳ４０３〜Ｓ４１１によって、ｔ＝１〜Ｔ、ｉ＝１〜Ｎについてｍｔ，ｋ，ｌを算出する。具体的には、ステップＳ４０３ではｔ＝１と設定し、ステップＳ４０４ではｉ＝１と設定し、ステップＳ４０５ではｋ＝ｚｔ−１，ｉと設定し、ステップＳ４０６ではｌ＝ｚｔ，ｉと設定し、ステップＳ４０７ではｍｔ，ｋ，ｌを「１」増加し、ステップＳ４０８ではｉを「１」増加し、ステップＳ４０９ではｉ＞Ｎの場合はステップＳ４１０でｔを「１」増加し、ステップＳ４１１ではｔ＞Ｔの場合にステップＳ４１２へ処理を移行する。なお、ステップＳ４０９でｉ≦Ｎの場合はステップＳ４０５へ戻り、ステップＳ４１１でｔ≦Ｔの場合はステップＳ４０４へ戻る。
【００７２】
次に、ステップＳ４１２では、ｔ＝１〜Ｔ、ｋ＝１〜Ｋに対して、式（24）を用いてπｔ，ｋをサンプリングし、記憶部３０のπｔ，ｋを更新し保存する。
【数１２】
2011048583.tif 000014
ここで、κは式（19）で使用されるκと同じであり、また、βｕ＝１−ΣＫｋ＝１βｋである。
【００７３】
（クラスタ間関連度（η）推定部）
図２に示すクラスタ間関連度（η）推定部２５は、ηｋ，ｌのサンプリングを行い、そのサンプリングした値によって記憶部３０に保存してあるηｋ，ｌを更新し記憶する。クラスタ間関連度（η）推定部２５の処理フローについて図５を用いて説明する。まず、ステップＳ５０１では、記憶部３０に現在保存してあるｚｔ，ｉ，ηｋ，ｌ，Ｋとハイパーパラメータξ，Ψと関係データ時系列Ｘとを取得する。
【００７４】
ここで、ｚｔ，ｉ＝ｋ，ｚｔ，ｊ＝ｌとなる（ｔ，ｉ，ｊ）の組の数をＮｋ，ｌ、そのうちｘｔ，ｉ，ｊ＝１となった観測値の数をｎｋ，ｌとし、それらのＮｋ，ｌおよびｎｋ，ｌをカウント変数とする。ステップＳ５０２では、各Ｎｋ，ｌおよび各ｎｋ，ｌの初期値は「０」に初期化（設定）する。
【００７５】
次に、ステップＳ５０３〜Ｓ５１４によって、ｔ＝１〜Ｔ、ｉ＝１〜ＮについてＮｋ，ｌおよびｎｋ，ｌを算出する。具体的には、ステップＳ５０３ではｔ＝１と設定し、ステップＳ５０４ではｉ＝１と設定し、ステップＳ５０５ではｋ＝ｚｔ，ｉと設定し、ステップＳ５０６ではｊ＝１と設定し、ステップＳ５０７ではｌ＝ｚｔ，ｉと設定し、ステップＳ５０８ではＮｋ，ｌを「１」増加するとともにｎｋ，ｌにｘｔ，ｉ，ｊを加算し、ステップＳ５０９ではｊを「１」増加し、ステップＳ５１０ではｊ＞Ｎの場合はステップＳ５１１でｉを「１」増加し、ステップＳ５１２ではｉ＞Ｎの場合はステップＳ５１３でｔを「１」増加し、ステップＳ５１４ではｔ＞Ｔの場合にステップＳ５１５へ処理を移行する。なお、ステップＳ５１０でｊ≦Ｎの場合はステップＳ５０７へ戻り、ステップＳ５１２でｉ≦Ｎの場合はステップＳ５０５へ戻り、ステップＳ５１４でｔ≦Ｔの場合はステップＳ５０４へ戻る。
【００７６】
次に、ステップＳ５１５では、ｋ＝１〜Ｋ、ｌ＝１〜Ｋに対して、式（25）を用いてηｋ，ｌをサンプリングし、記憶部３０のηｋ，ｌを更新し保存する。
【数１３】
2011048583.tif 000015
【００７７】
（混合比（β）推定部）
図２に示す混合比（β）推定部２６は、βのサンプリングを行い、そのサンプリングした値によって記憶部３０に保存してあるβを更新し記憶する。混合比（β）推定部２６の処理フローについて図６を用いて説明する。まず、ステップＳ６０１では、記憶部３０に現在保存してあるｚｔ，ｉ，β，Ｋとハイパーパラメータγ，α０，κとカウント変数ｍｔ，ｋ，ｌとを取得する。
【００７８】
βのサンプリングにおいては、３つの補助変数が必要となる。まず、ステップＳ６０２では、すべてのｔ＝１，〜Ｔ、ｋ＝１〜Ｋ、ｌ＝１〜Ｋに対して、式（５）を用いて補助変数Ｒｔ，ｋ，ｌ＝ｒ，  ｒ∈｛１，２，・・・，ｍｔ，ｋ，ｌ｝をサンプリングする。
【数１４】
2011048583.tif 000016
【００７９】
ここで、ｓ（ｘ，ａ）は第１種スターリング数（unsigned stirling number of the first kind）とよばれる関数であり、ｎ≧ａ≧０に対してｘ（ｘ＋１）（ｘ＋２）・・・（ｘ＋ｎ−１）のｘａの係数をｓ（ｘ，ａ）の値とする関数である。以下の漸化式で与えられる。
ｓ（ｎ，ｋ）＝ｓ（ｎ−１，ｋ−１）＋（ｎ−１）ｓ（ｎ−１，ｋ）
【００８０】
次に、ステップＳ６０３では、ｔ＝１〜Ｔ、ｋ＝１〜Ｋに対し、式（６）を用いて補助変数Ｏｔ，ｋをサンプリングする。
【数１５】
2011048583.tif 000017
ここで、Binomial（）は、二項分布を表す。
【００８１】
次に、ステップＳ６０４では、ｔ＝１〜Ｔ、ｋ＝１〜Ｋ、ｌ＝１〜Ｋに対して、式（７）を用いて補助変数＾Ｒｔ，ｋ，ｌをサンプリングする。
【数１６】
2011048583.tif 000018
【００８２】
最後に、ステップＳ６０５では、式（26）を用いてβをサンプリングし、記憶部３０のβを更新し保存する。
【数１７】
2011048583.tif 000019
【００８３】
（ハイパーパラメータ算出部）
図２に示すハイパーパラメータ算出部２７は、ハイパーパラメータを事前に与えない場合に用いられる。ハイパーパラメータ算出部２７におけるハイパーパラメータの推定演算は、「J.Van Gael, Y.Saatci, Y.W.Teh and Z.Ghahramani,“Beam Sampling for the Infinite Hidden Markov Model”, Proceedings of the 25th International Conference on Machine Learning (ICML), 2008」に記載の技術を用いて実行することが可能である。
【００８４】
（終了判定部）
図２に示す終了判定部２８は、所定の終了条件を満足しているか否かを判定し、その終了条件を満足している場合、演算結果を出力する。終了条件は、例えば、予め決めておいた所定の繰り返し回数になったことである。また、他の終了条件の例として、１つ前のサンプリング処理におけるパラメータの値との差分の絶対値が、予め設定しておいた所定の閾値以下になったこと、としても構わない。また、終了判定部２８は、演算結果の出力においては、所属クラスタｚｔ，ｉ，クラスタ時間変化πｔ，ｋ，クラスタ間関連度ηｋ，ｌ，混合比βの中から、入力装置１０を介してユーザによって出力するように指示された、必要な演算結果のみを表示装置４０に出力する。
【００８５】
（ｄＩＲＭの性能評価例）
本実施形態におけるｄＩＲＭの性能について、シミュレーションによって確認した結果を以下に示す。
【００８６】
クラスタリングの正解の分かっている人工データを作成して、その人工データを用いてｄＩＲＭの定量的な評価を行った。人工データは、全体の時間ステップはＴ＝５、オブジェクト数はＮ＝１６、クラスタ数はＫ＝４とした。オブジェクトのインデックスｉ＝１〜４はほぼ常にクラスタ１に所属し、オブジェクトのインデックスｉ＝５〜８はほぼ常にクラスタ２、オブジェクトのインデックスｉ＝９〜１２はほぼ常にクラスタ３に所属し、オブジェクトのインデックスｉ＝１３〜１６はほぼ常にクラスタ４に所属するように設定する。ただし、一部のオブジェクトは時間に応じてクラスタ間を遷移させた。クラスタ間関連度ηｋ，ｌは、ｐｏｓｉｔｉｖｅなクラスタ間ではη＝０．９、ｎｅｇａｔｉｖｅなクラスタ間ではη＝０．１の２種類を用いた。クラスタ間関連度ηがｐｏｓｉｔｉｖｅかｎｅｇａｔｉｖｅかも事前に設定し、与えられたηに従って各時刻の関係データｘｔ，ｉ，ｊを生成した。
【００８７】
シミュレーションでは、ｄＩＲＭおよび公知のＩＲＭに対して、前記の手続きに従って生成した関係データ時系列Ｘを用いて、オブジェクトのクラスタリングを実行した。ＩＲＭでは、式（13）においてσ＝０．５としてクラスタリングしたのち、ＩＲＭでのオブジェクトのインデックスｉのクラスタリング結果ｚｉを、各時刻ｔでのオブジェクトのインデックスｉのクラスタリングｚｔ，ｉとみなした。
【００８８】
シミュレーションでは、クラスタリングの評価尺度の一つである、ｒａｎｄ  ｉｎｄｅｘを利用した定量的評価を行った。ｒａｎｄ  ｉｎｄｅｘとは、あるデータに対して２つのクラスタリング結果が与えられた時、２つのクラスタリング結果の類似度を測る指標である。ｒａｎｄ  ｉｎｄｅｘの最大値は１で、このとき２つのクラスタリング結果は完全に一致していることを表す。シミュレーションでは、観測データ生成時に利用した正しいクラスタリング結果と、各モデルでの推定結果から得られたクラスタリング結果とのｒａｎｄ  ｉｎｄｅｘを計算した。
【００８９】
ｒａｎｄ  ｉｎｄｅｘの計算結果を図７（ａ）および図７（ｂ）に示す。図７（ａ）は、本実施形態のｄＩＲＭに対するｒａｎｄ  ｉｎｄｅｘの結果を示し、（ｂ）は、公知のＩＲＭに対するｒａｎｄ  ｉｎｄｅｘの結果を示している。図７（ａ）と図７（ｂ）とを比較すると、ｄＩＲＭの場合は、繰り返し回数を増加するにしたがってｒａｎｄ  ｉｎｄｅｘがほぼ１となった。それに対して、ＩＲＭの場合は、繰り返し回数を増加しても、ｒａｎｄ  ｉｎｄｅｘが１になることは無かった。結論として、ｄＩＲＭは、時間変化する関係データに対して、ＩＲＭに比較して、より良くモデル化できることが確認された。
【００９０】
以上、本実施形態のクラスタリング装置８０によれば、公知のＩＲＭを拡張して、クラスタの時間変化を表すパラメータπを導入し、所属クラスタｚの推定方法も時間に従って考慮可能なように改造したことによって、クラスタの時間変化を推定することが可能となった。
【００９１】
なお、本実施形態は、これに限定されるものではなく、その趣旨を変えない範囲で実施することができる。例えば、図１のステップＳ１０４〜Ｓ１０７の演算の順番は、任意の順番で構わない。また、図２の所属クラスタ（ｚ）推定部２３、時間変化情報（π）推定部２４、クラスタ間関連度（η）推定部２５、混合比（β）推定部２６の処理の順番は、任意の順番で構わない。
【００９２】
また、本実施形態において、クラスタリング装置８０（図２参照）の各部２３〜２６の処理は、クラスタリング装置８０をコンピュータで実現したときに搭載されるプログラムによって実現されてもよい。このプログラムは、通信回線を介して提供することもできるし、ＣＤ−ＲＯＭ等のコンピュータ読み取り可能な記録媒体に書き込んで配布することも可能である。
【符号の説明】
【００９３】
１０    入力装置
２０    処理部
２１    初期設定部
２２    推定演算部
２３    所属クラスタ（ｚ）推定部
２４    時間変化情報（π）推定部
２５    クラスタ間関連度（η）推定部
２６    混合比（β）推定部
２８    終了判定部
３０    記憶部
４０    表示装置
８０    クラスタリング装置
９０    クラスタリングシステム
【図１】
2011048583.tif 000024
【図２】
2011048583.tif 000025
【図３】
2011048583.tif 000026
【図４】
2011048583.tif 000027
【図５】
2011048583.tif 000028
【図６】
2011048583.tif 000029
【図７】
2011048583.tif 000030
【図８】
2011048583.tif 000031