(57)【特許請求の範囲】
【請求項１】
各時刻に観測された複数の特徴量からなる観測値の時系列データに基づいて、時刻のクラスタリングを実行するクラスタリング装置であって、
各時刻ｔにおける各特徴量ｉからなる観測値ｘi,t、潜在クラスタ数Ｋ、時刻ｔが所属するクラスタを示すクラスタ割り当て変数ｚt、各クラスタｋに対して特徴量ｉが重要であるか否かを示す特徴選択変数ｒk,i、あるクラスタに対して特徴量ｉが重要である確率を示す特徴選択パラメータλi、特徴量ｉが時刻ｔの所属するクラスタｋに対して重要である場合に観測値ｘi,tを生成するために用いられる特徴観測パラメータθk,i、特徴量ｉが時刻ｔの所属するクラスタに対して重要でない場合に観測値ｘi,tを生成するために用いられるノイズ観測パラメータφi、時刻ｔ−１がクラスタｋに所属し、かつ、次の時刻ｔがクラスタｌに所属する確率を示すクラスタ遷移パラメータπk,ｌ、各クラスタｋに所属する時刻の割合を示す混合比パラメータβk、補助変数ｕt、ｐt,k及びハイパーパラメータγ、α、ａi、ｂi、ｃi、ｄi、ｅi、ｆiを記憶する記憶部と、
前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴選択パラメータλi、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπk,i、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ｃi、ｄiを取得し、各時刻の前記補助変数ｕtを決定すると共に前記潜在クラスタ数Ｋを決定し、各時刻ｔ及び各クラスタｋの前記補助変数ｐt,kを算出して、各時刻ｔの前記クラスタ割り当て変数ｚｔをサンプリングにより算出し、前記決定した前記補助変数ｕt及び前記潜在クラスタ数Ｋ、前記算出した前記補助変数ｐt,k、並びに前記算出した前記クラスタ割り当て変数ｚtによって、前記記憶部に記憶してある前記補助変数ｕt、前記潜在クラスタ数Ｋ、前記補助変数ｐt,k、及び前記クラスタ割り当て変数ｚtを更新し記憶するクラスタ割り当て変数更新手段と、
前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタｋ及び各特徴量ｉの前記特徴選択変数ｒk,iをサンプリングにより算出し、前記算出した前記特徴選択変数ｒk,iによって前記記憶部に記憶してある前記特徴選択変数ｒk,iを更新し記憶する特徴選択変数更新手段と、
前記記憶部から前記潜在クラスタ数Ｋ、前記特徴選択変数ｒk,i、及び前記ハイパーパラメータei、fiを取得し、各特徴量ｉの前記特徴選択パラメータλiをサンプリングにより算出し、前記算出した前記特徴選択パラメータλiによって前記記憶部に記憶してある前記特徴選択パラメータλiを更新し記憶する特徴選択パラメータ更新手段と、
前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、及び前記ハイパーパラメータｃi、ｄiを取得し、各クラスタｋ及び各特徴量ｉの前記特徴観測パラメータθk,iをサンプリングにより算出し、前記算出した前記特徴観測パラメータθk,iによって前記記憶部に記憶してある前記特徴観測パラメータθk,iを更新し記憶する特徴観測パラメータ更新手段と、
前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、及び前記ハイパーパラメータａi、ｂiを取得し、各特徴量ｉの前記ノイズ観測パラメータφiをサンプリングにより算出し、前記算出した前記ノイズ観測パラメータφiによって前記記憶部に記憶してある前記ノイズ観測パラメータφiを更新し記憶するノイズ観測パラメータ更新手段と、
前記記憶部から前記潜在クラスタ数Ｋ、前記混合比パラメータβk、及び前記ハイパーパラメータαを取得し、各クラスタｋ、ｌの前記クラスタ遷移パラメータπk,lをサンプリングにより算出し、前記算出した前記クラスタ遷移パラメータπk,lによって前記記憶部に記憶してある前記クラスタ遷移パラメータπk,lを更新し記憶するクラスタ遷移パラメータ更新手段と、
前記記憶部から前記潜在クラスタ数Ｋ、及び前記ハイパーパラメータγを取得し、各クラスタｋの前記混合比パラメータβkをサンプリングにより算出し、前記算出した前記混合比パラメータβkによって前記記憶部に記憶してある前記混合比パラメータβkを更新し記憶する混合比パラメータ更新手段と、
前記クラスタ割り当て変数更新手段、前記特徴選択変数更新手段、前記特徴選択パラメータ更新手段、前記特徴観測パラメータ更新手段、前記ノイズ観測パラメータ更新手段、前記クラスタ遷移パラメータ更新手段、及び前記混合比パラメータ更新手段による更新を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す終了条件判定手段と、
を含むクラスタリング装置。
【請求項２】
前記クラスタ割り当て変数更新手段は、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴選択パラメータλi、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπk,l、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ｃi、ｄiを取得し、各時刻ｔの前記補助変数ｕtを決定すると共に、前記潜在クラスタ数Ｋを決定し、各時刻ｔについて、各特徴量ｉが時刻ｔの所属するクラスタｋに対して重要であるか否かに基づいて前記補助変数ｐt,kを算出し、前記補助変数ｕt、ｐt,kを用いてサンプリングして前記クラスタ割り当て変数ｚtの値を算出し、
前記特徴選択変数更新手段は、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタｋ及び各特徴量ｉについて、特徴量ｉがクラスタｋに対して重要でないと仮定した場合の適合度及び特徴量ｉがクラスタｋに対して重要であると仮定した場合の適合度を用いてサンプリングして前記特徴選択変数ｒk,iを算出する請求項１記載のクラスタリング装置。
【請求項３】
前記クラスタ割り当て変数更新手段は、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴選択パラメータλi、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπk,l、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ｃi、ｄiを取得し、各時刻の前記補助変数ｕtを決定すると共に、前記潜在クラスタ数Ｋを決定し、時刻ｔ＝１，２，．．．，Ｔの順番で、全てのｋ＝１，２，．．．，Ｋに対して以下の（Ｉ）式に示す漸化式に従って、前記補助変数ｐt,kを算出し、時刻ｔ＝Ｔ，Ｔ−１，．．．，１の順番で、以下の（ＩＩ）式に示す各クラスタｋに対する確率ｐ（ｚt＝ｋ｜ｚt+1，Ｘ、Ｒ、Φ、Θ、Π、Ｕ）からサンプリングして前記クラスタ割り当て変数ｚtの値を算出し、
前記特徴選択変数更新手段は、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタｋ及び各特徴量ｉについて、以下の（ＩＩＩ）式に示す確率ｐ（ｒk,i＝０｜Ｘ、Ｚ、Ｒ＼(k,i)、Λ、Φ、Θ）及び以下の（ＩＶ）式に示す確率ｐ（ｒk,i＝１｜Ｘ、Ｚ、Ｒ＼(k,i)、Λ、Φ、Θ）からサンプリングして前記特徴選択変数ｒk,iを算出する請求項１記載のクラスタリング装置。
【数１】
0005851205.tif 000036
ここでBernoulli(・) はベルヌーイ分布を表す。
【数２】
0005851205.tif 000037
I(c)は、条件cが成立した場合に1,そうでない場合には0を返す関数である。
【数３】
0005851205.tif 000038
【請求項４】
前記記憶部から前記潜在クラスタ数Ｋ、前記特徴選択変数ｒk,i、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、前記混合比パラメータβk、及び前記ハイパーパラメータγ、α、ａi、ｂi、ｃi、ｄi、ｅi、ｆiを取得し、ハイパーパラメータγをサンプリングにより算出し、ハイパーパラメータαをサンプリングにより算出し、各特徴量ｉのハイパーパラメータａiをサンプリングにより算出し、各特徴量ｉのハイパーパラメータｂiをサンプリングにより算出し、各特徴量ｉのハイパーパラメータｃiをサンプリングにより算出し、各特徴量ｉのハイパーパラメータｄiをサンプリングにより算出し、各特徴量ｉのハイパーパラメータｅiをサンプリングにより算出し、各特徴量ｉのハイパーパラメータｆiをサンプリングにより算出し、前記算出した前記ハイパーパラメータγ、α、ａi、ｂi、ｃi、ｄi、ｅi、ｆiによって前記記憶部に記憶してある前記ハイパーパラメータγ、α、ａi、ｂi、ｃi、ｄi、ｅi、ｆiを更新し記憶するハイパーパラメータ更新手段を更に含み、
前記終了条件判定手段は、前記クラスタ割り当て変数更新手段、前記特徴選択変数更新手段、前記特徴選択パラメータ更新手段、前記特徴観測パラメータ更新手段、前記ノイズ観測パラメータ更新手段、前記クラスタ遷移パラメータ更新手段、前記混合比パラメータ更新手段、及び前記ハイパーパラメータ更新手段による更新を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す請求項１〜請求項３の何れか１項記載のクラスタリング装置。
【請求項５】
各時刻に観測された複数の特徴量からなる観測値の時系列データに基づいて、時刻のクラスタリングを実行するクラスタリング装置におけるクラスタリング方法であって、
前記クラスタリング装置は、
各時刻ｔにおける各特徴量ｉからなる観測値ｘi,t、潜在クラスタ数Ｋ、時刻ｔが所属するクラスタを示すクラスタ割り当て変数ｚt、各クラスタｋに対して特徴量ｉが重要であるか否かを示す特徴選択変数ｒk,i、あるクラスタに対して特徴量ｉが重要である確率を示す特徴選択パラメータλi、特徴量ｉが時刻ｔの所属するクラスタｋに対して重要である場合に観測値ｘi,tを生成するために用いられる特徴観測パラメータθk,i、特徴量ｉが時刻ｔの所属するクラスタに対して重要でない場合に観測値ｘi,tを生成するために用いられるノイズ観測パラメータφi、時刻ｔ−１がクラスタｋに所属し、かつ、次の時刻ｔがクラスタｌに所属する確率を示すクラスタ遷移パラメータπk,ｌ、各クラスタｋに所属する時刻の割合を示す混合比パラメータβk、補助変数ｕt、ｐt,k及びハイパーパラメータγ、α、ａi、ｂi、ｃi、ｄi、ｅi、ｆiを記憶する記憶部と解析部とを備え、
前記解析部は、
前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴選択パラメータλi、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπk,i、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ｃi、ｄiを取得し、各時刻の前記補助変数ｕtを決定すると共に前記潜在クラスタ数Ｋを決定し、各時刻ｔ及び各クラスタｋの前記補助変数ｐt,kを算出して、各時刻ｔの前記クラスタ割り当て変数ｚｔをサンプリングにより算出し、前記決定した前記補助変数ｕt及び前記潜在クラスタ数Ｋ、前記算出した前記補助変数ｐt,k、並びに前記算出した前記クラスタ割り当て変数ｚtによって、前記記憶部に記憶してある前記補助変数ｕt、前記潜在クラスタ数Ｋ、前記補助変数ｐt,k、及び前記クラスタ割り当て変数ｚtを更新し記憶するクラスタ割り当て変数更新ステップと、
前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタｋ及び各特徴量ｉの前記特徴選択変数ｒk,iをサンプリングにより算出し、前記算出した前記特徴選択変数ｒk,iによって前記記憶部に記憶してある前記特徴選択変数ｒk,iを更新し記憶する特徴選択変数更新ステップと、
前記記憶部から前記潜在クラスタ数Ｋ、前記特徴選択変数ｒk,i、及び前記ハイパーパラメータei、fiを取得し、各特徴量ｉの前記特徴選択パラメータλiをサンプリングにより算出し、前記算出した前記特徴選択パラメータλiによって前記記憶部に記憶してある前記特徴選択パラメータλiを更新し記憶する特徴選択パラメータ更新ステップと、
前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、及び前記ハイパーパラメータｃi、ｄiを取得し、各クラスタｋ及び各特徴量ｉの前記特徴観測パラメータθk,iをサンプリングにより算出し、前記算出した前記特徴観測パラメータθk,iによって前記記憶部に記憶してある前記特徴観測パラメータθk,iを更新し記憶する特徴観測パラメータ更新ステップと、
前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、及び前記ハイパーパラメータａi、ｂiを取得し、各特徴量ｉの前記ノイズ観測パラメータφiをサンプリングにより算出し、前記算出した前記ノイズ観測パラメータφiによって前記記憶部に記憶してある前記ノイズ観測パラメータφiを更新し記憶するノイズ観測パラメータ更新ステップと、
前記記憶部から前記潜在クラスタ数Ｋ、前記混合比パラメータβk、及び前記ハイパーパラメータαを取得し、各クラスタｋ、ｌの前記クラスタ遷移パラメータπk,lをサンプリングにより算出し、前記算出した前記クラスタ遷移パラメータπk,lによって前記記憶部に記憶してある前記クラスタ遷移パラメータπk,lを更新し記憶するクラスタ遷移パラメータ更新ステップと、
前記記憶部から前記潜在クラスタ数Ｋ、及び前記ハイパーパラメータγを取得し、各クラスタｋの前記混合比パラメータβkをサンプリングにより算出し、前記算出した前記混合比パラメータβkによって前記記憶部に記憶してある前記混合比パラメータβkを更新し記憶する混合比パラメータ更新ステップと、
前記クラスタ割り当て変数更新ステップ、前記特徴選択変数更新ステップ、前記特徴選択パラメータ更新ステップ、前記特徴観測パラメータ更新ステップ、前記ノイズ観測パラメータ更新ステップ、前記クラスタ遷移パラメータ更新ステップ、及び前記混合比パラメータ更新ステップによる更新を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す終了条件判定ステップと、
を実行することを特徴とするクラスタリング方法。
【請求項６】
前記クラスタ割り当て変数更新ステップは、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴選択パラメータλi、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπk,l、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ｃi、ｄiを取得し、各時刻ｔの前記補助変数ｕtを決定すると共に、前記潜在クラスタ数Ｋを決定し、各時刻ｔについて、各特徴量ｉが時刻ｔの所属するクラスタｋに対して重要であるか否かに基づいて前記補助変数ｐt,kを算出し、前記補助変数ｕt、ｐt,kを用いてサンプリングして前記クラスタ割り当て変数ｚtの値を算出し、
前記特徴選択変数更新ステップは、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタｋ及び各特徴量ｉについて、特徴量ｉがクラスタｋに対して重要でないと仮定した場合の適合度及び特徴量ｉがクラスタｋに対して重要であると仮定した場合の適合度を用いてサンプリングして前記特徴選択変数ｒk,iを算出する請求項５記載のクラスタリング方法。
【請求項７】
前記記憶部から前記潜在クラスタ数Ｋ、前記特徴選択変数ｒk,i、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、前記混合比パラメータβk、及び前記ハイパーパラメータγ、α、ａi、ｂi、ｃi、ｄi、ｅi、ｆiを取得し、ハイパーパラメータγをサンプリングにより算出し、ハイパーパラメータαをサンプリングにより算出し、各特徴量ｉのハイパーパラメータａiをサンプリングにより算出し、各特徴量ｉのハイパーパラメータｂiをサンプリングにより算出し、各特徴量ｉのハイパーパラメータｃiをサンプリングにより算出し、各特徴量ｉのハイパーパラメータｄiをサンプリングにより算出し、各特徴量ｉのハイパーパラメータｅiをサンプリングにより算出し、各特徴量ｉのハイパーパラメータｆiをサンプリングにより算出し、前記算出した前記ハイパーパラメータγ、α、ａi、ｂi、ｃi、ｄi、ｅi、ｆiによって前記記憶部に記憶してある前記ハイパーパラメータγ、α、ａi、ｂi、ｃi、ｄi、ｅi、ｆiを更新し記憶するハイパーパラメータ更新ステップを更に含み、
前記終了条件判定ステップは、前記クラスタ割り当て変数更新ステップ、前記特徴選択変数更新ステップ、前記特徴選択パラメータ更新ステップ、前記特徴観測パラメータ更新ステップ、前記ノイズ観測パラメータ更新ステップ、前記クラスタ遷移パラメータ更新ステップ、前記混合比パラメータ更新ステップ、及び前記ハイパーパラメータ更新ステップによる更新を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す請求項５又は６記載のクラスタリング方法。
【請求項８】
請求項１〜請求項４の何れか１項に記載のクラスタリング装置を構成する各手段として、コンピュータを機能させることを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、クラスタリング装置、方法、及びプログラムに係り、特に、複数の特徴量からなる観測値の時系列データに基づいて、時刻のクラスタリングを実行するクラスタリング装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
時系列データの解析は古くから多くの分野で必要とされている技術である。メディアデータにおいては携帯電話上の音声認識アプリケーションや3D動作認識によるゲーム操作器、為替データではロボットプログラムによる市況解析と自動株式売買など、実応用にいても多数利用されている。
【０００３】
また、近年ではWeb 上に大量のデジタルデータが存在している。それらの中でも時系列データはそのリッチな情報と複雑さから解析技術が強く望まれているデータである。例えばWeb上の新聞やブログデータの単語の特徴の時間推移を測ることで話題の変遷を推測したり、あるいは様々な商品の購買履歴データから売れ筋の商品の変化と今後の予測を行う、などが応用として考えることができる。
【０００４】
これまで、時系列データの解析方法として様々な手法が開発されてきたが、近年最も成功を収めているのが隠れマルコフモデル(HMM)である。この方法は、各時刻のデータ変動を、少数の潜在クラスタ間の遷移で表現する。各時刻ステップは、複数あるうちの一つのクラスタへとアサインされる。そして、その時刻に観測された特徴量はそのクラスタの持つ観測パラメータに従って生成される。このモデルを利用すると、時系列データの伸縮を吸収しながら、各時刻ステップをクラスタリングすることができる。また、その潜在クラスタへの帰属パターンを識別することで、複雑な動作や発話の識別を高精度に行うことが可能である。
【０００５】
一方、Webなどの大規模データ解析では、得られる情報・特徴量が非常に多い。したがって、本来データの性質を規定する真の特徴量が、それ以外のノイズに対応するような情報に埋もれてしまうことがある。そこで、収集されたデータからデータのクラスタリングや識別に有用な特徴だけを選択する特徴選択手法が盛んに研究されている。最も良く利用されているのはサポートベクターマシン(ＳＶＭ)に代表される識別モデル学習の際にL1正則化項を入れる手法である。一方、生成モデルに対しては特徴ごとに有用であるか、あるいは有用でないかを表す潜在変数を明示的に割り当てるモデルがいくつか提案されている。それらの多くはデータのクラスタリングに利用されており、とくに、非特許文献１、２に記載の技術はクラスタ数を推定するとともにそれぞれのクラスタ固有の有用な特徴、不要な特徴を推定可能である。
【０００６】
また、高次元の観測量をもつ時系列データのクラスタリングとして、既存手法であるHMMに基づく手法の多くは、各時刻に観測される特徴に区別をつけないことが多い。したがって、各潜在クラスタは観測されるベクトルの単純な類似度に基づいていることになる。
【０００７】
また、HMMに基づく手法では、非常に多くの手法が事前に潜在クラスタの数を固定する必要がある。しかし、webで収集されるような巨大なデータに対して事前に適切なクラスタ数を設計することは難しい。近年、無限ＨＭＭ(ｉＨＭＭ)という手法が提案されている（非特許文献３、４）。この手法は与えられたデータから自動的にクラスタ数を推定することができる。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】K. Nowicki and T. A. B. Snijders. Estimation and prediction for stochastic blockstructures. Journal of the American Statistical Association, 96(455):1077−1087, 2001.
【非特許文献２】Guan, Y., Dy, J., and Jordan, M., “A Unified Probabilistic Model for Global and Local Unsupervised Feature Selection”, in Proceedings of the 27th International Conference on Machine Learning (ICML), 2011.
【非特許文献３】Teh, Y. W., Jordan, M. I., Beal, M. J., and Blei, D. M., “Hierarchical Dirichlet Process”, Journal of The American Statistical Association, Vol. 101, No. 476, po. 1566−1581, 2006.
【非特許文献４】Van Gael, J., Saatci, Y., Teh, Y. W., and Ghahramani, Z., “Beam Sampling for the Infinite Hidden Markov Model”, in Proceedings of the 25th International Conference on Machine Learning (ICML), 2008.
【発明の概要】
【発明が解決しようとする課題】
【０００９】
特に近年webなどで入手されるデジタルデータは非常に高次元かつスパースである。そのような状態では、各時刻で観測される特徴ベクトルのほとんどの要素はまったく同じ値をもち、ごく一部だけが変化するという状況に起こる。したがって、本当にクラスタに区別に有用な観測量のみを抜き出してクラスタリングをしなければ、多数の無益な観測量にクラスタ構造が隠されてしまう、という特徴選択の問題がある。上記の非特許文献３、４に記載の手法を用いても、上記の特徴選択の問題は解決されない。
【００１０】
次に特徴選択手法の問題であるが、特に時系列データに適用する際にはこれらの手法の仮定するデータの独立性が問題となる。これは、クラスタリングされる各観測量ベクトルは互いに独立かつ一様に分布する、という仮定である。多くの課題においてこの仮定は問題を起こさないが、時系列データはその例外の一つである。
【００１１】
時系列データは本質的に過去の情報の影響を受けるため、常に前の時刻との依存関係が発生する。したがって、この依存関係をモデル化しない手法は時系列データの特性を無視していることになる。
【００１２】
このように、既存の時系列データのクラスタリングモデルと特徴選択手法はいずれも一長一短であることが分かる。
【００１３】
本発明は、上記の事情を鑑みてなされたもので、時系列データが与えられたときに、クラスタ数を自動的に決定すると共に、時刻のクラスタリングに寄与する真の特徴量と、時刻のクラスリングに寄与しないノイズ特徴の影響を自動的に推定・分離して、時刻のクラスタリングを行うことができるクラスタリング装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１４】
上記の目的を達成するためにクラスタリング装置は、各時刻に観測された複数の特徴量からなる観測値の時系列データに基づいて、時刻のクラスタリングを実行するクラスタリング装置であって、各時刻ｔにおける各特徴量ｉからなる観測値ｘi,t、潜在クラスタ数Ｋ、時刻ｔが所属するクラスタを示すクラスタ割り当て変数ｚt、各クラスタｋに対して特徴量ｉが重要であるか否かを示す特徴選択変数ｒk,i、あるクラスタに対して特徴量ｉが重要である確率を示す特徴選択パラメータλi、特徴量ｉが時刻ｔの所属するクラスタｋに対して重要である場合に観測値ｘi,tを生成するために用いられる特徴観測パラメータθk,i、特徴量ｉが時刻ｔの所属するクラスタに対して重要でない場合に観測値ｘi,tを生成するために用いられるノイズ観測パラメータφi、時刻ｔ−１がクラスタｋに所属し、かつ、次の時刻ｔがクラスタｌに所属する確率を示すクラスタ遷移パラメータπk,ｌ、各クラスタｋに所属する時刻の割合を示す混合比パラメータβk、補助変数ｕt、ｐt,k及びハイパーパラメータγ、α、ａi、ｂi、ｃi、ｄi、ｅi、ｆiを記憶する記憶部と、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴選択パラメータλi、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπk,i、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ｃi、ｄiを取得し、各時刻の前記補助変数ｕtを決定すると共に前記潜在クラスタ数Ｋを決定し、各時刻ｔ及び各クラスタｋの前記補助変数ｐt,kを算出して、各時刻ｔの前記クラスタ割り当て変数ｚｔをサンプリングにより算出し、前記決定した前記補助変数ｕt及び前記潜在クラスタ数Ｋ、前記算出した前記補助変数ｐt,k、並びに前記算出した前記クラスタ割り当て変数ｚtによって、前記記憶部に記憶してある前記補助変数ｕt、前記潜在クラスタ数Ｋ、前記補助変数ｐt,k、及び前記クラスタ割り当て変数ｚtを更新し記憶するクラスタ割り当て変数更新手段と、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタｋ及び各特徴量ｉの前記特徴選択変数ｒk,iをサンプリングにより算出し、前記算出した前記特徴選択変数ｒk,iによって前記記憶部に記憶してある前記特徴選択変数ｒk,iを更新し記憶する特徴選択変数更新手段と、前記記憶部から前記潜在クラスタ数Ｋ、前記特徴選択変数ｒk,i、及び前記ハイパーパラメータei、fiを取得し、各特徴量ｉの前記特徴選択パラメータλiをサンプリングにより算出し、前記算出した前記特徴選択パラメータλiによって前記記憶部に記憶してある前記特徴選択パラメータλiを更新し記憶する特徴選択パラメータ更新手段と、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、及び前記ハイパーパラメータｃi、ｄiを取得し、各クラスタｋ及び各特徴量ｉの前記特徴観測パラメータθk,iをサンプリングにより算出し、前記算出した前記特徴観測パラメータθk,iによって前記記憶部に記憶してある前記特徴観測パラメータθk,iを更新し記憶する特徴観測パラメータ更新手段と、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、及び前記ハイパーパラメータａi、ｂiを取得し、各特徴量ｉの前記ノイズ観測パラメータφiをサンプリングにより算出し、前記算出した前記ノイズ観測パラメータφiによって前記記憶部に記憶してある前記ノイズ観測パラメータφiを更新し記憶するノイズ観測パラメータ更新手段と、前記記憶部から前記潜在クラスタ数Ｋ、前記混合比パラメータβk、及び前記ハイパーパラメータαを取得し、各クラスタｋ、ｌの前記クラスタ遷移パラメータπk,lをサンプリングにより算出し、前記算出した前記クラスタ遷移パラメータπk,lによって前記記憶部に記憶してある前記クラスタ遷移パラメータπk,lを更新し記憶するクラスタ遷移パラメータ更新手段と、前記記憶部から前記潜在クラスタ数Ｋ、及び前記ハイパーパラメータγを取得し、各クラスタｋの前記混合比パラメータβkをサンプリングにより算出し、前記算出した前記混合比パラメータβkによって前記記憶部に記憶してある前記混合比パラメータβkを更新し記憶する混合比パラメータ更新手段と、前記クラスタ割り当て変数更新手段、前記特徴選択変数更新手段、前記特徴選択パラメータ更新手段、前記特徴観測パラメータ更新手段、前記ノイズ観測パラメータ更新手段、前記クラスタ遷移パラメータ更新手段、及び前記混合比パラメータ更新手段による更新を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す終了条件判定手段と、を含んで構成されている。
【００１５】
本発明に係るクラスタリング方法は、各時刻に観測された複数の特徴量からなる観測値の時系列データに基づいて、時刻のクラスタリングを実行するクラスタリング装置におけるクラスタリング方法であって、前記クラスタリング装置は、各時刻ｔにおける各特徴量ｉからなる観測値ｘi,t、潜在クラスタ数Ｋ、時刻ｔが所属するクラスタを示すクラスタ割り当て変数ｚt、各クラスタｋに対して特徴量ｉが重要であるか否かを示す特徴選択変数ｒk,i、あるクラスタに対して特徴量ｉが重要である確率を示す特徴選択パラメータλi、特徴量ｉが時刻ｔの所属するクラスタｋに対して重要である場合に観測値ｘi,tを生成するために用いられる特徴観測パラメータθk,i、特徴量ｉが時刻ｔの所属するクラスタに対して重要でない場合に観測値ｘi,tを生成するために用いられるノイズ観測パラメータφi、時刻ｔ−１がクラスタｋに所属し、かつ、次の時刻ｔがクラスタｌに所属する確率を示すクラスタ遷移パラメータπk,ｌ、各クラスタｋに所属する時刻の割合を示す混合比パラメータβk、補助変数ｕt、ｐt,k及びハイパーパラメータγ、α、ａi、ｂi、ｃi、ｄi、ｅi、ｆiを記憶する記憶部と解析部とを備え、前記解析部は、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴選択パラメータλi、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、前記クラスタ遷移パラメータπk,i、前記混合比パラメータβk、前記ハイパーパラメータγ、α、ｃi、ｄiを取得し、各時刻の前記補助変数ｕtを決定すると共に前記潜在クラスタ数Ｋを決定し、各時刻ｔ及び各クラスタｋの前記補助変数ｐt,kを算出して、各時刻ｔの前記クラスタ割り当て変数ｚｔをサンプリングにより算出し、前記決定した前記補助変数ｕt及び前記潜在クラスタ数Ｋ、前記算出した前記補助変数ｐt,k、並びに前記算出した前記クラスタ割り当て変数ｚtによって、前記記憶部に記憶してある前記補助変数ｕt、前記潜在クラスタ数Ｋ、前記補助変数ｐt,k、及び前記クラスタ割り当て変数ｚtを更新し記憶するクラスタ割り当て変数更新ステップと、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、前記特徴観測パラメータθk,i、前記ノイズ観測パラメータφi、及び前記ハイパーパラメータei、fiを取得し、各クラスタｋ及び各特徴量ｉの前記特徴選択変数ｒk,iをサンプリングにより算出し、前記算出した前記特徴選択変数ｒk,iによって前記記憶部に記憶してある前記特徴選択変数ｒk,iを更新し記憶する特徴選択変数更新ステップと、前記記憶部から前記潜在クラスタ数Ｋ、前記特徴選択変数ｒk,i、及び前記ハイパーパラメータei、fiを取得し、各特徴量ｉの前記特徴選択パラメータλiをサンプリングにより算出し、前記算出した前記特徴選択パラメータλiによって前記記憶部に記憶してある前記特徴選択パラメータλiを更新し記憶する特徴選択パラメータ更新ステップと、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、及び前記ハイパーパラメータｃi、ｄiを取得し、各クラスタｋ及び各特徴量ｉの前記特徴観測パラメータθk,iをサンプリングにより算出し、前記算出した前記特���観測パラメータθk,iによって前記記憶部に記憶してある前記特徴観測パラメータθk,iを更新し記憶する特徴観測パラメータ更新ステップと、前記記憶部から前記観測値ｘi,t、前記潜在クラスタ数Ｋ、前記クラスタ割り当て変数ｚｔ、前記特徴選択変数ｒk,i、及び前記ハイパーパラメータａi、ｂiを取得し、各特徴量ｉの前記ノイズ観測パラメータφiをサンプリングにより算出し、前記算出した前記ノイズ観測パラメータφiによって前記記憶部に記憶してある前記ノイズ観測パラメータφiを更新し記憶するノイズ観測パラメータ更新ステップと、前記記憶部から前記潜在クラスタ数Ｋ、前記混合比パラメータβk、及び前記ハイパーパラメータαを取得し、各クラスタｋ、ｌの前記クラスタ遷移パラメータπk,lをサンプリングにより算出し、前記算出した前記クラスタ遷移パラメータπk,lによって前記記憶部に記憶してある前記クラスタ遷移パラメータπk,lを更新し記憶するクラスタ遷移パラメータ更新ステップと、前記記憶部から前記潜在クラスタ数Ｋ、及び前記ハイパーパラメータγを取得し、各クラスタｋの前記混合比パラメータβkをサンプリングにより算出し、前記算出した前記混合比パラメータβkによって前記記憶部に記憶してある前記混合比パラメータβkを更新し記憶する混合比パラメータ更新ステップと、前記クラスタ割り当て変数更新ステップ、前記特徴選択変数更新ステップ、前記特徴選択パラメータ更新ステップ、前記特徴観測パラメータ更新ステップ、前記ノイズ観測パラメータ更新ステップ、前記クラスタ遷移パラメータ更新ステップ、及び前記混合比パラメータ更新ステップによる更新を任意の順番で実行する過程を、所定の終了条件を満足するまで繰り返す終了条件判定ステップと、を実行することを特徴とする。
【００１６】
本発明によれば、クラスタ割り当て変数を更新するときにクラスタ数を決定すると共に、クラスタ割り当て変数の更新、特徴選択変数の更新、特徴選択パラメータの更新、特徴観測パラメータの更新、ノイズ観測パラメータの更新、クラスタ遷移パラメータの更新、及び混合比パラメータの更新を繰り返すことにより、時系列データが与えられたときに、クラスタ数を自動的に決定すると共に、時刻のクラスタリングに寄与する真の特徴量と、時刻のクラスリングに寄与しないノイズ特徴の影響を自動的に推定・分離して、時刻のクラスタリングを行うことができる。
【００１７】
本発明に係るプログラムは、上記のクラスタリング装置の各手段としてコンピュータを機能させるためのプログラムである。
【発明の効果】
【００１８】
以上説明したように、本発明のクラスタリング装置、方法、及びプログラムによれば、クラスタ割り当て変数を更新するときにクラスタ数を決定すると共に、クラスタ割り当て変数の更新、特徴選択変数の更新、特徴選択パラメータの更新、特徴観測パラメータの更新、ノイズ観測パラメータの更新、クラスタ遷移パラメータの更新、及び混合比パラメータの更新を繰り返すことにより、時系列データが与えられたときに、クラスタ数を自動的に決定すると共に、時刻のクラスタリングに寄与する真の特徴量と、時刻のクラスリングに寄与しないノイズ特徴の影響を自動的に推定・分離して、時刻のクラスタリングを行うことができる、という効果が得られる。
【図面の簡単な説明】
【００１９】
【図１】本発明のコンセプトを表わした模式図である。
【図２】本発明の実施の形態に係るクラスタリング装置の構成を示す概略図である。
【図３】本発明の実施の形態に係るクラスタリング装置におけるクラスタリング処理ルーチンの内容を示すフローチャートである。
【図４】本発明の実施の形態に係るクラスタリング装置のクラスタ割り当て変数更新部のフローチャートである。
【図５】クラスタ割り当て変数更新部による補助変数ｕのサンプリング法を示すフローチャートである。
【図６】クラスタ割り当て変数更新部による未使用クラスタの整理法を示すフローチャートである。
【図７】本発明の実施の形態に係るクラスタリング装置の特徴選択変数更新部のフローチャートである。
【図８】本発明の実施の形態に係るクラスタリング装置の特徴選択パラメータ更新部のフローチャートである。
【図９】本発明の実施の形態に係るクラスタリング装置の特徴観測パラメータ更新部のフローチャートである。
【図１０】本発明の実施の形態に係るクラスタリング装置のノイズ観測パラメータ更新部のフローチャートである。
【図１１】本発明の実施の形態に係るクラスタリング装置のクラスタ遷移パラメータ更新部のフローチャートである。
【図１２】本発明の実施の形態に係るクラスタリング装置の混合比パラメータ更新部のフローチャートである。
【図１３】本発明の実施の形態に係るクラスタリング装置のハイパーパラメータ更新部のフローチャートである。
【図１４】文書データの解析結果を示す図である。
【発明を実施するための形態】
【００２０】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００２１】
＜概要＞
まず、本発明の概要について説明する。
【００２２】
これまで、時系列データに対して各時刻を少数の潜在的なクラスタに帰属することでモデル化する手法が多数提案されている。この枠組みでは、一般に多次元ベクトルとして表現可能な観測量が各時刻で観測される。
【００２３】
そして、この特徴量の分布の類似度でもって各時刻を少数のパターン（クラスタ）に分類する。この隠れたクラスタの推定によって、例えば人間の動作や音声といった複雑な時系列データを精度良く解析・認識することが可能である。
【００２４】
本発明では、これらの手法のもつ欠点のうち、大きく２つの部分を改善した手法を提案する。まず一点目として、観測量のうち、時刻のクラスタリングに寄与しないノイズ特徴量の存在を仮定し、その影響を自動的に排除する。二点目として、クラスタの数を事前に決定することなく、与えられたデータを表現する上で最も適切なクラスタ数を自動的に決定する。
【００２５】
図１は提案するクラスタリングシステムのコンセプトを模式的に表した図である。
【００２６】
クラスタリングしたい観測値の時系列データをＸとする。観測値の時系列データはN種類の特徴量、すなわちN次元の特徴ベクトルが総時間ステップ数Tに渡って観測されたものとする。特徴量をi=1,2,...,N、時間ステップをt=1，2,...,Tでインデキシングする。時刻tにおけるi番目の特徴量の観測値をxt,iとする。すなわち、Ｘ＝{ｘｉ，ｔ},ｉ＝１,...,N,ｔ＝１,...,Nである。ｘｉ,ｔは0,１のバイナリデータ、あるいは自然数、実数、あるいはシンボルデータである。
【００２７】
この観測データと、ハイパーパラメータと呼ばれる定数とを入力することで、クラスタリングシステムから解析結果が出力される。主要な解析出力はクラスタ数K、クラスタ割り当て変数Z、特徴選択変数Rである。クラスタ数KはXに内在する特徴観測パターン、すなわち潜在クラスタ数である。Z={ｚｔ}，t=1,...,Tは、全時刻ステップに対する潜在クラスタの割り当てを表す。zt=k，k=1,...,Kは、時刻ステップkの特徴量がクラスタkから確率的に生成されたことを表す。また、文脈に応じてztは1-of-Kベクトルとしても利用する。すなわち、zt=kのとき、zt,k=1，ｚt,l≠k=0とする。特徴選択変数R={rk,i}，k=1,...,K，i=1,...,Nは全クラスタにおける特徴選択結果を表す。rk,i∈{0,1}は、i番目の特徴量がk番目のクラスタの特性を表す真の特徴(rk,i=1)か、そうでないノイズ特徴(rk,i=0)かを表す。図1ではrzt,i=1となる部分を濃く表示している。ほとんどの白い部分は重要でないノイズ特徴と判断されている。k=6のクラスタは、特徴的なものが存在しない、すなわちノイズのみが観測された期間であると解釈される。
【００２８】
＜システム構成＞
図２は、本発明の実施の形態に係るクラスタリング装置１００の模式図の一例である。クラスタリング装置１００は、解析したい観測データの入力を受け付ける入力部１と、クラスタリングを行う解析部２と、入力および計算された変数を記憶する記憶部３と、クラスタリング結果を出力する出力部４と、を備えている。
【００２９】
入力部１は、既知のビデオカメラ、マイクロフォン、記憶装置などの入力器により実現される。また、出力部４はディスプレイ、プリンタ、磁気ディスクなどで実装される。
【００３０】
解析部２は、初期設定部２１、変数計算部２２、及び解析終了条件判定部２１１を備えている。変数計算部２２は、クラスタ割り当て変数更新部２３、特徴選択変数更新部２４、特徴選択パラメータ更新部２５、特徴観測パラメータ更新部２６、ノイズ観測パラメータ更新部２７、クラスタ遷移パラメータ更新部２８、混合比パラメータ更新部２９、及びハイパーパラメータ更新部２１０を備えている。
【００３１】
記憶部３は、定数３１、観測値３２、潜在クラスタ数３３、クラスタ割り当て変数３４、特徴選択変数３５、特徴選択パラメータ３６、特徴観測パラメータ３７、ノイズ観測パラメータ３８、クラスタ遷移パラメータ３９、混合比パラメータ３１０、補助変数３１１、及びハイパーパラメータ３１２から構成される。
【００３２】
初期設定部２１は、入力部1より入力された情報を記憶部３の各部へと記憶するとともに、記憶部３の各部の値の初期化を行う。
【００３３】
変数計算部２２は、記憶部３に保存された情報を利用し、記憶部３の各部の情報を更新、保存する。
【００３４】
クラスタ割り当て変数更新部２３は記憶部３の情報を用いて潜在クラスタ数３３、クラスタ割り当て変数３４、特徴選択変数３５、特徴観測パラメータ３７、クラスタ遷移パラメータ３９、混合比パラメータ３１０、補助変数３１１を更新し、保存する。
【００３５】
特徴選択変数更新部24は、記憶部３の情報を用いて特徴選択変数３５を更新し、保存する。
【００３６】
特徴選択パラメータ更新部２５は、記憶部３の情報を用いて特徴選択パラメータ３６を更新し、保存する。特徴観測パラメータ更新部２６は記憶部３の情報を用いて特徴観測パラメータ３７を更新し、保存する。ノイズ観測パラメータ更新部２７は記憶部３の情報を用いてノイズ観測パラメータ３８を更新し、保存する。クラスタ遷移パラメータ更新部２８は、記憶部３の情報を用いてクラスタ遷移パラメータ３９を更新し、保存する。混合比パラメータ更新部２９は、記憶部３の情報を用いて混合比パラメータ３１０を更新し、保存する。ハイパーパラメータ更新部２１０は記憶部３の情報を用いてハイパーパラメータ３１２を更新し、保存する。
【００３７】
変数計算部２２の処理の流れは、クラスタ割り当て変数更新部２３、特徴選択変数更新部２４、特徴選択パラメータ更新部２５、特徴観測パラメータ更新部２６、ノイズ観測パラメータ更新部２７、クラスタ遷移パラメータ更新部２８、混合比パラメータ更新部２９、ハイパーパラメータ更新部２１０の順番に限らず、任意の順番で実施してよい。また、必要がなければハイパーパラメータ更新部２１０による更新は実施しなくてもよい。
【００３８】
解析終了条件判定部２１１は、変数計算部２２および記憶部３を監視し、解析終了条件を満たしたかどうかを計算する。
【００３９】
出力部４は解析終了条件判定部２１１を監視し、解析終了条件判定部２１１により解析終了条件が満たされたと判断された場合、クラスタリング結果の出力を実施する。
【００４０】
＜クラスタリング装置の動作＞
次に、本実施の形態に係るクラスタリング装置１００の作用について説明する。まず、各種の定数及び観測値の時系列データが入力部１を介してクラスタリング装置１００に入力されると、クラスタリング装置１００において、図３に示すクラスタリング処理ルーチンが実行される。
【００４１】
まず、ステップＳ１０１において、初期設定部２１によって、入力部１により入力された、定数及び観測値の時系列データを受け付けて、記憶部３に格納する。
【００４２】
そして、ステップＳ１０２において、初期設定部２１によって、記憶部３の各種パラメータを初期化する。ステップＳ１０３では、繰り返し回数を示す変数ｉｔｒに初期値１を設定する。
【００４３】
次のステップＳ１０４では、クラスタ割り当て変数更新部２３によって、記憶部３の情報に基づいて、クラスタ割り当て変数を更新して、記憶部３に記憶する。ステップＳ１０５では、特徴選択変数更新部２４によって、記憶部３の情報に基づいて、特徴選択変数を更新して、記憶部３に記憶する。ステップＳ１０６では、特徴選択パラメータ更新部２５によって、記憶部３の情報に基づいて、特徴選択パラメータを更新して、記憶部３に記憶する。
【００４４】
そして、ステップＳ１０７において、特徴観測パラメータ更新部２６によって、記憶部３の情報に基づいて、特徴観測パラメータを更新して、記憶部３に記憶する。ステップＳ１０８では、ノイズ観測パラメータ更新部２７によって、記憶部３の情報に基づいて、ノイズ観測パラメータを更新して、記憶部３に記憶する。ステップＳ１０９では、クラスタ遷移パラメータ更新部２８によって、記憶部３の情報に基づいて、クラスタ遷移パラメータを更新して、記憶部３に記憶する。
【００４５】
次のステップＳ１１０では、混合比パラメータ更新部２９によって、記憶部３の情報に基づいて、混合比パラメータを更新して、記憶部３に記憶する。ステップＳ１１１において、ハイパーパラメータ更新部２１０によって、記憶部３の情報に基づいて、ハイパーパラメータを更新して、記憶部３に記憶する。
【００４６】
そして、ステップＳ１１２において、解析終了条件判定部２１１によって、解析終了条件を満足したか否かを判定する。解析終了条件は、例えば、繰り返し変数ｉｔｒが、予め決めておいた繰り返し回数になったことである。解析終了条件を満足していない場合には、ステップＳ１１３において、変数ｉｔｒを１インクリメントして、上記ステップＳ１０４へ戻る。一方、終了条件を満足する場合には、ステップＳ１１４において、上記ステップＳ１０４で最終的に更新されたクラスタ割り当て変数に基づくクラスタリング結果を出力部４により出力して、クラスタリング処理ルーチンを終了する。
【００４７】
次に、図３の各ステップの処理を実行する、クラスタリング装置１００の構成および処理の詳細について説明する。
【００４８】
＜初期設定部２１＞
初期設定部２１は、入力部１を通じて、解析したい観測データの時系列データＸおよび複数の定数の入力を受け付ける。
【００４９】
入力データは以下のような形式であり、初期設定部２１でロードされる。観測値の時系列データは、Ｎ種類の特徴量、すなわちＮ次元の特徴ベクトルが総時間ステップ数Ｔに渡って観測されたものとする。特徴量をi=1,2,...,N、時間ステップをt=1,2,...,Tでインデキシングする。時刻tにおける観測値のi番目の特徴量をxi,tとする。
【００５０】
すなわち、X={xi,t}、i=1,...,N、t=1,...,Nである。xi,tは0,1のバイナリデータ、あるいは自然数、実数、あるいはシンボルデータである。
【００５１】
また、観測値の時系列データとは別に、以下の定数も別途入力が必要である。Tは観測データの最大時間ステップサイズである。Nは観測データの特徴量のサイズ、すなわち特徴ベクトルの次元数である。K0は潜在クラスタ数３３の初期値である。Lは解析終了条件判定部２１１に使用する終了条件である。これは例えば変数計算部２２の最大繰り返し回数、あるいは何らかの評価値(適切に定められた評価関数)の変化率の閾値など、ユーザの目的に応じて定める。
【００５２】
さらに、モデルに必要なハイパーパラメータと呼ばれる定数も同時に入力する。γは混合比パラメータ３１０の計算で必要となる。αはクラスタ遷移パラメータ３０９の計算で必要となる。a1,...,ai,...,aN，b1,..., bi,...,bNは主にノイズ観測パラメータ３８の計算で必要となる。c1,...,ci,...,cN，d1,...,di,...,dNは主に特徴観測パラメータ３７の計算で必要となる。e1,...,ei,...,eN，f 1,..., fi,...,fNは主に特徴選択パラメータ３６の計算で必要となる。
【００５３】
各ハイパーパラメータの定義域はxi,tの形式およびそれに応じたモデル化に従う。例えば最も単純なバイナリデータの例、すなわちxi,t∈{0,1}の場合、γは0より大きい実数、ai,bi,ci,di,ei,fiは1以上の実数が望ましい。
【００５４】
これらの入力を受け付けたのち、初期設定部２１は、記憶部３の各部の値を初期化する。定数３１には入力されたT,N,K0,Lを設定する。観測値３２には入力されたXを設定する。潜在クラスタ数３３にはK0を設定する。
【００５５】
潜在クラスタ割り当て変数３４には、1,...,Kの値を各ztにランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。
【００５６】
特徴選択変数３５には、設定した潜在クラスタ数Kに従って各rk,iに対し0または1の値をランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。
【００５７】
特徴選択パラメータ３６には、各λiに対し0以上1未満の値をランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。
【００５８】
特徴観測パラメータ３７には、xi,tが0または1のバイナリデータの場合、設定されたKに従って各θk,iに対し0以上1未満の値をランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。xi,tが別の形式の場合、適切に選択されたモデルに合わせて値を設定する。
【００５９】
ノイズ観測パラメータ３８には、xi,tが0または1のバイナリデータの場合、各φiに対し0以上1未満の値をランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。xi,tが別の形式の場合、適切に選択されたモデルに合わせて値を設定する。
【００６０】
クラスタ遷移パラメータ３９には、設定されたKに従って各πk,lに対し0以上1未満の値をランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。ただし、全てのk∈{1,...,K＋1}に対し、ΣK+1l=1πk,l=1が満たされる必要がある。
【００６１】
混合比パラメータ３１０には、設定されたKに従って各βｋに対し0以上1未満の値をランダムに設定するか、あるいは後述するモデルに基づき、確率的に値を設定する。ただし、ΣK+1k=1πｋ，ｌk=1が満たされる必要がある。
【００６２】
補助変数３１１は、潜在クラスタ数３３、クラスタ割り当て変数３４、クラスタ遷移パラメータ３９の設定後に、後述するモデルに従って設定する。
【００６３】
ハイパーパラメータ３１２には、入力されたハイパーパラメータを設定する。
【００６４】
なお、xi，tの形式に従って、変数計算部２２の具体的な実装および記憶部３の各部のメモリサイズは変更されうる。本実施の形態では、最も単純なバイナリデータの例、すなわちxi，t∈{0，1}の場合の例について説明するが、他のデータ形式の場合も一般的な教科書、例えば参考文献（ビショップ, C. M., “パターン認識と機械学習”, シュプリンガージャパン, 2007.）を参照することで簡単に実装することが可能である。
【００６５】
初期化処理の終了後、処理を変数計算部２２に移す。
【００６６】
＜モデルの定式化（記憶部３）＞
それでは主に記憶部３に保存される情報について、モデルにおけるそれらの意味とともに説明する。なお、ここではxi，tがバイナリデータ、すなわちxi，t∈{0,1}の場合を例に説明する。
【００６７】
本発明では、観測された時系列データXを少数の観測値のパターン、すなわち潜在クラスタに分解して、データの大まかな生成パターンをとらえやすくする。また、それによってノイズに対してロバストな予測モデルの構築を実現する。
【００６８】
まず、潜在クラスタ数３３はデータに内在するパターン、つまりクラスタの数を表す。クラスタ数Ｋは一般には未知であるが、本モデルでは真のクラスタ数がどのような値でも適切に推定・解析可能なディリクレ過程混合モデル(Dirichlet Process Mixture, DPM)、より具体的にはそのHMMモデル応用であるiHMM(infinite HMM（上記非特許文献３参照）) を利用する。
【００６９】
まず、全Ｔ個の時刻ステップのうち、クラスタｋによって説明されている時刻ステップの割合を表す混合比パラメータ３１０を、以下の（１）式に示す確率モデルでサンプリングする。
【００７０】
【数１】
0005851205.tif 000002
【００７１】
ここで、β＝（β１，．．．，βK，．．．）Tは無限次元のベクトルである。Stick(・)はStick breaking processと呼ばれる確率過程である。より具体的には各βkは以下の（２）式に示す確率を用いて、以下の（３）式に従って決定される。
【００７２】
【数２】
0005851205.tif 000003
【００７３】
ここでBetaはベータ分布を表す。上記（３）式によれば、加算無限個のクラスタ数を考えたときにベクトルβは全ての要素が0以上かつ、足して1となる。
【００７４】
さて、実際の解析時にはKは有限の値をとる。すなわち、データは有限なので必ず有限Ｋ個のクラスタで表現できるはずである。実際の実装上ではβはK+1次元ベクトルとなる。βｋ≦Kは、全Ｔ時刻のうち、クラスタkで説明される時刻がどれだけあるかを表す割合である。βK+1はunrepresented clusters、すなわち理論上は存在するが現在はデータのモデル化に不要なクラスタ集合のもつ重みである。βK+１は次の（４）式で決定する。
【００７５】
【数３】
0005851205.tif 000004
【００７６】
Stick breaking process の定義より、どのようなＫに対してもベクトルβは全ての要素が0以上かつ、足して1となる。混合比パラメータ３１０にはこのＫ＋１次元のベクトルが格納される。
【００７７】
次に、このモデルでは時系列データの時間発展をクラスタ割り当ての時間遷移で表現する。すなわち、複雑な時間変動をおこなう時系列データを、たかだかＫ種類のパターン間の遷移で表現する。この遷移は、クラスタ遷移パラメータ３９によって制御される。クラスタ遷移パラメータπｋ、ｌは時刻ｔ−１におけるクラスタがkだった場合、時刻t におけるクラスタ割り当てがlになる確率をあらわす。クラスタkからの遷移先の確率をπｋ＝（πｋ，１，πｋ，２，．．．，πｋ，ｌ，．．．，）Ｔで表現する。このベクトルも理論上は無限次元なので、以下の（５）式に示すように、クラスタ遷移パラメータの定義はディリクレ過程(Dirichlet Process, DP) を必要とする。DPは有限次元におけるDirichlet分布を拡張したもので, Stick-breaking processと同様無限次元のベクトルを生成する。
【００７８】
【数４】
0005851205.tif 000005
【００７９】
実際にはβはＫ＋１次元で有限なので、実装上は以下の（６）式で計算を行う。
【００８０】
【数５】
0005851205.tif 000006
【００８１】
Dirichlet(・)は有限次元のディリクレ分布である。得られるπkはK+1次元のベクトルであり、これをクラスタ遷移パラメータ３９として格納する。
【００８２】
このパラメータを利用して、実際のクラスタ遷移を実現する。この結果得られるクラスタ割り当ての時系列はクラスタ割り当て変数３４に記録される。クラスタ割り当て変数ｚｔ＝ｋは時刻t におけるクラスタ割り当てである。ｚｔは以下の（７）式で決定する。
【００８３】
【数６】
0005851205.tif 000007
【００８４】
ここでMultinomial(・) は多項分布を表す。なお、z1の決定については種々の方法があるが、ここでは便宜的にz0=1を仮定して上記（７）式を適用することとする。
【００８５】
これで、与えられた時系列データXを、時間軸方向にK種のクラスタへと分解することができる。次に、各クラスタについてＮ種の観測量のうち真に役に立つクラスタの特徴を表す観測量とノイズでしかない観測量とを分離するモデルを説明する。
【００８６】
まず、ｉ番目の観測値があるクラスタにとって重要な観測量、すなわち”特徴”である確率をλｉとする。この値は以下の（８）式で決定する。
【００８７】
【数７】
0005851205.tif 000008
【００８８】
上記（８）式に従ってサンプリングされるλｉの値を特徴選択パラメータ３６に格納する。
【００８９】
次に、各クラスタkにおいて、i番目の観測値を重要な”特徴”とみなすか、あるいは単なる”ノイズ”観測量とみなすかを表す特徴選択変数３５を、以下の（９）式で決定する。
【００９０】
【数８】
0005851205.tif 000009
【００９１】
ここでBernoulli(・)はベルヌーイ分布を表す。ｒｋ，ｉ＝１のとき、ｉ番目の観測値はクラスタｋにとって重要な” 特徴”であり、ｒｋ，ｉ＝０のときはただのノイズ観測量と考える。このｒｋ，ｉを全てのk,iについて確率的にサンプリングし、特徴選択変数３５に格納する。これによって、時刻tの観測値のi番目の特徴量ｘｉ，ｔを、その時刻ｔが所属するクラスタ番号k=ztと、特徴として選択されたかどうかを表すrk，iで分類することができる。
【００９２】
最後に、観測値の確率モデルを定義する。観測値ｘｉ，ｔがモデルから実際に生成される確率の計算には、zt，rk，iで指定される観測パラメータが必要である。zt=kかつrk，i=1の場合には、特徴観測パラメータ３７を使ってデータを生成する。一方zt=kかつrk，i=0の場合には、ノイズ観測パラメータ３８を使ってデータを生成する。
【００９３】
観測値がバイナリデータの場合、それぞれのパラメータは、以下の（１０）式、（１１）式のように決定される。
【００９４】
特徴観測パラメータはクラスタkごとに特有の値を持つと考えられるので、（１０）式で定式化される。
【００９５】
【数９】
0005851205.tif 000010
【００９６】
一方、ノイズ観測パラメータはクラスタの番号によらないと考え、（１１）式で定式化される。
【００９７】
【数１０】
0005851205.tif 000011
【００９８】
特徴観測パラメータ及びノイズ観測パラメータそれぞれのサンプリング結果を、記憶部３の対応する部分へと格納する。
【００９９】
最終的な観測値xi，tの生成確率は、以下の（１２）式のように評価される。
【０１００】
【数１１】
0005851205.tif 000012
【０１０１】
以上の観測値を除く全ての変数を何らかの評価尺度の意味で最適な値に決定することが本クラスタリングシステムの主な目的である。結果として、観測量xi，tだけを観察することで、時系列データの生成パターンを、パターンを区別する上で重要な特徴だけを利用して理解、あるいは上記図１のように可視化して出力することが可能となる。
【０１０２】
＜モデルのバリエーションについて＞
以上のモデルは観測値がバイナリデータの場合である。しかし、システムの構成をほとんど変えることなく、他の場合にも適用可能なシステムを実現できる。
【０１０３】
まず、観測値が非負整数の場合には、観測確率式である上記（１２）式の関数として例えばポアソン分布が利用できる。その場合、関数の共役性から、上記（１０）式、上記（１１）式はそれぞれガンマ分布がふさわしい。それに応じてハイパーパラメータの個数や定義域を適切に変更すれば非負整数観測値に対応できる。
【０１０４】
観測値が実数の場合には、観測確率式である上記（１２）式の関数として例えば正規分布が利用できる。その場合、関数の共役性から、上記（１０）式、（１１）式はそれぞれガウス-ガンマ分布あるいはガウス-逆ガンマ分布がふさわしい。それに応じてハイパーパラメータの個数や定義域を適切に変更すれば実数観測値に対応できる。
【０１０５】
観測値がシンボルの場合には、観測確率式である上記（１２）式の関数として例えば多項分布が利用できる。その場合、関数の共役性から上記（１０）式、（１１）式はそれぞれディリクレ分布がふさわしい。それに応じてハイパーパラメータの個数や定義域を適切に変更すればシンボル観測値に対応できる。以上の各関数やパラメータの設定については標準的な教科書、例えば、上記参考文献などに記載されている。
【０１０６】
＜変数の最適化と出力(変数計算部２２)＞
変数計算部２２の目的は、観測データＸが与えられた際に、記憶部３に格納された変数やパラメータを何らかの意味で最適な値に決定し、解析結果として出力することである。
【０１０７】
本発明では、上述のバイナリ観測値による定式化の場合を例にとって、効率の良い解析が可能とされているスライスサンプリング法（上記非特許文献４を参照）を利用した最適化法の実装方法を説明する。
【０１０８】
なお、各変数やパラメータの最適化には、モデルのバリエーションによらず他にも最尤推定法や事後確率最大化法、ギブスサンプリング法、変分ベイズ法、マルコフ連鎖モンテカルロ法、あるいは種々の数値最適化法などの様々な標準的な手法（上記参考文献を参照）を用いることも可能である。
【０１０９】
変数計算部２２は、記憶部３に定義されている推定可能な変数やパラメータに対応する形で、各変数とパラメータの更新部分に分離されている。変数・パラメータの更新の順番は任意であるが、通常は以下に説明する順番に一度ずつ更新を行い、一回の変数更新処理を完了する。
【０１１０】
＜クラスタ割り当て変数更新部２３＞
クラスタ割り当て変数更新部２３では、クラスタ割り当て変数３４を更新する。また、その過程で補助変数３１１も更新する。また、クラスタ割り当て変数Zの更新では潜在クラスタ数３３も更新される。それに従って、Kに依存する変数やパラメータである特徴選択変数３５、特徴観測パラメータ３７、クラスタ遷移パラメータ３９、混合比パラメータ３１０も更新される可能性がある。クラスタ割り当て変数更新部２３のフローチャートを図４に示す。
【０１１１】
まず、ステップＳ１２１において、記憶部３から、T, N, K, x, z, r, θ, φ, π, λ, β, γ, α, c, dを読み込む。ステップＳ１２２において、所定の条件を満たすまで、補助変数ｕｔ，ｔ＝１，．．．，Ｔをサンプリングする。
【０１１２】
上記ステップＳ１２２は、図５に示すフローチャートにより実現される。
【０１１３】
ステップＳ１５１において、時刻ｔに初期値１を設定し、ステップＳ１５２で、以下の（１３）式に示す確率でサンプリングして、補助変数ｕｔを決定する。
【０１１４】
【数１２】
0005851205.tif 000013
【０１１５】
ここでUniform(・) は一様分布を表す。
【０１１６】
次のステップＳ１５３では、時刻ｔを１インクリメントとして、ステップＳ１５４において、時刻ｔがＴより大きいか否かを判定する。時刻ｔがＴ以下である場合には、上記ステップＳ１５２へ戻る。一方、時刻ｔがＴより大きい場合には、ステップＳ１５５へ進む。
【０１１７】
ステップＳ１５５では、時刻ｔに初期値１を設定し、ステップＳ１５６において、補助変数utについて、以下の条件を満たすか否かを判定する。
【０１１８】
uｔ≧πｚ４−１，Ｋ＋１
【０１１９】
上記の条件を満たす場合には、ステップＳ１５７において、時刻ｔを１インクリメントして、ステップＳ１５８において、時刻ｔがＴより大きいか否かを判定する。時刻ｔがＴ以下である場合には、上記ステップＳ１５６へ戻る。一方、時刻ｔがＴより大きい場合には、全てのtについてuｔ≧πｚ４−１，Ｋ＋１であることを確認できたとして、上記ステップＳ１２２の処理を終了する。
【０１２０】
一方、上記の条件を満たさない場合には、次に示す手順でKを一つ増加させ、再度uをサンプリングする。
【０１２１】
ステップＳ１５９では、ｕのサンプリング結果を破棄し、ステップＳ１６０において、以下の式に従って、変数ｂをサンプリングする。
【０１２２】
b〜Beta(1，γ)
【０１２３】
なお、b〜Beta(1，γ)で、ベータ分布Beta(1，γ)より独立かつ一様にbをサンプリングすることを表す。
【０１２４】
そして、ステップＳ１６１において、以下の式に従って、β＊、β＊＊を計算して、βを拡張し更新する。
【０１２５】
β＊←βK+1＊b，
β＊＊←βK+1−β＊
β←（β１，β２，．．．，βＫ，β＊，β＊＊）
【０１２６】
次のステップＳ１６２では、変数ｉに初期値１を設定し、ステップＳ１６３において、上記（１０）式に従ってθＫ＋１，ｉをサンプリングして更新する。ステップＳ１６４では、変数ｉを１インクリメントして、ステップＳ１６５において、変数ｉがＮより大きいか否かを判定する。変数ｉがＮ以下である場合には、上記ステップＳ１６３へ戻る。一方、変数ｉがＮより大きい場合には、全てのｉについてθＫ＋１，ｉをサンプリングしたとして、ステップＳ１６６へ進む。
【０１２７】
次のステップＳ１６６では、変数ｉに初期値１を設定し、ステップＳ１６７において、上記（９）式に従ってｒＫ＋１，ｉをサンプリングして更新する。ステップＳ１６８では、変数ｉを１インクリメントして、ステップＳ１６９において、変数ｉがＮより大きいか否かを判定する。変数ｉがＮ以下である場合には、上記ステップＳ１６７へ戻る。一方、変数ｉがＮより大きい場合には、全てのｉについてｒＫ＋１，ｉをサンプリングしたとして、ステップＳ１７０へ進む。
【０１２８】
ステップＳ１７０では、潜在クラスタ数Ｋを、Ｋ＋１に更新する。そして、ステップＳ１７１において、クラスタ遷移パラメータ更新部２８を用いてクラスタ遷移パラメータπk，lを更新して、上記ステップＳ１５１へ戻る。
【０１２９】
上記のように、utのサンプリングの後、全てのtについてut≧πzｔ−１，K+1であることを確認した時点で補助変数utが決定する。
【０１３０】
この補助変数utは、理論上無限個存在する潜在クラスタを「足切り」して、有限のクラスタに限定する効果を持つ。より正確には、πｋによるクラスタ遷移を可能性の高い遷移だけに限定している。
【０１３１】
また、上記図４のフローチャートのステップＳ１２３では、時刻ｔに初期値１を設定し、ステップＳ１２４で、変数ｋを初期値１に設定する。
【０１３２】
そして、ステップＳ１２５において、以下の（１４）式に示す補助変数pt，kを計算する。
【０１３３】
【数１３】
0005851205.tif 000014
【０１３４】
ただし、ｘ1:N,1:tは、｛ｘi,t'｝i=1,...,N,ｔ'=1,...,tを要素とするベクトルを表わす。この補助変数pt，kは、時刻tまでの観測値の情報を得たときに、時刻tがクラスタk に割り振られる確率を表す。
【０１３５】
この補助変数pt，kは、以下の（１５）式に示す漸化式に従って計算される。
【０１３６】
【数１４】
0005851205.tif 000015
【０１３７】
ただし、∝は比例を表す。上記（１５）式の右辺第一項は時刻tにおける全観測量の尤度をzt=kと仮定したときに計算した量である。すなわち、上記（１２）式をt，zt=kに固定してi=1，．．．，Ｎまで掛け込んだ値である。これによって、各特徴量ｉが時刻ｔの所属するクラスタｋに対して重要であるか否かに基づいて、補助変数ｐt,kが算出される。
【０１３８】
ステップＳ１２６では、ｋを１インクリメントして、ステップＳ１２７において、変数ｋがＫより大きいか否かを判定する。変数ｋがＫ以下である場合には、上記ステップＳ１２５へ戻る。一方、変数ｋがＫより大きい場合には、全てのｋについてｐt,kを計算したとして、ステップＳ１２８へ進む。
【０１３９】
ステップＳ１２８では、上記ステップＳ１２５で計算されたpt，kをΣkpt，k=1となるように正規化する。ステップＳ１２９では、時刻ｔを１インクリメントして、ステップＳ１３０において、時刻ｔがＴより大きいか否かを判定する。時刻ｔがＴ以下である場合には、上記ステップＳ１２４へ戻る。一方、時刻ｔがＴより大きい場合には、この補助変数をｔ＝１，２，．．．，Ｔの順番で、全てのｋ＝１，２，．．．，Ｋに対して計算したとして、ステップＳ１３１へ進む。
【０１４０】
ステップステップＳ１３１において、時刻ｔをＴに設定する。そして、ステップＳ１３２において、変数ｋを初期値１に設定する。ステップＳ１３３では、ｚt＝ｋとなる確率を以下の（１６）式に従って計算する。
【０１４１】
【数１５】
0005851205.tif 000016
【０１４２】
上記（１６）式の意味するところは、全時刻の観測量と1時刻後のクラスタ割り当てzt+1が与えられた時に時刻tのクラスタ割り当てがkとなる確率である。I(c)は、条件cが成立した場合に1,そうでない場合には0を返す関数である。なお、最初にt=Tを決定する際には、zT+1の値を便宜的に固定しておくか、あるいは以下の（１７）式で計算する。
【０１４３】
【数１６】
0005851205.tif 000017
【０１４４】
次のステップＳ１３４では、変数ｋを１インクリメントして、ステップＳ１３５において、変数ｋがＫより大きいか否かを判定する。変数ｋがＫ以下である場合には、上記ステップＳ１３３へ戻る。一方、変数ｋがＫより大きい場合には、全てのｋについてｚt＝ｋとなる確率を計算したとして、ステップＳ１３６へ進む。
【０１４５】
そして、ステップＳ１３６では、上記ステップＳ１３３で計算された全てのｋについてｚt＝ｋとなる確率（上記（１６）式の右辺の量に比例する確率）でztの値をk∈{１，２，．．．，Ｋ}より一つ選択するようにサンプリングする。次のステップＳ１３７では、時刻ｔを１デクリメントして、ステップＳ１３８において、時刻ｔが１より小さいか否かを判定する。時刻ｔが１以上である場合には、上記ステップＳ１３２へ戻る。一方、時刻ｔが１より小さい場合には、全てのｔについてｚｒtをサンプリングしたとして、ステップＳ１３９へ進む。
【０１４６】
このように、t=T，T-1，．．．，1の順番でztの値を決定する。
【０１４７】
続いて、ステップＳ１３９において、サンプリングによって不要となったクラスタの整理を行う。z1:Tをサンプリングした結果、zt=kとなる時刻インデックスt が一つもないクラスタが発生しうる。これは、今回の推論の結果、クラスタkはデータを表現する上で不要であると判断されたことを意味する。そこで、そのようなクラスタは全てメモリから消去し、同時にＫに依存する変数を全て更新する。
【０１４８】
例えば、utのサンプリング後にK=5となり、ztのサンプリング後にk=3となるztが一つもなかったとする。このとき、以下の処理を行い、潜在クラスタ数を一つ減らす。
【０１４９】
１．クラスタ割り当て変数Zのうち、zt>3となる全てのztの値を1小さくして更新
２．特徴選択変数Rのうち、k>3となるrk，iの値をrk-1，i に上書きして更新し、rK，i はメモリから消去する
３．特徴観測パラメータθk，iのうち、k>3となるθk，iの値をθk-1，iに上書きして更新し、θK，iはメモリから消去する
４．クラスタ遷移パラメータπk，lのうち、k>3となるkπ，lの値をπk-1，lに上書きして更新し、l>3となるπk，l-１の値をπk，l_-1に上書きして更新し、θK+1，l,θk，K+1は全てメモリから消去する
５．混合比パラメータβkのうち、k>3となるβkの値をβk-1に上書きして更新し、βＫ←１−ΣK-1k=1βkと更新し、βＫK+1はメモリから消去する
６．K←K-1
【０１５０】
以上と同等の処理を、使われていないクラスタ全てに対して繰り返す。
【０１５１】
上記ステップＳ１３９におけるクラスタの整理手順のフローチャートを図６に示す。
【０１５２】
ステップＳ１８０において、変数ｋをＫに設定し、ステップＳ１８１において、時刻ｔを初期値１に設定する。ステップＳ１８２では、ｚt＝ｋであるか否かを判定する。ｚt＝ｋである場合には、クラスタｋは使われていると判断し、後述するステップＳ２１１へ移行する。一方、ｚt＝ｋでない場合には、ステップＳ１８３において、時刻ｔを１インクリメントして、ステップＳ１８４において、時刻ｔがＴより大きいか否かを判定する。時刻ｔがＴ以下である場合には、上記ステップＳ１８２へ戻る。一方、時刻ｔがＴより大きい場合には、クラスタｋが使われないと判断し、ステップＳ１８５へ進む。
【０１５３】
ステップＳ１８５では、ｋ*にｋを設定し、ステップＳ１８６において、時刻ｔに１を設定する。ステップＳ１８７において、ｚｔ＞ｋ*であるか否かを判定し、ｚｔ＞ｋ*であれば、ステップＳ１８８で、zt>ｋ*となるztの値を1小さくして更新する。一方、ｚｔ＞ｋ*でない場合には、ステップＳ１８９へ移行する。
【０１５４】
ステップＳ１８９では、時刻ｔを１インクリメントして、ステップＳ１９０において、時刻ｔが、時刻ｔがＴより大きいか否かを判定する。時刻ｔがＴ以下である場合には、上記ステップＳ１８６へ戻る。一方、時刻ｔがＴより大きい場合には、zt>ｋ*となるztの値を更新したとして、ステップＳ１９１へ進む。
【０１５５】
ステップＳ１９１では、変数ｋにｋ*を設定し、ステップＳ１９２において、ｉ＝１，..．，Ｎに対しrk，iの値をrk-1，i に上書きして更新する。ステップＳ１９３では、ｉ＝１，..．，Ｎに対しθk，iの値をθk-1，i に上書きして更新する。
【０１５６】
そして、ステップＳ１９４において、変数ｋを１インクリメントして、ステップＳ１９６で、変数ｋがＫに達したか否かを判定する。変数ｋがＫに達していない場合には、ステップＳ１９２へ戻るが、一方、変数ｋがＫに達した場合には、k>k*となるrk，iの値及びθk，iの値を更新したとして、ステップＳ１９７へ移行する。
【０１５７】
ステップＳ１９７では、ｉ＝１，..．，Ｎに対しrK，iをメモリから消去する。次のステップＳ１９８では、ｉ＝１，..．，Ｎに対し、θK，iをメモリから消去する。
【０１５８】
そして、ステップＳ１９９では、変数ｋにｋ*を設定し、ステップＳ２００において、ｌ＝１，..．，Ｋ＋１に対しπk，lの値をπk-1，lに上書きして更新する。ステップＳ２０１では、ｌ＝１，..．，Ｋ＋１に対しπl，kの値をπl-1，kに上書きして更新する。
【０１５９】
ステップＳ２０２において、βkの値をβk-1に上書きして更新する。ステップＳ２０３では、ｋを１インクリメントして、ステップＳ２０４において、変数ｋがＫ＋１に達したか否かを判定する。変数ｋがＫ＋１に達していない場合には、ステップＳ２００へ戻るが、一方、変数ｋがＫ＋１に達した場合には、k>k*となるπl-１，kの値、πk，l-１の値、及びβk-1の値を更新したとして、ステップＳ２０５へ移行する。
【０１６０】
ステップＳ２０５では、ｌ＝１，..．，Ｋ＋１に対しπK+1，lをメモリから消去する。次のステップＳ２０６では、ｌ＝１，..．，Ｋ＋１に対し、πl，K+1をメモリから消去する。ステップＳ２０７では、βＫ←１−ΣK-1k=1βkと更新し、ステップＳ２０８において、βＫK+1はメモリから消去する。
【０１６１】
次のステップＳ２０９において、潜在クラスタ数Ｋを１つ減らして更新し、ステップＳ２０１において、ｋにｋ*を設定して、ステップＳ２１１へ移行する。
【０１６２】
ステップＳ２１１では、変数ｋを１デクリメントして、ステップＳ２１２において、変数ｋが、１より小さいか否かを判定する。変数ｋが１以上である場合には、上記ステップＳ１８１へ戻る。一方、変数ｋが１より小さい場合には、使われていないクラスタ全てを整理したとして、上記ステップＳ１３９の処理を終了する。
【０１６３】
上記図４のフローチャートのステップＳ１４０では、変数ｋを１に設定し、ステップＳ１４１において、変数ｌを１に設定する。そして、ステップＳ１４２において、以下の（１８）式に従って、補助変数mk，lを計算する。なお、mk，lは全時刻を通じてクラスタkからクラスタlに遷移した回数を表す。
【０１６４】
【数１７】
0005851205.tif 000018
【０１６５】
そして、ステップＳ１４３では、補助変数Sk，lを確率的に以下の（１９）式でサンプリングする。なお、Sk，lは混合比パラメータ３１０の更新などに必要となる。
【０１６６】
【数１８】
0005851205.tif 000019
【０１６７】
ここでs(・，・) はStirling number of the first kindと呼ばれる関数である．
【０１６８】
次のステップＳ１４４では、変数ｌを１インクリメントして、ステップＳ１４５において、変数ｌがＫより大きいか否かを判定する。変数ｌがＫ以下である場合には、上記ステップＳ１４２へ戻る。一方、変数ｌがＫより大きい場合には、ステップＳ１４６において、変数ｋを１インクリメントする。そして、ステップＳ１４７において、変数ｋがＫより大きいか否かを判定する。変数ｋがＫ以下である場合には、上記ステップＳ１４１へ戻る。変数ｋがＫより大きい場合には、1≦k，l≦Kについて、補助変数mk，l、Sk，lを計算したとして、ステップＳ１４８へ進む。
【０１６９】
ステップＳ１４８では、更新された全ての変数、パラメータ、補助変数を記憶部３の対応する部分に設定して、クラスタ割り当て変数の更新が終了する。
【０１７０】
＜特徴選択変数更新部２４＞
特徴選択変数更新部２４では、特徴選択変数３５を更新する。
【０１７１】
特徴選択変数更新部２４の更新手順のフローチャートを図７に示す。
【０１７２】
まず、ステップＳ２２２において、記憶部３により、T, N, K, x, z, r, φ, θ, e, fをロードする。そして、ステップＳ２２３において、変数ｋを１に設定し、ステップＳ２２４において、変数ｉを１に設定する。
【０１７３】
そして、ステップＳ２２５において、以下の（２０）式、（２１）式を利用してrk，iの値を0か1かに更新する。
【０１７４】
【数１９】
0005851205.tif 000020
【０１７５】
上記の式で利用されている記号は以下のとおり。
【０１７６】
【数２０】
0005851205.tif 000021
【０１７７】
上記（２０）式の右辺は、特徴量ｉがクラスタｋに対して重要でないと仮定した場合の適合度に対応し、上記（２１）式の右辺は、特徴量ｉがクラスタｋに対して重要であると仮定した場合の適合度に対応している。
【０１７８】
そして、ステップＳ２２６において、変数ｉを１インクリメントし、ステップＳ２２７において、変数ｉがＮより大きいか否かを判定する。変数ｉがＮ以下である場合には、上記ステップＳ２２５へ戻る。一方、変数ｉがＮより大きい場合には、全てのｉについてｒk，ｉをサンプリングしたとして、ステップＳ２２８へ進む。
【０１７９】
ステップＳ２２８では、変数ｋを１インクリメントし、ステップＳ２２９において、変数ｋがＫより大きいか否かを判定する。変数ｋがＫ以下である場合には、上記ステップＳ２２４へ戻る。一方、変数ｋがＫより大きい場合には、k=1，．．．，K、i=1，．．．，Nの全てについてｒk，ｉをサンプリングしたとして、ステップＳ２３０へ進む。
【０１８０】
ステップＳ２３０では、更新されたrk，iを記憶部３の特徴選択変数３５に格納し、更新を終了する。
【０１８１】
＜特徴選択パラメータ更新部２５＞
特徴選択パラメータ更新部２５では、特徴選択パラメータ３６を更新する。
【０１８２】
特徴選択パラメータ更新部２５の更新手順のフローチャートを図８に示す。
【０１８３】
まず、ステップＳ２３１において、記憶部３により、N, K, r, e, fをロードする。そして、ステップＳ２３２において、変数ｉを１に設定し、ステップＳ２３３において以下の（２６）式を利用してλiの値をサンプリングして更新する。
【０１８４】
【数２１】
0005851205.tif 000022
【０１８５】
そして、ステップＳ２３４において、変数ｉを１インクリメントし、ステップＳ２３５において、変数ｉがＮより大きいか否かを判定する。変数ｉがＮ以下である場合には、上記ステップＳ２３３へ戻る。一方、変数ｉがＮより大きい場合には、i=1，．．．，Nの全てについてλｉの値を更新したとして、ステップＳ２３６へ進む。
【０１８６】
ステップＳ２３６では、更新されたλiを記憶部３の特徴選択パラメータ３６に格納し、更新を終了する。
【０１８７】
＜特徴観測パラメータ更新部２６＞
特徴観測パラメータ更新部２６では、特徴観測パラメータ３７を更新する。
【０１８８】
特徴観測パラメータ更新部２６の更新手順のフローチャートを図９に示す。
【０１８９】
まず、ステップＳ２４１において、記憶部３により、T, N, K, x, z, r, c, dをロードする。そして、ステップＳ２４２において、変数ｋを１に設定し、ステップＳ２４３において、変数ｉを１に設定する。
【０１９０】
そして、ステップＳ２４４において、以下の（２７）式を利用してθk，iの値をサンプリングして更新する。
【０１９１】
【数２２】
0005851205.tif 000023
【０１９２】
ただし、nk，i，Nk，iはそれぞれ上記（２４）式、（２５）式で定義されたものと同様である。
【０１９３】
そして、ステップＳ２４５において、変数ｉを１インクリメントし、ステップＳ２４６において、変数ｉがＮより大きいか否かを判定する。変数ｉがＮ以下である場合には、上記ステップＳ２４４へ戻る。一方、変数ｉがＮより大きい場合には、全てのｉについてθk，ｉをサンプリングしたとして、ステップＳ２４７へ進む。
【０１９４】
ステップＳ２４７では、変数ｋを１インクリメントし、ステップＳ２４８において、変数ｋがＫより大きいか否かを判定する。変数ｋがＫ以下である場合には、上記ステップＳ２４３へ戻る。一方、変数ｋがＫより大きい場合には、k=1，．．．，K、i=1，．．．，Nの全てについてθk，ｉをサンプリングしたとして、ステップＳ２４９へ進む。
【０１９５】
ステップＳ２４９では、更新されたθk，iを記憶部３の特徴観測パラメータ３７に格納し、更新を終了する。
【０１９６】
＜ノイズ観測パラメータ更新部２７＞
ノイズ観測パラメータ更新部２７では、ノイズ観測パラメータ３８を更新する。
【０１９７】
ノイズ観測パラメータ更新部２７の更新手順のフローチャートを図１０に示す。
【０１９８】
まず、ステップＳ２５１において、記憶部３により、N, K, x, z, r, a, bをロードする。そして、ステップＳ２５２において、変数ｉを１に設定し、ステップＳ２５３において以下の（２８）式〜（３０)式を利用して、φiの値をサンプリングして更新する。
【０１９９】
【数２３】
0005851205.tif 000024
【０２００】
ただし、nk，i，Nk，iはそれぞれ上記（２４）式、（２５）式で定義されたものと同様である。
【０２０１】
そして、ステップＳ２５４において、変数ｉを１インクリメントし、ステップＳ２５５において、変数ｉがＮより大きいか否かを判定する。変数ｉがＮ以下である場合には、上記ステップＳ２５３へ戻る。一方、変数ｉがＮより大きい場合には、i=1，．．．，Nの全てについてφｉの値を更新したとして、ステップＳ２５６へ進む。
【０２０２】
ステップＳ２５６では、更新されたφiを記憶部３のノイズ観測パラメータ３８に格納し、更新を終了する。
【０２０３】
＜クラスタ遷移パラメータ更新部２８＞
クラスタ遷移パラメータ更新部２８では、クラスタ遷移パラメータ３９を更新する。
【０２０４】
クラスタ遷移パラメータ更新部２８の更新手順のフローチャートを図１１に示す。
【０２０５】
まず、ステップＳ２６１において、記憶部３により、K, α, β, mをロードする。そして、ステップＳ２６２において、変数ｋを１に設定し、ステップＳ２６３において、補助変数mk，lを利用して、以下の（３１）式にしたがって、πkの値をサンプリングして更新する。
【０２０６】
【数２４】
0005851205.tif 000025
【０２０７】
そして、ステップＳ２６４において、変数ｋを１インクリメントし、ステップＳ２６５において、変数ｋがＫ＋１より大きいか否かを判定する。変数ｋがＫ＋１以下である場合には、上記ステップＳ２６３へ戻る。一方、変数ｋがＫ＋１より大きい場合には、ｋ=1，．．．，Ｋ＋１の全てについてπkの値を更新したとして、ステップＳ２６６へ進む。
【０２０８】
ステップＳ２６６では、更新されたπkを記憶部３のクラスタ遷移パラメータ３９に格納し、更新を終了する。
【０２０９】
＜混合比パラメータ更新部２９＞
混合比パラメータ更新部２９では、混合比パラメータ３１０を更新する。
【０２１０】
混合比パラメータ更新部２９の更新手順のフローチャートを図１２に示す。
【０２１１】
まず、ステップＳ２７１において、記憶部３により、K, γ,Ｓをロードする。そして、ステップＳ２７２において、補助変数Sk，lを利用して、以下の（３３）式に従って、K+1次元ベクトルをサンプリングして更新する。
【０２１２】
【数２５】
0005851205.tif 000026
【０２１３】
ステップＳ２７３では、更新されたβを記憶部３の混合比パラメータ３１０に格納し、更新を終了する。
【０２１４】
＜ハイパーパラメータ更新部２１０＞
ハイパーパラメータ更新部２１０では、ハイパーパラメータ３１２の各要素を更新する。
【０２１５】
ここでは、超事前分布としてガンマ分布を仮定し、その事後分布からハイパーパラメータをサンプリングする方法を記載する。なお、以下は簡単のため超事前分布のガンマ分布のパラメータは省略しGamma(・)とのみ表記する。
【０２１６】
ハイパーパラメータ更新部２１０の更新手順のフローチャートを図１３に示す。
【０２１７】
まず、ステップＳ２８１において、記憶部３により、T, K, N, r, θ, φ, γ, α, β, a, b, c, d, e, f, Sをロードする。そして、ステップＳ２８２において、以下の（３４)式を利用して、γの値をサンプリングして更新する。
【０２１８】
【数２６】
0005851205.tif 000027
【０２１９】
次のステップＳ２８３では、以下の（３５)式を利用して、αの値をサンプリングして更新する。
【０２２０】
【数２７】
0005851205.tif 000028
【０２２１】
そして、ステップＳ２８４において、変数ｉを１に設定し、ステップＳ２８５において以下の（３６）式を利用してａiの値をサンプリングして更新する。
【０２２２】
【数２８】
0005851205.tif 000029
【０２２３】
ここでB(・，・) はベータ関数である。
【０２２４】
次のステップＳ２８６において以下の（３７）式を利用してｂiの値をサンプリングして更新する。
【０２２５】
【数２９】
0005851205.tif 000030
【０２２６】
ステップＳ２８７では、以下の（３８）式を利用してｃiの値をサンプリングして更新する。
【０２２７】
【数３０】
0005851205.tif 000031
【０２２８】
ステップＳ２８８では、以下の（３９）式を利用してｄiの値をサンプリングして更新する。
【０２２９】
【数３１】
0005851205.tif 000032
【０２３０】
ステップＳ２８９では、以下の（４０）式を利用してｅiの値をサンプリングして更新する。
【０２３１】
【数３２】
0005851205.tif 000033
【０２３２】
ステップＳ２９０では、以下の（４１）式を利用してｆiの値をサンプリングして更新する。
【０２３３】
【数３３】
0005851205.tif 000034
【０２３４】
そして、ステップＳ２９１において、変数ｉを１インクリメントし、ステップＳ２９２において、変数ｉがＮより大きいか否かを判定する。変数ｉがＮ以下である場合には、上記ステップＳ２８５へ戻る。一方、変数ｉがＮより大きい場合には、i=1，．．．，Nの全てについてａｉ，ｂｉ，ｃｉ，ｄｉ，ｅｉ，ｆｉの値を更新したとして、ステップＳ２９３へ進む。
【０２３５】
ステップＳ２９３では、更新されたハイパーパラメータを記憶部３のハイパーパラメータ３１２に格納し、更新を終了する。
【０２３６】
なお、ハイパーパラメータの更新には、他にも様々な方法を利用することができる。また、ハイパーパラメータは必ずしも更新しなくても良い。
【０２３７】
＜解析終了条件判定部２１１＞
解析終了条件判定部２１１は、変数計算部22を監視し、変数の更新処理を続けるか終了するかを判定する。
【０２３８】
一般に、変数の更新は繰り返し処理により最適値へ近づく。しかし、いつ繰り返しを止めるかについては様々な方法がある。例えば、変数更新後に記憶部３に格納された変数やパラメータを使って、与えられたデータの尤度を計測し, それが一定値あるいは一定割合以上の変化をしなくなったら終了とする、あるいは単純に変数計算部２２の繰り返し回数が一定値以上になったら更新を終了するなどである。
【０２３９】
以下の実験結果では、最も簡単に、Ｌを繰り返し処理の最大値として入力し、変数計算部の更新処理をＬ回繰り返した時点で解析処理を終了することとした。実際にはどのような終了条件を決定しても構わない。
【０２４０】
設定した解析終了条件が満たされていない場合、再び変数計算部２２に処理を戻す。設定した解析終了条件が満たされると、解析結果はユーザが規程した所望の形で出力部４へと出力される。
【０２４１】
＜実験例＞
次に、本実施の形態で提案する時系列データのクラスタリング手法を用いた実験の結果について説明する。実験では、３種類、計５つのデータセットに対して実験を行った。
【０２４２】
データの種類は以下の３種類である。
【０２４３】
１．人工データ(synth1, synth2, synth3): 提案法の想定する時系列データの生成過程に従い、人工的に二値のバイナリデータを生成した。”synth1”はN=200、T=50、”synth2”はN=100、T=50、“synth3”はN=100、T=70であり、それぞれ想定する構造に違いがある。
【０２４４】
２.文書データ(PNAS):Prceedings of the National Academy of Sciences(PNAS)と呼ばれる学術論文誌のデータを利用した。1915年から2005年までの医学・薬学分野の論文タイトルに出現する単語の頻度を年度ごとに集計し行列データとして解析する。年度tで50回以上出現した単語iの観測量をxi，t=1とし、50回に満たない単語はxi，t=0とする。T=91年に渡って一度でもxi，t=1となるN=118種の単語の0/1データである。
【０２４５】
３．音楽データ(instrument):およそ11秒の合成音楽の波形を解析した。音楽波形に対して移動窓を使った短時間フーリエ解析を行い、T=742個フレームにおいてN=513の周波数帯におけるスペクトル値を算出する。xi，tとしては、フレームt、周波数帯iにおけるスペクトル値の自然対数(log)を計算する。
【０２４６】
対数スペクトル値が一定の閾値を越えていればxi，t=1、そうでなければxi，t=0とした。
【０２４７】
定量評価として、テストデータ対数尤度を用いた。与えられたN×T個の観測データのうちおよそ5%の観測データを未知のテストデータとしてランダムに選択肢、変数更新部の処理から除外する。解析終了後、これらテストデータの対数尤度を計算し、それをテストデータ数で平均した。
【０２４８】
比較対象として、特徴選択変数rを含まない、通常のiHMM（上記非特許文献３、４を参照）を利用した。提案法および比較対象であるiHMMともに複数のハイパーパラメータの組み合わせでそれぞれ２０回試行を行い、最も数値の良かったハイパーパラメータを選択した。
【０２４９】
実験結果を以下の表1にまとめた。
【０２５０】
【表１】
0005851205.tif 000035
対数尤度は数字が大きいほど良い。この表１より分かるように、特徴選択の機能が加わった提案法は大部分のデータで既存手法より正確に時系列データのモデル化が可能であることが示された。
【０２５１】
最後に定性的な評価としてPNASデータの解析によって得られた、年度ごとの” 特徴” 単語を図１４に示す。
【０２５２】
各コラムはクラスタを表わし、クラスタ上部の数字はクラスタの対応する年度（時間ステップt)であり、下部の単語列は特徴として選ばれた単語である。
【０２５３】
これによれば、例えば“vitro” という単語は1980年ごろ(左から2番目のコラム)に現れ1996年ごろ(右から2番目のコラム) には論文のタイトルとして現れなくなった。それと入れかわって、“vivo”というという単語が1991年ごろ(左から4番目のコラム)に現れた。これらは”in vitro” と”in vivo” という生化学・医学研究における重要な概念の出現と変化に符合している。また、分子生物学を表すと思われる“molecular” という単語は1985年ごろ(左から3番目のコラム) に現れさらに遺伝子学や遺伝子治療に関連する。“genetic” という単語は1996年ごろ(右から2番目のコラム) に現れる。このように、各年代で大きく話題になったトピックの単語を正しく各クラスタの特徴として抽出できている．これにより、例えば医学・生理学分野の技術トレンドを少数のクラスタの選択された特徴単語から簡単に調査することができる。
【０２５４】
以上説明したように、本発明の実施の形態に係るクラスタリング装置によれば、クラスタ割り当て変数ｚtを更新するときに潜在クラスタ数Ｋを決定すると共に、クラスタ割り当て変数ｚtの更新、特徴選択変数ｒk,iの更新、特徴選択パラメータの更新λi、特徴観測パラメータの更新θk,i、ノイズ観測パラメータの更新φi、クラスタ遷移パラメータπk,lの更新、混合比パラメータβk、及びハイパーパラメータの更新を繰り返すことにより、観測値の時系列データＸが与えられたときに、クラスタ数を自動的に決定すると共に、時刻のクラスタリングに寄与する真の特徴量と、時刻のクラスリングに寄与しないノイズ特徴の影響を自動的に推定・分離して、時刻のクラスタリングを行うことができる。
【０２５５】
また、高次元の特徴データを複数の似た特性を示す時間帯へとクラスタリングし、複雑な時間変化を少数のパターンの変化として直感的に解析し、データのトレンドを簡単に調査・可視化できる。
【０２５６】
また、既存の特徴選択を含まない時系列データモデリングよりも高いデータ表現精度を持つことができる。
【０２５７】
また、与えられた観測量の中から、時刻のクラスタリングに寄与する真の特徴量と、クラスタ化に寄与しないノイズ特徴量の存在をモデル化する。これによって、与えられた観測値の時系列データをクラスタリングする上で意味のある情報と意味の無い情報を分離することができるので、クラスタの推定時にノイズ特徴の影響を排除することが可能となる。
【０２５８】
また、真の特徴量とノイズ特徴量は各クラスタごとに固有に定義される。このことで、動的に特徴の分布が変化する時系列データに対してより精度の良いモデル化が可能となる。また、各時刻を特徴づける特徴量を自動的に選択・抽出するので、与えられた時系列データが高次元の特徴量を持っていてもその定性的な性質を容易に理解することが可能となる。
【０２５９】
また、観測量を重要な特徴、ノイズ観測量へと分類することができるので、各時刻において本質的に違いを生み出した情報を自動的に抽出して提供可能である。
【０２６０】
また、クラスタの数は自動的にデータに応じて決定される。これまで多くの手法では探索的に決定することが多かった。このため、手法利用時の工数を減らすことができる。
【０２６１】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０２６２】
例えば、上述のクラスタリング装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【０２６３】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【０２６４】
１００ クラスタリング装置
１     入力部
２     解析部
３     記憶部
２１   初期設定部
２２   変数計算部
２３   変数更新部
２４   特徴選択変数更新部
２５   特徴選択パラメータ更新部
２６   特徴観測パラメータ更新部
２７   ノイズ観測パラメータ更新部
２８   クラスタ遷移パラメータ更新部
２９   混合比パラメータ更新部
２１０ ハイパーパラメータ更新部
２１１ 解析終了条件判定部
【図１】
0005851205.tif 000039
【図２】
0005851205.tif 000040
【図３】
0005851205.tif 000041
【図４】
0005851205.tif 000042
【図５】
0005851205.tif 000043
【図６】
0005851205.tif 000044
【図７】
0005851205.tif 000045
【図８】
0005851205.tif 000046
【図９】
0005851205.tif 000047
【図１０】
0005851205.tif 000048
【図１１】
0005851205.tif 000049
【図１２】
0005851205.tif 000050
【図１３】
0005851205.tif 000051
【図１４】
0005851205.tif 000052