"各時刻に観測された複数の特徴量からなる観測値の時系列データに基づいて、時刻のクラスタリングを実行するクラスタリング装置であって、各時刻tにおける各特徴量iからなる観測値xit、潜在クラスタ数K、時刻tが所属するクラスタを示すクラスタ割り当て変数zt、各クラスタkに対して特徴量iが重要であるか否かを示す特徴選択変数rki、あるクラスタに対して特徴量iが重要である確率を示す特徴選択パラメータλi、特徴量iが時刻tの所属するクラスタkに対して重要である場合に観測値xitを生成するために用いられる特徴観測パラメータθki、特徴量iが時刻tの所属するクラスタに対して重要でない場合に観測値xitを生成するために用いられるノイズ観測パラメータφi、時刻t-1においてクラスタkに所属していたオブジェクトが次の時刻tにクラスタlに所属する確率を示すクラスタ遷移パラメータπkl、各クラスタkに所属する時刻の割合を示す混合比パラメータβk、補助変数ut、ptk及びハイパーパラメータγ、α、ai、bi、ci、di、ei、fiを記憶する記憶部と、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴選択パラメータλi、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπki、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ci、diを取得し、各時刻の前記補助変数utを決定すると共に前記潜在クラスタ数Kを決定し、各時刻t及び各クラスタkの前記補助変数ptkを算出して、各時刻tの前記クラスタ割り当て変数ztをサンプリングにより算出し、前記決定した前記補助変数ut及び前記潜在クラスタ数K、前記算出した前記補助変数ptk、並びに前記算出した前記クラスタ割り当て変数ztによって、前記記憶部に記憶してある前記補助変数ut、前記潜在クラスタ数K、前記補助変数ptk、及び前記クラスタ割り当て変数ztを更新し記憶するクラスタ割り当て変数更新手段と、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタk及び各特徴量iの前記特徴選択変数rkiをサンプリングにより算出し、前記算出した前記特徴選択変数rkiによって前記記憶部に記憶してある前記特徴選択変数rkiを更新し記憶する特徴選択変数更新手段と、前記記憶部から前記潜在クラスタ数K、前記特徴選択変数rki、及び前記ハイパーパラメータei、fiを取得し、各特徴量iの前記特徴選択パラメータλiをサンプリングにより算出し、前記算出した前記特徴選択パラメータλiによって前記記憶部に記憶してある前記特徴選択パラメータλiを更新し記憶する特徴選択パラメータ更新手段と、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、及び前記ハイパーパラメータci、diを取得し、各クラスタk及び各特徴量iの前記特徴観測パラメータθkiをサンプリングにより算出し、前記算出した前記特徴観測パラメータθkiによって前記記憶部に記憶してある前記特徴観測パラメータθkiを更新し記憶する特徴観測パラメータ更新手段と、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、及び前記ハイパーパラメータai、biを取得し、各特徴量iの前記ノイズ観測パラメータφiをサンプリングにより算出し、前記算出した前記ノイズ観測パラメータφiによって前記記憶部に記憶してある前記ノイズ観測パラメータφiを更新し記憶するノイズ観測パラメータ更新手段と、前記記憶部から前記潜在クラスタ数K、前記混合比パラメータβk、及び前記ハイパーパラメータαを取得し、各クラスタk、lの前記クラスタ遷移パラメータπklをサンプリングにより算出し、前記算出した前記クラスタ遷移パラメータπklによって前記記憶部に記憶してある前記クラスタ遷移パラメータπklを更新し記憶するクラスタ遷移パラメータ更新手段と、前記記憶部から前記潜在クラスタ数K、及び前記ハイパーパラメータγを取得し、各クラスタkの前記混合比パラメータβkをサンプリングにより算出し、前記算出した前記混合比パラメータβkによって前記記憶部に記憶してある前記混合比パラメータβkを更新し記憶する混合比パラメータ更新手段と、前記クラスタ割り当て変数更新手段、前記特徴選択変数更新手段、前記特徴選択パラメータ更新手段、前記特徴観測パラメータ更新手段、前記ノイズ観測パラメータ更新手段、前記クラスタ遷移パラメータ更新手段、及び前記混合比パラメータ更新手段による更新を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す終了条件判定手段と、を含むクラスタリング装置。前記クラスタ割り当て変数更新手段は、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴選択パラメータλi、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπkl、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ci、diを取得し、各時刻tの前記補助変数utを決定すると共に、前記潜在クラスタ数Kを決定し、各時刻tについて、各特徴量iが時刻tの所属するクラスタkに対して重要であるか否かに基づいて前記補助変数ptkを算出し、前記補助変数ut、ptkを用いてサンプリングして前記クラスタ割り当て変数ztの値を算出し、前記特徴選択変数更新手段は、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタk及び各特徴量iについて、特徴量iがクラスタkに対して重要でないと仮定した場合の適合度及び特徴量iがクラスタkに対して重要であると仮定した場合の適合度を用いてサンプリングして前記特徴選択変数rkiを算出する請求項1記載のクラスタリング装置。前記クラスタ割り当て変数更新手段は、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴選択パラメータλi、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπkl、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ci、diを取得し、各時刻の前記補助変数utを決定すると共に、前記潜在クラスタ数Kを決定し、時刻t=1,2,...,Tの順番で、全てのk=1,2,...,Kに対して以下のI式に示す漸化式に従って、前記補助変数ptkを算出し、時刻t=T,T-1,...,1の順番で、以下のII式に示す各クラスタkに対する確率pzt=k|zt1,X、R、Φ、Θ、Π、Uからサンプリングして前記クラスタ割り当て変数ztの値を算出し、前記特徴選択変数更新手段は、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタk及び各特徴量iについて、以下のIII式に示す確率prki=0|X、Z、R\、Λ、Φ、Θ及び以下のIV式に示す確率prki=1|X、Z、R\、Λ、Φ、Θからサンプリングして前記特徴選択変数rkiを算出する請求項1記載のクラスタリング装置。ここでBernoulliはベルヌーイ分布を表す。Iは、条件cが成立した場合に1そうでない場合には0を返す関数である。前記記憶部から前記潜在クラスタ数K、前記特徴選択変数rki、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、前記混合比パラメータβk、及び前記ハイパーパラメータγ、α、ai、bi、ci、di、ei、fiを取得し、ハイパーパラメータγをサンプリングにより算出し、ハイパーパラメータαをサンプリングにより算出し、各特徴量iのハイパーパラメータaiをサンプリングにより算出し、各特徴量iのハイパーパラメータbiをサンプリングにより算出し、各特徴量iのハイパーパラメータciをサンプリングにより算出し、各特徴量iのハイパーパラメータdiをサンプリングにより算出し、各特徴量iのハイパーパラメータeiをサンプリングにより算出し、各特徴量iのハイパーパラメータfiをサンプリングにより算出し、前記算出した前記ハイパーパラメータγ、α、ai、bi、ci、di、ei、fiによって前記記憶部に記憶してある前記ハイパーパラメータγ、α、ai、bi、ci、di、ei、fiを更新し記憶するハイパーパラメータ更新手段を更に含み、前記終了条件判定手段は、前記クラスタ割り当て変数更新手段、前記特徴選択変数更新手段、前記特徴選択パラメータ更新手段、前記特徴観測パラメータ更新手段、前記ノイズ観測パラメータ更新手段、前記クラスタ遷移パラメータ更新手段、前記混合比パラメータ更新手段、及び前記ハイパーパラメータ更新手段による更新を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す請求項1~請求項3の何れか1項記載のクラスタリング装置。各時刻に観測された複数の特徴量からなる観測値の時系列データに基づいて、時刻のクラスタリングを実行するクラスタリング装置におけるクラスタリング方法であって、前記クラスタリング装置は、各時刻tにおける各特徴量iからなる観測値xit、潜在クラスタ数K、時刻tが所属するクラスタを示すクラスタ割り当て変数zt、各クラスタkに対して特徴量iが重要であるか否かを示す特徴選択変数rki、あるクラスタに対して特徴量iが重要である確率を示す特徴選択パラメータλi、特徴量iが時刻tの所属するクラスタkに対して重要である場合に観測値xitを生成するために用いられる特徴観測パラメータθki、特徴量iが時刻tの所属するクラスタに対して重要でない場合に観測値xitを生成するために用いられるノイズ観測パラメータφi、時刻t-1においてクラスタkに所属していたオブジェクトが次の時刻tにクラスタlに所属する確率を示すクラスタ遷移パラメータπkl、各クラスタkに所属する時刻の割合を示す混合比パラメータβk、補助変数ut、ptk及びハイパーパラメータγ、α、ai、bi、ci、di、ei、fiを記憶する記憶部と解析部とを備え、前記解析部は、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴選択パラメータλi、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπki、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ci、diを取得し、各時刻の前記補助変数utを決定すると共に前記潜在クラスタ数Kを決定し、各時刻t及び各クラスタkの前記補助変数ptkを算出して、各時刻tの前記クラスタ割り当て変数ztをサンプリングにより算出し、前記決定した前記補助変数ut及び前記潜在クラスタ数K、前記算出した前記補助変数ptk、並びに前記算出した前記クラスタ割り当て変数ztによって、前記記憶部に記憶してある前記補助変数ut、前記潜在クラスタ数K、前記補助変数ptk、及び前記クラスタ割り当て変数ztを更新し記憶するクラスタ割り当て変数更新ステップと、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタk及び各特徴量iの前記特徴選択変数rkiをサンプリングにより算出し、前記算出した前記特徴選択変数rkiによって前記記憶部に記憶してある前記特徴選択変数rkiを更新し記憶する特徴選択変数更新ステップと、前記記憶部から前記潜在クラスタ数K、前記特徴選択変数rki、及び前記ハイパーパラメータei、fiを取得し、各特徴量iの前記特徴選択パラメータλiをサンプリングにより算出し、前記算出した前記特徴選択パラメータλiによって前記記憶部に記憶してある前記特徴選択パラメータλiを更新し記憶する特徴選択パラメータ更新ステップと、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、及び前記ハイパーパラメータci、diを取得し、各クラスタk及び各特徴量iの前記特徴観測パラメータθkiをサンプリングにより算出し、前記算出した前記特徴観測パラメータθkiによって前記記憶部に記憶してある前記特徴観測パラメータθkiを更新し記憶する特徴観測パラメータ更新ステップと、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、及び前記ハイパーパラメータai、biを取得し、各特徴量iの前記ノイズ観測パラメータφiをサンプリングにより算出し、前記算出した前記ノイズ観測パラメータφiによって前記記憶部に記憶してある前記ノイズ観測パラメータφiを更新し記憶するノイズ観測パラメータ更新ステップと、前記記憶部から前記潜在クラスタ数K、前記混合比パラメータβk、及び前記ハイパーパラメータαを取得し、各クラスタk、lの前記クラスタ遷移パラメータπklをサンプリングにより算出し、前記算出した前記クラスタ遷移パラメータπklによって前記記憶部に記憶してある前記クラスタ遷移パラメータπklを更新し記憶するクラスタ遷移パラメータ更新ステップと、前記記憶部から前記潜在クラスタ数K、及び前記ハイパーパラメータγを取得し、各クラスタkの前記混合比パラメータβkをサンプリングにより算出し、前記算出した前記混合比パラメータβkによって前記記憶部に記憶してある前記混合比パラメータβkを更新し記憶する混合比パラメータ更新ステップと、前記クラスタ割り当て変数更新ステップ、前記特徴選択変数更新ステップ、前記特徴選択パラメータ更新ステップ、前記特徴観測パラメータ更新ステップ、前記ノイズ観測パラメータ更新ステップ、前記クラスタ遷移パラメータ更新ステップ、及び前記混合比パラメータ更新ステップによる更新を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す終了条件判定ステップと、を実行することを特徴とするクラスタリング方法。前記クラスタ割り当て変数更新ステップは、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴選択パラメータλi、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπkl、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ci、diを取得し、各時刻tの前記補助変数utを決定すると共に、前記潜在クラスタ数Kを決定し、各時刻tについて、各特徴量iが時刻tの所属するクラスタkに対して重要であるか否かに基づいて前記補助変数ptkを算出し、前記補助変数ut、ptkを用いてサンプリングして前記クラスタ割り当て変数ztの値を算出し、前記特徴選択変数更新ステップは、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタk及び各特徴量iについて、特徴量iがクラスタkに対して重要でないと仮定した場合の適合度及び特徴量iがクラスタkに対して重要であると仮定した場合の適合度を用いてサンプリングして前記特徴選択変数rkiを算出する請求項5記載のクラスタリング方法。前記記憶部から前記潜在クラスタ数K、前記特徴選択変数rki、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、前記混合比パラメータβk、及び前記ハイパーパラメータγ、α、ai、bi、ci、di、ei、fiを取得し、ハイパーパラメータγをサンプリングにより算出し、ハイパーパラメータαをサンプリングにより算出し、各特徴量iのハイパーパラメータaiをサンプリングにより算出し、各特徴量iのハイパーパラメータbiをサンプリングにより算出し、各特徴量iのハイパーパラメータciをサンプリングにより算出し、各特徴量iのハイパーパラメータdiをサンプリングにより算出し、各特徴量iのハイパーパラメータeiをサンプリングにより算出し、各特徴量iのハイパーパラメータfiをサンプリングにより算出し、前記算出した前記ハイパーパラメータγ、α、ai、bi、ci、di、ei、fiによって前記記憶部に記憶してある前記ハイパーパラメータγ、α、ai、bi、ci、di、ei、fiを更新し記憶するハイパーパラメータ更新ステップを更に含み、前記終了条件判定ステップは、前記クラスタ割り当て変数更新ステップ、前記特徴選択変数更新ステップ、前記特徴選択パラメータ更新ステップ、前記特徴観測パラメータ更新ステップ、前記ノイズ観測パラメータ更新ステップ、前記クラスタ遷移パラメータ更新ステップ、前記混合比パラメータ更新ステップ、及び前記ハイパーパラメータ更新ステップによる更新を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す請求項5又は6記載のクラスタリング方法。請求項1~請求項4の何れか1項に記載のクラスタリング装置を構成する各手段として、コンピュータを機能させることを特徴とするプログラム。特に近年webなどで入手されるデジタルデータは非常に高次元かつスパースである。そのような状態では、各時刻で観測される特徴ベクトルのほとんどの要素はまったく同じ値をもち、ごく一部だけが変化するという状況に起こる。したがって、本当にクラスタに区別に有用な観測量のみを抜き出してクラスタリングをしなければ、多数の無益な観測量にクラスタ構造が隠されてしまう、という特徴選択の問題がある。上記の非特許文献3、4に記載の手法を用いても、上記の特徴選択の問題は解決されない。次に特徴選択手法の問題であるが、特に時系列データに適用する際にはこれらの手法の仮定するデータの独立性が問題となる。これは、クラスタリングされる各観測量ベクトルは互いに独立かつ一様に分布する、という仮定である。多くの課題においてこの仮定は問題を起こさないが、時系列データはその例外の一つである。時系列データは本質的に過去の情報の影響を受けるため、常に前の時刻との依存関係が発生する。したがって、この依存関係をモデル化しない手法は時系列データの特性を無視していることになる。このように、既存の時系列データのクラスタリングモデルと特徴選択手法はいずれも一長一短であることが分かる。本発明は、上記の事情を鑑みてなされたもので、時系列データが与えられたときに、クラスタ数を自動的に決定すると共に、時刻のクラスタリングに寄与する真の特徴量と、時刻のクラスリングに寄与しないノイズ特徴の影響を自動的に推定分離して、時刻のクラスタリングを行うことができるクラスタリング装置、方法、及びプログラムを提供することを目的とする。上記の目的を達成するためにクラスタリング装置は、各時刻に観測された複数の特徴量からなる観測値の時系列データに基づいて、時刻のクラスタリングを実行するクラスタリング装置であって、各時刻tにおける各特徴量iからなる観測値xit、潜在クラスタ数K、時刻tが所属するクラスタを示すクラスタ割り当て変数zt、各クラスタkに対して特徴量iが重要であるか否かを示す特徴選択変数rki、あるクラスタに対して特徴量iが重要である確率を示す特徴選択パラメータλi、特徴量iが時刻tの所属するクラスタkに対して重要である場合に観測値xitを生成するために用いられる特徴観測パラメータθki、特徴量iが時刻tの所属するクラスタに対して重要でない場合に観測値xitを生成するために用いられるノイズ観測パラメータφi、時刻t-1においてクラスタkに所属していたオブジェクトが次の時刻tにクラスタlに所属する確率を示すクラスタ遷移パラメータπkl、各クラスタkに所属する時刻の割合を示す混合比パラメータβk、補助変数ut、ptk及びハイパーパラメータγ、α、ai、bi、ci、di、ei、fiを記憶する記憶部と、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴選択パラメータλi、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπki、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ci、diを取得し、各時刻の前記補助変数utを決定すると共に前記潜在クラスタ数Kを決定し、各時刻t及び各クラスタkの前記補助変数ptkを算出して、各時刻tの前記クラスタ割り当て変数ztをサンプリングにより算出し、前記決定した前記補助変数ut及び前記潜在クラスタ数K、前記算出した前記補助変数ptk、並びに前記算出した前記クラスタ割り当て変数ztによって、前記記憶部に記憶してある前記補助変数ut、前記潜在クラスタ数K、前記補助変数ptk、及び前記クラスタ割り当て変数ztを更新し記憶するクラスタ割り当て変数更新手段と、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタk及び各特徴量iの前記特徴選択変数rkiをサンプリングにより算出し、前記算出した前記特徴選択変数rkiによって前記記憶部に記憶してある前記特徴選択変数rkiを更新し記憶する特徴選択変数更新手段と、前記記憶部から前記潜在クラスタ数K、前記特徴選択変数rki、及び前記ハイパーパラメータei、fiを取得し、各特徴量iの前記特徴選択パラメータλiをサンプリングにより算出し、前記算出した前記特徴選択パラメータλiによって前記記憶部に記憶してある前記特徴選択パラメータλiを更新し記憶する特徴選択パラメータ更新手段と、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、及び前記ハイパーパラメータci、diを取得し、各クラスタk及び各特徴量iの前記特徴観測パラメータθkiをサンプリングにより算出し、前記算出した前記特徴観測パラメータθkiによって前記記憶部に記憶してある前記特徴観測パラメータθkiを更新し記憶する特徴観測パラメータ更新手段と、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、及び前記ハイパーパラメータai、biを取得し、各特徴量iの前記ノイズ観測パラメータφiをサンプリングにより算出し、前記算出した前記ノイズ観測パラメータφiによって前記記憶部に記憶してある前記ノイズ観測パラメータφiを更新し記憶するノイズ観測パラメータ更新手段と、前記記憶部から前記潜在クラスタ数K、前記混合比パラメータβk、及び前記ハイパーパラメータαを取得し、各クラスタk、lの前記クラスタ遷移パラメータπklをサンプリングにより算出し、前記算出した前記クラスタ遷移パラメータπklによって前記記憶部に記憶してある前記クラスタ遷移パラメータπklを更新し記憶するクラスタ遷移パラメータ更新手段と、前記記憶部から前記潜在クラスタ数K、及び前記ハイパーパラメータγを取得し、各クラスタkの前記混合比パラメータβkをサンプリングにより算出し、前記算出した前記混合比パラメータβkによって前記記憶部に記憶してある前記混合比パラメータβkを更新し記憶する混合比パラメータ更新手段と、前記クラスタ割り当て変数更新手段、前記特徴選択変数更新手段、前記特徴選択パラメータ更新手段、前記特徴観測パラメータ更新手段、前記ノイズ観測パラメータ更新手段、前記クラスタ遷移パラメータ更新手段、及び前記混合比パラメータ更新手段による更新を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す終了条件判定手段と、を含んで構成されている。本発明に係るクラスタリング方法は、各時刻に観測された複数の特徴量からなる観測値の時系列データに基づいて、時刻のクラスタリングを実行するクラスタリング装置におけるクラスタリング方法であって、前記クラスタリング装置は、各時刻tにおける各特徴量iからなる観測値xit、潜在クラスタ数K、時刻tが所属するクラスタを示すクラスタ割り当て変数zt、各クラスタkに対して特徴量iが重要であるか否かを示す特徴選択変数rki、あるクラスタに対して特徴量iが重要である確率を示す特徴選択パラメータλi、特徴量iが時刻tの所属するクラスタkに対して重要である場合に観測値xitを生成するために用いられる特徴観測パラメータθki、特徴量iが時刻tの所属するクラスタに対して重要でない場合に観測値xitを生成するために用いられるノイズ観測パラメータφi、時刻t-1においてクラスタkに所属していたオブジェクトが次の時刻tにクラスタlに所属する確率を示すクラスタ遷移パラメータπkl、各クラスタkに所属する時刻の割合を示す混合比パラメータβk、補助変数ut、ptk及びハイパーパラメータγ、α、ai、bi、ci、di、ei、fiを記憶する記憶部と解析部とを備え、前記解析部は、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴選択パラメータλi、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπki、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ci、diを取得し、各時刻の前記補助変数utを決定すると共に前記潜在クラスタ数Kを決定し、各時刻t及び各クラスタkの前記補助変数ptkを算出して、各時刻tの前記クラスタ割り当て変数ztをサンプリングにより算出し、前記決定した前記補助変数ut及び前記潜在クラスタ数K、前記算出した前記補助変数ptk、並びに前記算出した前記クラスタ割り当て変数ztによって、前記記憶部に記憶してある前記補助変数ut、前記潜在クラスタ数K、前記補助変数ptk、及び前記クラスタ割り当て変数ztを更新し記憶するクラスタ割り当て変数更新ステップと、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、前記特徴観測パラメータθki、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタk及び各特徴量iの前記特徴選択変数rkiをサンプリングにより算出し、前記算出した前記特徴選択変数rkiによって前記記憶部に記憶してある前記特徴選択変数rkiを更新し記憶する特徴選択変数更新ステップと、前記記憶部から前記潜在クラスタ数K、前記特徴選択変数rki、及び前記ハイパーパラメータei、fiを取得し、各特徴量iの前記特徴選択パラメータλiをサンプリングにより算出し、前記算出した前記特徴選択パラメータλiによって前記記憶部に記憶してある前記特徴選択パラメータλiを更新し記憶する特徴選択パラメータ更新ステップと、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、及び前記ハイパーパラメータci、diを取得し、各クラスタk及び各特徴量iの前記特徴観測パラメータθkiをサンプリングにより算出し、前記算出した前記特徴観測パラメータθkiによって前記記憶部に記憶してある前記特徴観測パラメータθkiを更新し記憶する特徴観測パラメータ更新ステップと、前記記憶部から前記観測値xit、前記潜在クラスタ数K、前記クラスタ割り当て変数zt、前記特徴選択変数rki、及び前記ハイパーパラメータai、biを取得し、各特徴量iの前記ノイズ観測パラメータφiをサンプリングにより算出し、前記算出した前記ノイズ観測パラメータφiによって前記記憶部に記憶してある前記ノイズ観測パラメータφiを更新し記憶するノイズ観測パラメータ更新ステップと、前記記憶部から前記潜在クラスタ数K、前記混合比パラメータβk、及び前記ハイパーパラメータαを取得し、各クラスタk、lの前記クラスタ遷移パラメータπklをサンプリングにより算出し、前記算出した前記クラスタ遷移パラメータπklによって前記記憶部に記憶してある前記クラスタ遷移パラメータπklを更新し記憶するクラスタ遷移パラメータ更新ステップと、前記記憶部から前記潜在クラスタ数K、及び前記ハイパーパラメータγを取得し、各クラスタkの前記混合比パラメータβkをサンプリングにより算出し、前記算出した前記混合比パラメータβkによって前記記憶部に記憶してある前記混合比パラメータβkを更新し記憶する混合比パラメータ更新ステップと、前記クラスタ割り当て変数更新ステップ、前記特徴選択変数更新ステップ、前記特徴選択パラメータ更新ステップ、前記特徴観測パラメータ更新ステップ、前記ノイズ観測パラメータ更新ステップ、前記クラスタ遷移パラメータ更新ステップ、及び前記混合比パラメータ更新ステップによる更新を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す終了条件判定ステップと、を実行することを特徴とする。本発明によれば、クラスタ割り当て変数を更新するときにクラスタ数を決定すると共に、クラスタ割り当て変数の更新、特徴選択変数の更新、特徴選択パラメータの更新、特徴観測パラメータの更新、ノイズ観測パラメータの更新、クラスタ遷移パラメータの更新、及び混合比パラメータの更新を繰り返すことにより、時系列データが与えられたときに、クラスタ数を自動的に決定すると共に、時刻のクラスタリングに寄与する真の特徴量と、時刻のクラスリングに寄与しないノイズ特徴の影響を自動的に推定分離して、時刻のクラスタリングを行うことができる。本発明に係るプログラムは、上記のクラスタリング装置の各手段としてコンピュータを機能させるためのプログラムである。以上説明したように、本発明のクラスタリング装置、方法、及びプログラムによれば、クラスタ割り当て変数を更新するときにクラスタ数を決定すると共に、クラスタ割り当て変数の更新、特徴選択変数の更新、特徴選択パラメータの更新、特徴観測パラメータの更新、ノイズ観測パラメータの更新、クラスタ遷移パラメータの更新、及び混合比パラメータの更新を繰り返すことにより、時系列データが与えられたときに、クラスタ数を自動的に決定すると共に、時刻のクラスタリングに寄与する真の特徴量と、時刻のクラスリングに寄与しないノイズ特徴の影響を自動的に推定分離して、時刻のクラスタリングを行うことができる、という効果が得られる。以下、図面を参照して本発明の実施の形態を詳細に説明する。<概要>まず、本発明の概要について説明する。これまで、時系列データに対して各時刻を少数の潜在的なクラスタに帰属することでモデル化する手法が多数提案されている。この枠組みでは、一般に多次元ベクトルとして表現可能な観測量が各時刻で観測される。そして、この特徴量の分布の類似度でもって各時刻を少数のパターンクラスタに分類する。この隠れたクラスタの推定によって、例えば人間の動作や音声といった複雑な時系列データを精度良く解析認識することが可能である。本発明では、これらの手法のもつ欠点のうち、大きく2つの部分を改善した手法を提案する。まず一点目として、観測量のうち、時刻のクラスタリングに寄与しないノイズ特徴量の存在を仮定し、その影響を自動的に排除する。二点目として、クラスタの数を事前に決定することなく、与えられたデータを表現する上で最も適切なクラスタ数を自動的に決定する。図1は提案するクラスタリングシステムのコンセプトを模式的に表した図である。クラスタリングしたい観測値の時系列データをXとする。観測値の時系列データはN種類の特徴量、すなわちN次元の特徴ベクトルが総時間ステップ数Tに渡って観測されたものとする。特徴量をi12N、時間ステップをt1,2Tでインデキシングする。時刻tにおけるi番目の特徴量の観測値をxtiとする。すなわち、X=xi,ti=1Nt=1Nである。xitは01のバイナリデータ、あるいは自然数、実数、あるいはシンボルデータである。この観測データと、ハイパーパラメータと呼ばれる定数とを入力することで、クラスタリングシステムから解析結果が出力される。主要な解析出力はクラスタ数K、クラスタ割り当て変数Z、特徴選択変数Rである。クラスタ数KはXに内在する特徴観測パターン、すなわち潜在クラスタ数である。Zzt,t1Tは、全時刻ステップに対する潜在クラスタの割り当てを表す。ztk,k1Kは、時刻ステップkの特徴量がクラスタkから確率的に生成されたことを表す。また、文脈に応じてztは1ofKベクトルとしても利用する。すなわち、ztkのとき、ztk1,ztl≠k0とする。特徴選択変数Rrki,k1K,i1Nは全クラスタにおける特徴選択結果を表す。rki∈01は、i番目の特徴量がk番目のクラスタの特性を表す真の特徴か、そうでないノイズ特徴かを表す。図1ではrzti1となる部分を濃く表示している。ほとんどの白い部分は重要でないノイズ特徴と判断されている。k6のクラスタは、特徴的なものが存在しない、すなわちノイズのみが観測された期間であると解釈される。<システム構成>図2は、本発明の実施の形態に係るクラスタリング装置100の模式図の一例である。クラスタリング装置100は、解析したい観測データの入力を受け付ける入力部1と、クラスタリングを行う解析部2と、入力および計算された変数を記憶する記憶部3と、クラスタリング結果を出力する出力部4と、を備えている。入力部1は、既知のビデオカメラ、マイクロフォン、記憶装置などの入力器により実現される。また、出力部4はディスプレイ、プリンタ、磁気ディスクなどで実装される。解析部2は、初期設定部21、変数計算部22、及び解析終了条件判定部211を備えている。変数計算部22は、クラスタ割り当て変数更新部23、特徴選択変数更新部24、特徴選択パラメータ更新部25、特徴観測パラメータ更新部26、ノイズ観測パラメータ更新部27、クラスタ遷移パラメータ更新部28、混合比パラメータ更新部29、及びハイパーパラメータ更新部210を備えている。記憶部3は、定数31、観測値32、潜在クラスタ数33、クラスタ割り当て変数34、特徴選択変数35、特徴選択パラメータ36、特徴観測パラメータ37、ノイズ観測パラメータ38、クラスタ遷移パラメータ39、混合比パラメータ310、補助変数311、及びハイパーパラメータ312から構成される。初期設定部21は、入力部1より入力された情報を記憶部3の各部へと記憶するとともに、記憶部3の各部の値の初期化を行う。変数計算部22は、記憶部3に保存された情報を利用し、記憶部3の各部の情報を更新、保存する。クラスタ割り当て変数更新部23は記憶部3の情報を用いて潜在クラスタ数33、クラスタ割り当て変数34、特徴選択変数35、特徴観測パラメータ37、クラスタ遷移パラメータ39、混合比パラメータ310、補助変数311を更新し、保存する。特徴選択変数更新部24は、記憶部3の情報を用いて特徴選択変数35を更新し、保存する。特徴選択パラメータ更新部25は、記憶部3の情報を用いて特徴選択パラメータ36を更新し、保存する。特徴観測パラメータ更新部26は記憶部3の情報を用いて特徴観測パラメータ37を更新し、保存する。ノイズ観測パラメータ更新部27は記憶部3の情報を用いてノイズ観測パラメータ38を更新し、保存する。クラスタ遷移パラメータ更新部28は、記憶部3の情報を用いてクラスタ遷移パラメータ39を更新し、保存する。混合比パラメータ更新部29は、記憶部3の情報を用いて混合比パラメータ310を更新し、保存する。ハイパーパラメータ更新部210は記憶部3の情報を用いてハイパーパラメータ312を更新し、保存する。変数計算部22の処理の流れは、クラスタ割り当て変数更新部23、特徴選択変数更新部24、特徴選択パラメータ更新部25、特徴観測パラメータ更新部26、ノイズ観測パラメータ更新部27、クラスタ遷移パラメータ更新部28、混合比パラメータ更新部29、ハイパーパラメータ更新部210の順番に限らず、任意の順番で実施してよい。また、必要がなければハイパーパラメータ更新部210による更新は実施しなくてもよい。解析終了条件判定部211は、変数計算部22および記憶部3を監視し、解析終了条件を満たしたかどうかを計算する。出力部4は解析終了条件判定部211を監視し、解析終了条件判定部211により解析終了条件が満たされたと判断された場合、クラスタリング結果の出力を実施する。<クラスタリング装置の動作>次に、本実施の形態に係るクラスタリング装置100の作用について説明する。まず、各種の定数及び観測値の時系列データが入力部1を介してクラスタリング装置100に入力されると、クラスタリング装置100において、図3に示すクラスタリング処理ルーチンが実行される。まず、ステップS101において、初期設定部21によって、入力部1により入力された、定数及び観測値の時系列データを受け付けて、記憶部3に格納する。そして、ステップS102において、初期設定部21によって、記憶部3の各種パラメータを初期化する。ステップS103では、繰り返し回数を示す変数itrに初期値1を設定する。次のステップS104では、クラスタ割り当て変数更新部23によって、記憶部3の情報に基づいて、クラスタ割り当て変数を更新して、記憶部3に記憶する。ステップS105では、特徴選択変数更新部24によって、記憶部3の情報に基づいて、特徴選択変数を更新して、記憶部3に記憶する。ステップS106では、特徴選択パラメータ更新部25によって、記憶部3の情報に基づいて、特徴選択パラメータを更新して、記憶部3に記憶する。そして、ステップS107において、特徴観測パラメータ更新部26によって、記憶部3の情報に基づいて、特徴観測パラメータを更新して、記憶部3に記憶する。ステップS108では、ノイズ観測パラメータ更新部27によって、記憶部3の情報に基づいて、ノイズ観測パラメータを更新して、記憶部3に記憶する。ステップS109では、クラスタ遷移パラメータ更新部28によって、記憶部3の情報に基づいて、クラスタ遷移パラメータを更新して、記憶部3に記憶する。次のステップS110では、混合比パラメータ更新部29によって、記憶部3の情報に基づいて、混合比パラメータを更新して、記憶部3に記憶する。ステップS111において、ハイパーパラメータ更新部210によって、記憶部3の情報に基づいて、ハイパーパラメータを更新して、記憶部3に記憶する。そして、ステップS112において、解析終了条件判定部211によって、解析終了条件を満足したか否かを判定する。解析終了条件は、例えば、繰り返し変数itrが、予め決めておいた繰り返し回数になったことである。解析終了条件を満足していない場合には、ステップS113において、変数itrを1インクリメントして、上記ステップS104へ戻る。一方、終了条件を満足する場合には、ステップS114において、上記ステップS104で最終的に更新されたクラスタ割り当て変数に基づくクラスタリング結果を出力部4により出力して、クラスタリング処理ルーチンを終了する。次に、図3の各ステップの処理を実行する、クラスタリング装置100の構成および処理の詳細について説明する。<初期設定部21>初期設定部21は、入力部1を通じて、解析したい観測データの時系列データXおよび複数の定数の入力を受け付ける。入力データは以下のような形式であり、初期設定部21でロードされる。観測値の時系列データは、N種類の特徴量、すなわちN次元の特徴ベクトルが総時間ステップ数Tに渡って観測されたものとする。特徴量をi12N、時間ステップをt12Tでインデキシングする。時刻tにおける観測値のi番目の特徴量をxitとする。すなわち、Xxit、i1N、t1Nである。xitは01のバイナリデータ、あるいは自然数、実数、あるいはシンボルデータである。また、観測値の時系列データとは別に、以下の定数も別途入力が必要である。Tは観測データの最大時間ステップサイズである。Nは観測データの特徴量のサイズ、すなわち特徴ベクトルの次元数である。K0は潜在クラスタ数33の初期値である。Lは解析終了条件判定部211に使用する終了条件である。これは例えば変数計算部22の最大繰り返し回数、あるいは何らかの評価値の変化率の閾値など、ユーザの目的に応じて定める。さらに、モデルに必要なハイパーパラメータと呼ばれる定数も同時に入力する。γは混合比パラメータ310の計算で必要となる。αはクラスタ遷移パラメータ309の計算で必要となる。a1aiaN,b1bibNは主にノイズ観測パラメータ38の計算で必要となる。c1cicN,d1didNは主に特徴観測パラメータ37の計算で必要となる。e1eieN,f1fifNは主に特徴選択パラメータ36の計算で必要となる。各ハイパーパラメータの定義域はxitの形式およびそれに応じたモデル化に従う。例えば最も単純なバイナリデータの例、すなわちxit∈01の場合、γは0より大きい実数、aibicidieifiは1以上の実数が望ましい。これらの入力を受け付けたのち、初期設定部21は、記憶部3の各部の値を初期化する。定数31には入力されたTNK0Lを設定する。観測値32には入力されたXを設定する。潜在クラスタ数33にはK0を設定する。潜在クラスタ割り当て変数34には、1Kの値を各ztにランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。特徴選択変数35には、設定した潜在クラスタ数Kに従って各rkiに対し0または1の値をランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。特徴選択パラメータ36には、各λiに対し0以上1未満の値をランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。特徴観測パラメータ37には、xitが0または1のバイナリデータの場合、設定されたKに従って各θkiに対し0以上1未満の値をランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。xitが別の形式の場合、適切に選択されたモデルに合わせて値を設定する。ノイズ観測パラメータ38には、xitが0または1のバイナリデータの場合、各φiに対し0以上1未満の値をランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。xitが別の形式の場合、適切に選択されたモデルに合わせて値を設定する。クラスタ遷移パラメータ39には、設定されたKに従って各πklに対し0以上1未満の値をランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。ただし、全てのk∈1K1に対し、ΣK1l1πkl1が満たされる必要がある。混合比パラメータ310には、設定されたKに従って各βkに対し0以上1未満の値をランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。ただし、ΣK1k1πk,lk1が満たされる必要がある。補助変数311は、潜在クラスタ数33、クラスタ割り当て変数34、クラスタ遷移パラメータ39の設定後に、後述するモデルに従って設定する。ハイパーパラメータ312には、入力されたハイパーパラメータを設定する。なお、xi,tの形式に従って、変数計算部22の具体的な実装および記憶部3の各部のメモリサイズは変更されうる。本実施の形態では、最も単純なバイナリデータの例、すなわちxi,t∈0,1の場合の例について説明するが、他のデータ形式の場合も一般的な教科書、例えば参考文献ビショップCMパターン認識と機械学習シュプリンガージャパン2007を参照することで簡単に実装することが可能である。初期化処理の終了後、処理を変数計算部22に移す。<モデルの定式化記憶部3>それでは主に記憶部3に保存される情報について、モデルにおけるそれらの意味とともに説明する。なお、ここではxi,tがバイナリデータ、すなわちxi,t∈01の場合を例に説明する。本発明では、観測された時系列データXを少数の観測値のパターン、すなわち潜在クラスタに分解して、データの大まかな生成パターンをとらえやすくする。また、それによってノイズに対してロバストな予測モデルの構築を実現する。まず、潜在クラスタ数33はデータに内在するパターン、つまりクラスタの数を表す。クラスタ数Kは一般には未知であるが、本モデルでは真のクラスタ数がどのような値でも適切に推定解析可能なディリクレ過程混合モデル、より具体的にはそのHMMモデル応用であるiHMMを利用する。まず、全T個の時刻ステップのうち、クラスタkによって説明されている時刻ステップの割合を表す混合比パラメータ310を、以下の1式に示す確率モデルでサンプリングする。ここで、β=β1,...,βK,...Tは無限次元のベクトルである。StickはStickbreakingprocessと呼ばれる確率過程である。より具体的には各βkは以下の2式に示す確率を用いて、以下の3式に従って決定される。ここでBetaはベータ分布を表す。上記3式によれば、加算無限個のクラスタ数を考えたときにベクトルβは全ての要素が0以上かつ、足して1となる。さて、実際の解析時にはKは有限の値をとる。すなわち、データは有限なので必ず有限K個のクラスタで表現できるはずである。実際の実装上ではβはK1次元ベクトルとなる。βk≦Kは、全T時刻のうち、クラスタkで説明される時刻がどれだけあるかを表す割合である。βK1はunrepresentedclusters、すなわち理論上は存在するが現在はデータのモデル化に不要なクラスタ集合のもつ重みである。βK1は次の4式で決定する。Stickbreakingprocessの定義より、どのようなKに対してもベクトルβは全ての要素が0以上かつ、足して1となる。混合比パラメータ310にはこのK1次元のベクトルが格納される。次に、このモデルでは時系列データの時間発展をクラスタ割り当ての時間遷移で表現する。すなわち、複雑な時間変動をおこなう時系列データを、たかだかK種類のパターン間の遷移で表現する。この遷移は、クラスタ遷移パラメータ39によって制御される。クラスタ遷移パラメータπk、lは時刻t-1におけるクラスタがkだった場合、時刻tにおけるクラスタ割り当てがlになる確率をあらわす。クラスタkからの遷移先の確率をπk=πk,1,πk,2,...,πk,l,...,Tで表現する。このベクトルも理論上は無限次元なので、以下の5式に示すように、クラスタ遷移パラメータの定義はディリクレ過程を必要とする。DPは有限次元におけるDirichlet分布を拡張したものでStickbreakingprocessと同様無限次元のベクトルを生成する。実際にはβはK1次元で有限なので、実装上は以下の6式で計算を行う。Dirichletは有限次元のディリクレ分布である。得られるπkはK1次元のベクトルであり、これをクラスタ遷移パラメータ39として格納する。このパラメータを利用して、実際のクラスタ遷移を実現する。この結果得られるクラスタ割り当ての時系列はクラスタ割り当て変数34に記録される。クラスタ割り当て変数zt=kは時刻tにおけるクラスタ割り当てである。ztは以下の7式で決定する。ここでMultinomialは多項分布を表す。なお、z1の決定については種々の方法があるが、ここでは便宜的にz01を仮定して上記7式を適用することとする。これで、与えられた時系列データXを、時間軸方向にK種のクラスタへと分解することができる。次に、各クラスタについてN種の観測量のうち真に役に立つクラスタの特徴を表す観測量とノイズでしかない観測量とを分離するモデルを説明する。まず、i番目の観測値があるクラスタにとって重要な観測量、すなわち特徴である確率をλiとする。この値は以下の8式で決定する。上記8式に従ってサンプリングされるλiの値を特徴選択パラメータ36に格納する。次に、各クラスタkにおいて、i番目の観測値を重要な特徴とみなすか、あるいは単なるノイズ観測量とみなすかを表す特徴選択変数35を、以下の9式で決定する。ここでBernoulliはベルヌーイ分布を表す。rk,i=1のとき、i番目の観測値はクラスタkにとって重要な特徴であり、rk,i=0のときはただのノイズ観測量と考える。このrk,iを全てのkiについて確率的にサンプリングし、特徴選択変数35に格納する。これによって、時刻tの観測値のi番目の特徴量xi,tを、その時刻tが所属するクラスタ番号kztと、特徴として選択されたかどうかを表すrk,iで分類することができる。最後に、観測値の確率モデルを定義する。観測値xi,tがモデルから実際に生成される確率の計算には、zt,rk,iで指定される観測パラメータが必要である。ztkかつrk,i1の場合には、特徴観測パラメータ37を使ってデータを生成する。一方ztkかつrk,i0の場合には、ノイズ観測パラメータ38を使ってデータを生成する。観測値がバイナリデータの場合、それぞれのパラメータは、以下の10式、11式のように決定される。特徴観測パラメータはクラスタkごとに特有の値を持つと考えられるので、10式で定式化される。一方、ノイズ観測パラメータはクラスタの番号によらないと考え、11式で定式化される。特徴観測パラメータ及びノイズ観測パラメータそれぞれのサンプリング結果を、記憶部3の対応する部分へと格納する。最終的な観測値xi,tの生成確率は、以下の12式のように評価される。以上の観測値を除く全ての変数を何らかの評価尺度の意味で最適な値に決定することが本クラスタリングシステムの主な目的である。結果として、観測量xi,tだけを観察することで、時系列データの生成パターンを、パターンを区別する上で重要な特徴だけを利用して理解、あるいは上記図1のように可視化して出力することが可能となる。<モデルのバリエーションについて>以上のモデルは観測値がバイナリデータの場合である。しかし、システムの構成をほとんど変えることなく、他の場合にも適用可能なシステムを実現できる。まず、観測値が非負整数の場合には、観測確率式である上記12式の関数として例えばポアソン分布が利用できる。その場合、関数の共役性から、上記10式、上記11式はそれぞれガンマ分布がふさわしい。それに応じてハイパーパラメータの個数や定義域を適切に変更すれば非負整数観測値に対応できる。観測値が実数の場合には、観測確率式である上記12式の関数として例えば正規分布が利用できる。その場合、関数の共役性から、上記10式、11式はそれぞれガウスガンマ分布あるいはガウス逆ガンマ分布がふさわしい。それに応じてハイパーパラメータの個数や定義域を適切に変更すれば実数観測値に対応できる。観測値がシンボルの場合には、観測確率式である上記12式の関数として例えば多項分布が利用できる。その場合、関数の共役性から上記10式、11式はそれぞれディリクレ分布がふさわしい。それに応じてハイパーパラメータの個数や定義域を適切に変更すればシンボル観測値に対応できる。以上の各関数やパラメータの設定については標準的な教科書、例えば、上記参考文献などに記載されている。<変数の最適化と出力>変数計算部22の目的は、観測データXが与えられた際に、記憶部3に格納された変数やパラメータを何らかの意味で最適な値に決定し、解析結果として出力することである。本発明では、上述のバイナリ観測値による定式化の場合を例にとって、効率の良い解析が可能とされているスライスサンプリング法上記非特許文献4を参照を利用した最適化法の実装方法を説明する。なお、各変数やパラメータの最適化には、モデルのバリエーションによらず他にも最尤推定法や事後確率最大化法、ギブスサンプリング法、変分ベイズ法、マルコフ連鎖モンテカルロ法、あるいは種々の数値最適化法などの様々な標準的な手法上記参考文献を参照を用いることも可能である。変数計算部22は、記憶部3に定義されている推定可能な変数やパラメータに対応する形で、各変数とパラメータの更新部分に分離されている。変数パラメータの更新の順番は任意であるが、通常は以下に説明する順番に一度ずつ更新を行い、一回の変数更新処理を完了する。<クラスタ割り当て変数更新部23>クラスタ割り当て変数更新部23では、クラスタ割り当て変数34を更新する。また、その過程で補助変数311も更新する。また、クラスタ割り当て変数Zの更新では潜在クラスタ数33も更新される。それに従って、Kに依存する変数やパラメータである特徴選択変数35、特徴観測パラメータ37、クラスタ遷移パラメータ39、混合比パラメータ310も更新される可能性がある。クラスタ割り当て変数更新部23のフローチャートを図4に示す。まず、ステップS121において、記憶部3から、TNKxzrθφπλβγαcdを読み込む。ステップS122において、所定の条件を満たすまで、補助変数ut,t=1,...,Tをサンプリングする。上記ステップS122は、図5に示すフローチャートにより実現される。ステップS151において、時刻tに初期値1を設定し、ステップS152で、以下の13式に示す確率でサンプリングして、補助変数utを決定する。ここでUniformは一様分布を表す。次のステップS153では、時刻tを1インクリメントとして、ステップS154において、時刻tがTより大きいか否かを判定する。時刻tがT以下である場合には、上記ステップS152へ戻る。一方、時刻tがTより大きい場合には、ステップS155へ進む。ステップS155では、時刻tに初期値1を設定し、ステップS156において、補助変数utについて、以下の条件を満たすか否かを判定する。ut≧πz4-1,K1上記の条件を満たす場合には、ステップS157において、時刻tを1インクリメントして、ステップS158において、時刻tがTより大きいか否かを判定する。時刻tがT以下である場合には、上記ステップS156へ戻る。一方、時刻tがTより大きい場合には、全てのtについてut≧πz4-1,K1であることを確認できたとして、上記ステップS122の処理を終了する。一方、上記の条件を満たさない場合には、次に示す手順でKを一つ増加させ、再度uをサンプリングする。ステップS159では、uのサンプリング結果を破棄し、ステップS160において、以下の式に従って、変数bをサンプリングする。b~Betaなお、b~Betaで、ベータ分布Betaより独立かつ一様にbをサンプリングすることを表す。そして、ステップS161において、以下の式に従って、β、βを計算して、βを拡張し更新する。β←βK1b,β←βK1-ββ←β1,β2,...,βK,β,β次のステップS162では、変数iに初期値1を設定し、ステップS163において、上記10式に従ってθK1,iをサンプリングして更新する。ステップS164では、変数iを1インクリメントして、ステップS165において、変数iがNより大きいか否かを判定する。変数iがN以下である場合には、上記ステップS163へ戻る。一方、変数iがNより大きい場合には、全てのiについてθK1,iをサンプリングしたとして、ステップS166へ進む。次のステップS166では、変数iに初期値1を設定し、ステップS167において、上記9式に従ってrK1,iをサンプリングして更新する。ステップS168では、変数iを1インクリメントして、ステップS169において、変数iがNより大きいか否かを判定する。変数iがN以下である場合には、上記ステップS167へ戻る。一方、変数iがNより大きい場合には、全てのiについてrK1,iをサンプリングしたとして、ステップS170へ進む。ステップS170では、潜在クラスタ数Kを、K1に更新する。そして、ステップS171において、クラスタ遷移パラメータ更新部28を用いてクラスタ遷移パラメータπk,lを更新して、上記ステップS151へ戻る。上記のように、utのサンプリングの後、全てのtについてut≧πzt-1,K1であることを確認した時点で補助変数utが決定する。この補助変数utは、理論上無限個存在する潜在クラスタを足切りして、有限のクラスタに限定する効果を持つ。より正確には、πkによるクラスタ遷移を可能性の高い遷移だけに限定している。また、上記図4のフローチャートのステップS123では、時刻tに初期値1を設定し、ステップS124で、変数kを初期値1に設定する。そして、ステップS125において、以下の14式に示す補助変数pt,kを計算する。ただし、x1N1tは、{xit}i1Nt1tを要素とするベクトルを表わす。この補助変数pt,kは、時刻tまでの観測値の情報を得たときに、時刻tがクラスタkに割り振られる確率を表す。この補助変数pt,kは、以下の15式に示す漸化式に従って計算される。ただし、∝は比例を表す。上記15式の右辺第一項は時刻tにおける全観測量の尤度をztkと仮定したときに計算した量である。すなわち、上記12式をt,ztkに固定してi1,...,Nまで掛け込んだ値である。これによって、各特徴量iが時刻tの所属するクラスタkに対して重要であるか否かに基づいて、補助変数ptkが算出される。ステップS126では、kを1インクリメントして、ステップS127において、変数kがKより大きいか否かを判定する。変数kがK以下である場合には、上記ステップS125へ戻る。一方、変数kがKより大きい場合には、全てのkについてptkを計算したとして、ステップS128へ進む。ステップS128では、上記ステップS125で計算されたpt,kをΣkpt,k1となるように正規化する。ステップS129では、時刻tを1インクリメントして、ステップS130において、時刻tがTより大きいか否かを判定する。時刻tがT以下である場合には、上記ステップS124へ戻る。一方、時刻tがTより大きい場合には、この補助変数をt=1,2,...,Tの順番で、全てのk=1,2,...,Kに対して計算したとして、ステップS131へ進む。ステップステップS131において、時刻tをTに設定する。そして、ステップS132において、変数kを初期値1に設定する。ステップS133では、zt=kとなる確率を以下の16式に従って計算する。上記16式の意味するところは、全時刻の観測量と1時刻後のクラスタ割り当てzt1が与えられた時に時刻tのクラスタ割り当てがkとなる確率である。Iは、条件cが成立した場合に1そうでない場合には0を返す関数である。なお、最初にtTを決定する際には、zT1の値を便宜的に固定しておくか、あるいは以下の17式で計算する。次のステップS134では、変数kを1インクリメントして、ステップS135において、変数kがKより大きいか否かを判定する。変数kがK以下である場合には、上記ステップS133へ戻る。一方、変数kがKより大きい場合には、全てのkについてzt=kとなる確率を計算したとして、ステップS136へ進む。そして、ステップS136では、上記ステップS133で計算された全てのkについてzt=kとなる確率上記16式の右辺の量に比例する確率でztの値をk∈1,2,...,Kより一つ選択するようにサンプリングする。次のステップS137では、時刻tを1デクリメントして、ステップS138において、時刻tが1より小さいか否かを判定する。時刻tが1以上である場合には、上記ステップS132へ戻る。一方、時刻tが1より小さい場合には、全てのtについてzrtをサンプリングしたとして、ステップS139へ進む。このように、tT,T1,...,1の順番でztの値を決定する。続いて、ステップS139において、サンプリングによって不要となったクラスタの整理を行う。z1Tをサンプリングした結果、ztkとなる時刻インデックスtが一つもないクラスタが発生しうる。これは、今回の推論の結果、クラスタkはデータを表現する上で不要であると判断されたことを意味する。そこで、そのようなクラスタは全てメモリから消去し、同時にKに依存する変数を全て更新する。例えば、utのサンプリング後にK5となり、ztのサンプリング後にk3となるztが一つもなかったとする。このとき、以下の処理を行い、潜在クラスタ数を一つ減らす。1.クラスタ割り当て変数Zのうち、zt3となる全てのztの値を1小さくして更新2.特徴選択変数Rのうち、k3となるrk,iの値をrk1,iに上書きして更新し、rK,iはメモリから消去する3.特徴観測パラメータθk,iのうち、k3となるθk,iの値をθk1,iに上書きして更新し、θK,iはメモリから消去する4.クラスタ遷移パラメータπk,lのうち、k3となるkπ,lの値をπk1,lに上書きして更新し、l3となるπk,l1の値をπk,l1に上書きして更新し、θK1,lθk,K1は全てメモリから消去する5.混合比パラメータβkのうち、k3となるβkの値をβk1に上書きして更新し、βK←1-ΣK1k1βkと更新し、βKK1はメモリから消去する6.K←K1以上と同等の処理を、使われていないクラスタ全てに対して繰り返す。上記ステップS139におけるクラスタの整理手順のフローチャートを図6に示す。ステップS180において、変数kをKに設定し、ステップS181において、時刻tを初期値1に設定する。ステップS182では、zt=kであるか否かを判定する。zt=kである場合には、クラスタkは使われていると判断し、後述するステップS211へ移行する。一方、zt=kでない場合には、ステップS183において、時刻tを1インクリメントして、ステップS184において、時刻tがTより大きいか否かを判定する。時刻tがT以下である場合には、上記ステップS182へ戻る。一方、時刻tがTより大きい場合には、クラスタkが使われないと判断し、ステップS185へ進む。ステップS185では、kにkを設定し、ステップS186において、時刻tに1を設定する。ステップS187において、zt>kであるか否かを判定し、zt>kであれば、ステップS188で、ztkとなるztの値を1小さくして更新する。一方、zt>kでない場合には、ステップS189へ移行する。ステップS189では、時刻tを1インクリメントして、ステップS190において、時刻tが、時刻tがTより大きいか否かを判定する。時刻tがT以下である場合には、上記ステップS186へ戻る。一方、時刻tがTより大きい場合には、ztkとなるztの値を更新したとして、ステップS191へ進む。ステップS191では、変数kにkを設定し、ステップS192において、i=1,.,Nに対しrk,iの値をrk1,iに上書きして更新する。ステップS193では、i=1,.,Nに対しθk,iの値をθk1,iに上書きして更新する。そして、ステップS194において、変数kを1インクリメントして、ステップS196で、変数kがKに達したか否かを判定する。変数kがKに達していない場合には、ステップS192へ戻るが、一方、変数kがKに達した場合には、kkとなるrk,iの値及びθk,iの値を更新したとして、ステップS197へ移行する。ステップS197では、i=1,.,Nに対しrK,iをメモリから消去する。次のステップS198では、i=1,.,Nに対し、θK,iをメモリから消去する。そして、ステップS199では、変数kにkを設定し、ステップS200において、l=1,.,K1に対しπk,lの値をπk1,lに上書きして更新する。ステップS201では、l=1,.,K1に対しπl,kの値をπl1,kに上書きして更新する。ステップS202において、βkの値をβk1に上書きして更新する。ステップS203では、kを1インクリメントして、ステップS204において、変数kがK1に達したか否かを判定する。変数kがK1に達していない場合には、ステップS200へ戻るが、一方、変数kがK1に達した場合には、kkとなるπl1,kの値、πk,l1の値、及びβk1の値を更新したとして、ステップS205へ移行する。ステップS205では、l=1,.,K1に対しπK1,lをメモリから消去する。次のステップS206では、l=1,.,K1に対し、πl,K1をメモリから消去する。ステップS207では、βK←1-ΣK1k1βkと更新し、ステップS208において、βKK1はメモリから消去する。次のステップS209において、潜在クラスタ数Kを1つ減らして更新し、ステップS201において、kにkを設定して、ステップS211へ移行する。ステップS211では、変数kを1デクリメントして、ステップS212において、変数kが、1より小さいか否かを判定する。変数kが1以上である場合には、上記ステップS181へ戻る。一方、変数kが1より小さい場合には、使われていないクラスタ全てを整理したとして、上記ステップS139の処理を終了する。上記図4のフローチャートのステップS140では、変数kを1に設定し、ステップS141において、変数lを1に設定する。そして、ステップS142において、以下の18式に従って、補助変数mk,lを計算する。なお、mk,lは全時刻を通じてクラスタkからクラスタlに遷移した回数を表す。そして、ステップS143では、補助変数Sk,lを確率的に以下の19式でサンプリングする。なお、Sk,lは混合比パラメータ310の更新などに必要となる。ここでsはStirlingnumberofthefirstkindと呼ばれる関数である.次のステップS144では、変数lを1インクリメントして、ステップS145において、変数lがKより大きいか否かを判定する。変数lがK以下である場合には、上記ステップS142へ戻る。一方、変数lがKより大きい場合には、ステップS146において、変数kを1インクリメントする。そして、ステップS147において、変数kがKより大きいか否かを判定する。変数kがK以下である場合には、上記ステップS141へ戻る。変数kがKより大きい場合には、1≦k,l≦Kについて、補助変数mk,l、Sk,lを計算したとして、ステップS148へ進む。ステップS148では、更新された全ての変数、パラメータ、補助変数を記憶部3の対応する部分に設定して、クラスタ割り当て変数の更新が終了する。<特徴選択変数更新部24>特徴選択変数更新部24では、特徴選択変数35を更新する。特徴選択変数更新部24の更新手順のフローチャートを図7に示す。まず、ステップS222において、記憶部3により、TNKxzrφθefをロードする。そして、ステップS223において、変数kを1に設定し、ステップS224において、変数iを1に設定する。そして、ステップS225において、以下の20式、21式を利用してrk,iの値を0か1かに更新する。上記の式で利用されている記号は以下のとおり。上記20式の右辺は、特徴量iがクラスタkに対して重要でないと仮定した場合の適合度に対応し、上記21式の右辺は、特徴量iがクラスタkに対して重要であると仮定した場合の適合度に対応している。そして、ステップS226において、変数iを1インクリメントし、ステップS227において、変数iがNより大きいか否かを判定する。変数iがN以下である場合には、上記ステップS225へ戻る。一方、変数iがNより大きい場合には、全てのiについてrk,iをサンプリングしたとして、ステップS228へ進む。ステップS228では、変数kを1インクリメントし、ステップS229において、変数kがKより大きいか否かを判定する。変数kがK以下である場合には、上記ステップS224へ戻る。一方、変数kがKより大きい場合には、k1,...,K、i1,...,Nの全てについてrk,iをサンプリングしたとして、ステップS230へ進む。ステップS230では、更新されたrk,iを記憶部3の特徴選択変数35に格納し、更新を終了する。<特徴選択パラメータ更新部25>特徴選択パラメータ更新部25では、特徴選択パラメータ36を更新する。特徴選択パラメータ更新部25の更新手順のフローチャートを図8に示す。まず、ステップS231において、記憶部3により、NKrefをロードする。そして、ステップS232において、変数iを1に設定し、ステップS233において以下の26式を利用してλiの値をサンプリングして更新する。そして、ステップS234において、変数iを1インクリメントし、ステップS235において、変数iがNより大きいか否かを判定する。変数iがN以下である場合には、上記ステップS233へ戻る。一方、変数iがNより大きい場合には、i1,...,Nの全てについてλiの値を更新したとして、ステップS236へ進む。ステップS236では、更新されたλiを記憶部3の特徴選択パラメータ36に格納し、更新を終了する。<特徴観測パラメータ更新部26>特徴観測パラメータ更新部26では、特徴観測パラメータ37を更新する。特徴観測パラメータ更新部26の更新手順のフローチャートを図9に示す。まず、ステップS241において、記憶部3により、TNKxzrcdをロードする。そして、ステップS242において、変数kを1に設定し、ステップS243において、変数iを1に設定する。そして、ステップS244において、以下の27式を利用してθk,iの値をサンプリングして更新する。ただし、nk,i,Nk,iはそれぞれ上記24式、25式で定義されたものと同様である。そして、ステップS245において、変数iを1インクリメントし、ステップS246において、変数iがNより大きいか否かを判定する。変数iがN以下である場合には、上記ステップS244へ戻る。一方、変数iがNより大きい場合には、全てのiについてθk,iをサンプリングしたとして、ステップS247へ進む。ステップS247では、変数kを1インクリメントし、ステップS248において、変数kがKより大きいか否かを判定する。変数kがK以下である場合には、上記ステップS243へ戻る。一方、変数kがKより大きい場合には、k1,...,K、i1,...,Nの全てについてθk,iをサンプリングしたとして、ステップS249へ進む。ステップS249では、更新されたθk,iを記憶部3の特徴観測パラメータ37に格納し、更新を終了する。<ノイズ観測パラメータ更新部27>ノイズ観測パラメータ更新部27では、ノイズ観測パラメータ38を更新する。ノイズ観測パラメータ更新部27の更新手順のフローチャートを図10に示す。まず、ステップS251において、記憶部3により、NKxzrabをロードする。そして、ステップS252において、変数iを1に設定し、ステップS253において以下の28式~30式を利用して、φiの値をサンプリングして更新する。ただし、nk,i,Nk,iはそれぞれ上記24式、25式で定義されたものと同様である。そして、ステップS254において、変数iを1インクリメントし、ステップS255において、変数iがNより大きいか否かを判定する。変数iがN以下である場合には、上記ステップS253へ戻る。一方、変数iがNより大きい場合には、i1,...,Nの全てについてφiの値を更新したとして、ステップS256へ進む。ステップS256では、更新されたφiを記憶部3のノイズ観測パラメータ38に格納し、更新を終了する。<クラスタ遷移パラメータ更新部28>クラスタ遷移パラメータ更新部28では、クラスタ遷移パラメータ39を更新する。クラスタ遷移パラメータ更新部28の更新手順のフローチャートを図11に示す。まず、ステップS261において、記憶部3により、Kαβmをロードする。そして、ステップS262において、変数kを1に設定し、ステップS263において、補助変数mk,lを利用して、以下の31式にしたがって、πkの値をサンプリングして更新する。そして、ステップS264において、変数kを1インクリメントし、ステップS265において、変数kがK1より大きいか否かを判定する。変数kがK1以下である場合には、上記ステップS263へ戻る。一方、変数kがK1より大きい場合には、k1,...,K1の全てについてπkの値を更新したとして、ステップS266へ進む。ステップS266では、更新されたπkを記憶部3のクラスタ遷移パラメータ39に格納し、更新を終了する。<混合比パラメータ更新部29>混合比パラメータ更新部29では、混合比パラメータ310を更新する。混合比パラメータ更新部29の更新手順のフローチャートを図12に示す。まず、ステップS271において、記憶部3により、KγSをロードする。そして、ステップS272において、補助変数Sk,lを利用して、以下の33式に従って、K1次元ベクトルをサンプリングして更新する。ステップS273では、更新されたβを記憶部3の混合比パラメータ310に格納し、更新を終了する。<ハイパーパラメータ更新部210>ハイパーパラメータ更新部210では、ハイパーパラメータ312の各要素を更新する。ここでは、超事前分布としてガンマ分布を仮定し、その事後分布からハイパーパラメータをサンプリングする方法を記載する。なお、以下は簡単のため超事前分布のガンマ分布のパラメータは省略しGammaとのみ表記する。ハイパーパラメータ更新部210の更新手順のフローチャートを図13に示す。まず、ステップS281において、記憶部3により、TKNrθφγαβabcdefSをロードする。そして、ステップS282において、以下の34式を利用して、γの値をサンプリングして更新する。次のステップS283では、以下の35式を利用して、αの値をサンプリングして更新する。そして、ステップS284において、変数iを1に設定し、ステップS285において以下の36式を利用してaiの値をサンプリングして更新する。ここでBはベータ関数である。次のステップS286において以下の37式を利用してbiの値をサンプリングして更新する。ステップS287では、以下の38式を利用してciの値をサンプリングして更新する。ステップS288では、以下の39式を利用してdiの値をサンプリングして更新する。ステップS289では、以下の40式を利用してeiの値をサンプリングして更新する。ステップS290では、以下の41式を利用してfiの値をサンプリングして更新する。そして、ステップS291において、変数iを1インクリメントし、ステップS292において、変数iがNより大きいか否かを判定する。変数iがN以下である場合には、上記ステップS285へ戻る。一方、変数iがNより大きい場合には、i1,...,Nの全てについてai,bi,ci,di,ei,fiの値を更新したとして、ステップS293へ進む。ステップS293では、更新されたハイパーパラメータを記憶部3のハイパーパラメータ312に格納し、更新を終了する。なお、ハイパーパラメータの更新には、他にも様々な方法を利用することができる。また、ハイパーパラメータは必ずしも更新しなくても良い。<解析終了条件判定部211>解析終了条件判定部211は、変数計算部22を監視し、変数の更新処理を続けるか終了するかを判定する。一般に、変数の更新は繰り返し処理により最適値へ近づく。しかし、いつ繰り返しを止めるかについては様々な方法がある。例えば、変数更新後に記憶部3に格納された変数やパラメータを使って、与えられたデータの尤度を計測しそれが一定値あるいは一定割合以上の変化をしなくなったら終了とする、あるいは単純に変数計算部22の繰り返し回数が一定値以上になったら更新を終了するなどである。以下の実験結果では、最も簡単に、Lを繰り返し処理の最大値として入力し、変数計算部の更新処理をL回繰り返した時点で解析処理を終了することとした。実際にはどのような終了条件を決定しても構わない。設定した解析終了条件が満たされていない場合、再び変数計算部22に処理を戻す。設定した解析終了条件が満たされると、解析結果はユーザが規程した所望の形で出力部4へと出力される。<実験例>次に、本実施の形態で提案する時系列データのクラスタリング手法を用いた実験の結果について説明する。実験では、3種類、計5つのデータセットに対して実験を行った。データの種類は以下の3種類である。1.人工データ提案法の想定する時系列データの生成過程に従い、人工的に二値のバイナリデータを生成した。synth1はN200、T50、synth2はN100、T50、synth3はN100、T70であり、それぞれ想定する構造に違いがある。2文書データPrceedingsoftheNationalAcademyofSciencesと呼ばれる学術論文誌のデータを利用した。1915年から2005年までの医学薬学分野の論文タイトルに出現する単語の頻度を年度ごとに集計し行列データとして解析する。年度tで50回以上出現した単語iの観測量をxi,t1とし、50回に満たない単語はxi,t0とする。T91年に渡って一度でもxi,t1となるN118種の単語の01データである。3.音楽データおよそ11秒の合成音楽の波形を解析した。音楽波形に対して移動窓を使った短時間フーリエ解析を行い、T742個フレームにおいてN513の周波数帯におけるスペクトル値を算出する。xi,tとしては、フレームt、周波数帯iにおけるスペクトル値の自然対数を計算する。対数スペクトル値が一定の閾値を越えていればxi,t1、そうでなければxi,t0とした。定量評価として、テストデータ対数尤度を用いた。与えられたN×T個の観測データのうちおよそ5の観測データを未知のテストデータとしてランダムに選択肢、変数更新部の処理から除外する。解析終了後、これらテストデータの対数尤度を計算し、それをテストデータ数で平均した。比較対象として、特徴選択変数rを含まない、通常のiHMM上記非特許文献3、4を参照を利用した。提案法および比較対象であるiHMMともに複数のハイパーパラメータの組み合わせでそれぞれ20回試行を行い、最も数値の良かったハイパーパラメータを選択した。実験結果を以下の表1にまとめた。対数尤度は数字が大きいほど良い。この表1より分かるように、特徴選択の機能が加わった提案法は大部分のデータで既存手法より正確に時系列データのモデル化が可能であることが示された。最後に定性的な評価としてPNASデータの解析によって得られた、年度ごとの特徴単語を図14に示す。各コラムはクラスタを表わし、クラスタ上部の数字はクラスタの対応する年度時間ステップtであり、下部の単語列は特徴として選ばれた単語である。これによれば、例えばvitroという単語は1980年ごろに現れ1996年ごろには論文のタイトルとして現れなくなった。それと入れかわって、vivoというという単語が1991年ごろに現れた。これらはinvitroとinvivoという生化学医学研究における重要な概念の出現と変化に符合している。また、分子生物学を表すと思われるmolecularという単語は1985年ごろに現れさらに遺伝子学や遺伝子治療に関連する。geneticという単語は1996年ごろに現れる。このように、各年代で大きく話題になったトピックの単語を正しく各クラスタの特徴として抽出できている.これにより、例えば医学生理学分野の技術トレンドを少数のクラスタの選択された特徴単語から簡単に調査することができる。以上説明したように、本発明の実施の形態に係るクラスタリング装置によれば、クラスタ割り当て変数ztを更新するときに潜在クラスタ数Kを決定すると共に、クラスタ割り当て変数ztの更新、特徴選択変数rkiの更新、特徴選択パラメータの更新λi、特徴観測パラメータの更新θki、ノイズ観測パラメータの更新φi、クラスタ遷移パラメータπklの更新、混合比パラメータβk、及びハイパーパラメータの更新を繰り返すことにより、観測値の時系列データXが与えられたときに、クラスタ数を自動的に決定すると共に、時刻のクラスタリングに寄与する真の特徴量と、時刻のクラスリングに寄与しないノイズ特徴の影響を自動的に推定分離して、時刻のクラスタリングを行うことができる。また、高次元の特徴データを複数の似た特性を示す時間帯へとクラスタリングし、複雑な時間変化を少数のパターンの変化として直感的に解析し、データのトレンドを簡単に調査可視化できる。また、既存の特徴選択を含まない時系列データモデリングよりも高いデータ表現精度を持つことができる。また、与えられた観測量の中から、時刻のクラスタリングに寄与する真の特徴量と、クラスタ化に寄与しないノイズ特徴量の存在をモデル化する。これによって、与えられた観測値の時系列データをクラスタリングする上で意味のある情報と意味の無い情報を分離することができるので、クラスタの推定時にノイズ特徴の影響を排除することが可能となる。また、真の特徴量とノイズ特徴量は各クラスタごとに固有に定義される。このことで、動的に特徴の分布が変化する時系列データに対してより精度の良いモデル化が可能となる。また、各時刻を特徴づける特徴量を自動的に選択抽出するので、与えられた時系列データが高次元の特徴量を持っていてもその定性的な性質を容易に理解することが可能となる。また、観測量を重要な特徴、ノイズ観測量へと分類することができるので、各時刻において本質的に違いを生み出した情報を自動的に抽出して提供可能である。また、クラスタの数は自動的にデータに応じて決定される。これまで多くの手法では探索的に決定することが多かった。このため、手法利用時の工数を減らすことができる。なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。例えば、上述のクラスタリング装置は、内部にコンピュータシステムを有しているが、コンピュータシステムは、WWWシステムを利用している場合であれば、ホームページ提供環境あるいは表示環境も含むものとする。また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。"