クラスタリング計算装置、クラスタリング計算方法、クラスタリング計算プログラム並びにそのプログラムを記録したコンピュータ読み取り可能な記録媒体
【課題】ダイアライゼーションにおいて、話者数と各話者を特徴づけるパラメータを正確に推定するクラスタリング技術を提供する。
【解決手段】クラスタリング計算装置３は、録音された会話データから抽出されたパワーベクトルを読み込んでｄＨＤＰ近似モデルに対応したベクトルの観測量に変換する観測量生成部２６と、変換された観測量の集合データを蓄積記憶する記憶手段１０と、観測量の集合データから複数のクラスタをｄＨＤＰ近似モデルにより生成するときの複数のパラメータの事後分布の値をＥＭアルゴリズムによりそれぞれ推定し、記憶手段１０に順次蓄積および更新する変分事後分布推論部３０と、事前に設定された終了条件が成立したときに、記憶手段１０に記憶されている複数のパラメータの事後分布の最新の推定値を出力する出力制御部２８とを備える。
【解決手段】クラスタリング計算装置３は、録音された会話データから抽出されたパワーベクトルを読み込んでｄＨＤＰ近似モデルに対応したベクトルの観測量に変換する観測量生成部２６と、変換された観測量の集合データを蓄積記憶する記憶手段１０と、観測量の集合データから複数のクラスタをｄＨＤＰ近似モデルにより生成するときの複数のパラメータの事後分布の値をＥＭアルゴリズムによりそれぞれ推定し、記憶手段１０に順次蓄積および更新する変分事後分布推論部３０と、事前に設定された終了条件が成立したときに、記憶手段１０に記憶されている複数のパラメータの事後分布の最新の推定値を出力する出力制御部２８とを備える。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、話者数が未知である会話の録音データから話者数を推定する技術に関する。
【背景技術】
【０００２】
従来、複数話者の会話からなる録音データから、その会話に参加した話者数とさらに各話者の発話したタイミングを推定する問題が知られている。この問題やこの問題を解決する技術は、ダイアライゼーション（diarization）と呼ばれている。ダイアライゼーションは、簡単には「いつ、誰が発話したか」を自動推定する技術といえる。この技術の応用としては、会議データへのアノテーション（annotation）やそれを用いた検索、音声強調など幅広い用途が期待されている。
【０００３】
既存のダイアライゼーションの解法は、大きく２つに分けることができる。第１の方法は、話者固有の声質を推定して話者を区別する方法であり、第２の方法は、話者の位置を推定する方法である。このうち、第１の方法（話者声質を利用する方法）では、各話者の声の特徴を抽出することで現在発話している話者を識別する。この第１の方法では、話者が移動したとしても話者を識別できるという利点があるが、複数の話者が同時に発話した場合には話者の識別が困難となるという欠点がある。
【０００４】
第２の方法（話者位置に関する情報を利用する方法）は、話者の位置を推定することで話者数とその位置を推定する方法である（例えば、非特許文献１〜３参照）。非特許文献１〜３の手法では、マイクロホンアレーを用いて各話者の位置を推定することで話者の識別を行う。従って、非特許文献１〜３の手法は、話者が移動すると同一話者として識別することができないという欠点があるが、複数話者が同時に発話する場合でも各話者の発話行動を識別することができるという利点がある。
【０００５】
ダイアライゼーションでは、一般に、録音された音声データ中の発話者数は未知であり、推定しなければならない。また、第１の方法（話者声質を利用する方法）によるもの、第２の方法（話者位置に関する情報を利用する方法）によるものに関わらず、各話者を識別するためには、各話者を特徴づける量（パラメータ）を推定する必要がある。これは、いわゆるクラスタリングの問題に相当する。クラスタリングの問題とは、観測されたデータを適切な数のクラスタに分割し、各クラスタのパラメータを推定する問題である。ダイアライゼーションを単純にクラスタリングの問題に置き換えようとしても、そのクラスタ数とデータの分割は未知なので、パラメータと共に学習する必要がある。ダイアライゼーションの問題では、クラスタ数は「話者数」を表し、各クラスタのパラメータは「話者の特徴量」を表し、結果得られるデータの分割結果から「各話者の発話タイミング」が示唆されることになる。
【０００６】
非特許文献１ではleader-followerアルゴリズムと呼ばれるクラスタリング手法を用いている。これは、逐次的なクラスタリング手法であり、少しずつサンプルを入力していきながらクラスタリングを行うものである。このクラスタリング手法では、クラスタリングの際に、新しく入力されたサンプルが既存のどのクラスタからも一定距離以上離れていた場合に、そのサンプルを中心として新しいクラスタを生成する。なお、leader-followerアルゴリズムについては、「R. O. Duda, P. E. Hart and D. G. Stork, “Pattern Classification”, John Wiley & Sons, 2001.」に記載されている。
【０００７】
また、非特許文献３では、ＢＩＣ基準と呼ばれる評価値を最大化するようなクラスタ数を選択する方法を提案している。この方法では、クラスタ数をＫに設定してから実際にクラスタリングを行って評価値を計算し、次はクラスタ数をＫ＋１に設定して評価値を計算する、ということを繰り返えして評価値を最大化するクラスタ数を探索する。
【０００８】
ところで、クラスタ数未知のデータに対するクラスタリングモデルとして、ノンパラメトリックベイズモデルが、近年、多くの場面で利用されるようになってきた。例えば、ノンパラメトリックベイズモデルの１種であるDirichlet Process Mixture（ＤＰＭ）は、クラスタ数と各サンプルのクラスタリングを確率的に同時に最適化することができる。従って、ＤＰＭは、既存のクラスタリングモデルのように、クラスタ数の最適化を簡単に実現できる点が大きな特徴である。このＤＰＭにおいて、確率分布の連続的な時間変化をモデル化した拡張モデルとしては、dynamic Hierarchical Dirichlet Process（ｄＨＤＰ）と呼ばれるモデルが知られている（非特許文献４参照）。
【先行技術文献】
【非特許文献】
【０００９】
【非特許文献１】S. Araki, M. Fujimoto, K. Ishizuka, H. Sawada and S. Makino, ”A DOA Based Speaker Diarization System for Real Meetings”, Proceedings of the Joint Workshop on Hands-Free Speech Communication and Microphone Arrays, p.29-32, 2008.
【非特許文献２】荒木章子、藤本雅清、石塚健太郎、澤田宏、牧野昭二、“音声区間検出と方向情報を用いた会議音声話者識別システムとその評価”、日本音響学会講演論文集、春季、vol. 1-10-1, p. 1-4, 2008.
【非特許文献３】J. M. Pardo, X. Anguera and C. Wooters, “Speaker Diarization for Multi-Microphone Meetings Using Only Between-Channel Differences”, Proceedings of the Third Joint Workshop on Multimodal Interaction and Related machine Learning Algorithms”, p. 257-264, 2008.
【非特許文献４】L. Ren, D. B. Dunson and L. Carin, “The Dynamic Hierarchical Dirichlet process”, Proceedings of International Conference on Machine Learning, p. 824-831, 2008.
【発明の概要】
【発明が解決しようとする課題】
【００１０】
従来のダイアライゼーションの研究において、例えば、非特許文献１に記載の手法は、計算が簡単で高速に実行できるが、新しいクラスタを生成するための距離閾値の設定が必要である。この閾値の設定によって、最終的に得られるクラスタリングとクラスタ数が決定される。一方で、クラスタ数やクラスタリング結果の推定値を真の値に近づけるという意味で、この閾値を最適化することは困難である。
【００１１】
また、例えば、非特許文献３に記載の方法は、実際には不適切なクラスタ数のもとでもクラスタリングを行う必要があり、処理の過程において、計算量や時間の観点からは大きな無駄が発生する。
【００１２】
したがって、従来のダイアライゼーションの研究では、話者数と話者の特徴を推定する部分、つまりクラスタリングの問題に改良の余地があった。また、ダイアライゼーションに対して、処理負荷の低減や処理の高速化が要望されている。さらに、従来のダイアライゼーションの研究では、ｄＨＤＰやＤＰＭなどのノンパラメトリックベイズモデルは、クラスタリングモデルとして採用されておらず、その適用方法が知られていなかった。
【００１３】
そこで、本発明では、前記した問題を解決し、ダイアライゼーションにおいて、話者数と各話者を特徴づけるパラメータを正確に推定するクラスタリング技術を提供することを目的とする。
【課題を解決するための手段】
【００１４】
前記目的を達成するために、本願発明者らは、ダイアライゼーションにおいて、話者数と各話者を特徴づけるパラメータを推定するクラスタリングにおいて種々検討を行った。その結果、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化できるモデルとしてノンパラメトリックベイズモデルを採用したときに、話者数を正確に推定できることを見出した。
【００１５】
そこで、本発明に係るクラスタリング計算装置は、話者数が未知である会話の録音データから前記会話の話者数を推定するために、各話者を特徴付ける特徴量を抽出する特徴量抽出手段と、前記抽出された特徴量から前記各話者に対応する複数のクラスタを生成するときの複数の未知パラメータの値をそれぞれ推定するクラスタリング計算装置と、前記推定された複数のパラメータ値により前記会話の各話者を識別する識別手段とを有したダイアライゼーションシステムの前記クラスタリング計算装置であって、前記抽出された特徴量を読み込む読込手段と、前記読み込んだ特徴量をノンパラメトリックベイズモデルに対応したベクトルの観測量に変換する観測量生成手段と、前記変換された観測量の集合データを蓄積記憶する観測量記憶手段と、前記観測量の集合データから複数のクラスタをノンパラメトリックベイズモデルにより生成するときの複数のパラメータの事後分布の値をＥＭアルゴリズムによりそれぞれ推定および更新する事後分布推論手段と、前記推定および更新された複数のパラメータの事後分布の値を蓄積記憶する推定値記憶手段と、事前に設定された終了条件が成立したときに前記推定値記憶手段に記憶されている前記複数のパラメータの事後分布の最新の推定値を出力する出力制御手段と、を備えることを特徴とする。
【００１６】
かかる構成によれば、クラスタリング計算装置は、会話の録音データから抽出された特徴量をノンパラメトリックベイズモデルに適応できるようにベクトルの観測量に変換することで生成した観測量の集合データを用いて、ノンパラメトリックベイズモデルに従った演算により、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化する。そして、クラスタリング計算装置は、この最適化において、ＥＭアルゴリズムにより事後分布の推定および更新を繰り返し、蓄積記憶する。ここで、ＥＭアルゴリズムは、局所最適解の計算アルゴリズムであるため、計算を繰り返すことで、ある１つの解に収束させることができる。また、ノンパラメトリックベイズモデルでは、予め準備した最大クラスタ数の個数のクラスタが最適化処理の過程で、有効なクラスタ数の個数に収束し、他のクラスタは、クラスタの混合比がほぼ０となる。このクラスタリング結果が得られれば、各クラスタへ帰属するサンプル数（観測量の個数）が計算できる。また、この各クラスタへ帰属するサンプル数を用いて各クラスタの混合比を計算可能である。さらに、この混合比により、有効なクラスタ数を決定することができる。ここで、クラスタ数は、会話の録音データ中の話者に対応しているので、話者数が決定できるようになる。また、これらクラスタリング結果と共に、それを用いた各種変数の推定値も同様に蓄積記憶することができる。そして、クラスタリング計算装置は、終了条件が成立したときに蓄積記憶されている最新の推定値を出力する。
【００１７】
また、本発明に係るクラスタリング計算装置は、前記事後分布推論手段が、前記ＥＭアルゴリズムのＥステップの処理として、ｄＨＤＰ（dynamic Hierarchical Dirichlet Process）モデルにおいて予め定められた事前分布、観測分布および最大クラスタ数と、ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＭステップまでに推定された隠れ変数の事後分布の推定値とを読み込んで、過去分を含めた演算対象時刻ごと、かつ、クラスタごと、かつ、演算対象時刻別の全データごとに、前記クラスタと、混合比と、前記クラスタ分布の時間変化の程度を表す重みとに関するパラメータの事後分布の値を演算するＥステップ用計算手段と、前記ＥＭアルゴリズムのＭステップの処理として、前記ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＥステップまでに推定されたパラメータの事後分布の推定値とを読み込んで、過去分を含めた演算対象時刻ごと、かつ、演算対象時刻別の全データごとに、２種類の隠れ変数の事後分布の値を推定し、前記２種類の隠れ変数のうち、前記クラスタ分布の時間変化の程度を表す重みに関連付けられた第１隠れ変数の事後分布の値についてはクラスタごとに演算し、前記第１隠れ変数および前記混合比に関連付けられた第２隠れ変数の事後分布の値については、演算対象時刻から過去に遡及した時刻ごとに演算するＭステップ用計算手段と、前記Ｅステップの処理と前記Ｍステップの処理とを交互に予め定められた回数だけ繰り返し実行させる制御を行う収束判定手段と、を備えることを特徴とする。
【００１８】
かかる構成によれば、クラスタリング計算装置は、ノンパラメトリックベイズモデルのうちｄＨＤＰモデルに従った演算により、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化するためにＥＭアルゴリズムを用いる。また、ｄＨＤＰモデルでは、ある時刻ｔにおけるパラメータ分布を、前の時刻ｔ−１でのパラメータ分布と、当該時刻ｔでのパラメータ分布の変化量とを、時刻ｔにおける変化量の変化の割合（重み）で重ね合わせた分布として表現している。すなわち、ｄＨＤＰでは、時刻に依存して分布が少しずつ変化するというモデルとなっている。クラスタリング計算装置のアルゴリズムは、ある時刻ｔにおけるパラメータ分布を、時刻ｔに実際に発話した発話者を構成しているものとして構築した。そのため、ｄＨＤＰモデルは、話者の入れ替わり（turn-taking）によって観測量（サンプル）の分布が時間変化するダイアライゼーションのタスクにとって好都合の確率モデルとなった。したがって、話者の入れ替わりがある会話録音データから、未知の話者数と各話者を特徴づけるパラメータを正確に推定することができる。
【００１９】
また、本発明に係るクラスタリング計算装置は、前記Ｅステップ用計算手段が、観測量の集合データと、隠れ変数の事後分布の推定値とについて、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込み、推定対象の時刻から事前に設定された設定値分だけ遡及した過去の時刻までを演算対象時刻として、前記クラスタと、前記混合比と、前記重みとに関するパラメータの事後分布の値を演算し、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、事前に設定された再推定すべきクラスタの判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合に、前記混合比に関するパラメータの事後分布の値を演算し、前記Ｍステップ用計算手段が、観測量の集合データと、パラメータの事後分布の推定値とについて、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込み、推定対象の時刻から事前に設定された設定値分まで遡及した過去の時刻までを演算対象時刻として、前記第１隠れ変数および第２隠れ変数の事後分布の値を演算し、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、前記判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合に、前記第１隠れ変数の事後分布の値を演算することを特徴とする。
【００２０】
かかる構成によれば、クラスタリング計算装置は、ＥステップおよびＭステップにおいて、推定対象の時刻から事前に定められた過去の一時点まで遡った観測量の集合データと、その演算に必要な推定値とを読み込み、当該推定対象の時刻から事前に定められた設定値分だけ過去に遡った過去の時刻までを演算対象時刻として、推定値を求める演算を行う。したがって、事前に過去の一時点や設定値を定めることなく、時間ステップの進展とともに推定すべき変数の個数が増加していくだけである場合と比べて、処理負荷の低減や処理の高速化を実現できる。また、クラスタリング計算装置は、ＥステップおよびＭステップにおいて、クラスタごとの演算のためにクラスタを更新したときに、必要な場合にだけそのクラスタにおける推定値の再推定処理を行う。したがって、ｄＨＤＰモデルにおいて予め設定される最大クラスタ数の個数のクラスタに関して、推定値の再推定処理を毎回実行する場合に比べて、処理負荷の低減や処理の高速化を実現できる。
【００２１】
また、課題を解決するため、本発明に係るクラスタリング計算方法は、話者数が未知である会話の録音データから前記会話の話者数を推定するダイアライゼーションシステムにおいて、記憶手段と、読込手段と、観測量生成手段と、事後分布推論手段と、出力制御手段とを備えて、前記録音データから抽出された各話者を特徴付ける特徴量から前記各話者に対応する複数のクラスタを生成するときの複数の未知パラメータの値をそれぞれ推定するクラスタリング計算装置のクラスタリング計算方法であって、前記読込手段によって、前記抽出された特徴量を読み込む特徴量読込ステップと、前記観測量生成手段によって、前記読み込んだ特徴量をノンパラメトリックベイズモデルに対応したベクトルの観測量に変換し、前記変換された観測量の集合データを記憶手段に順次蓄積する観測量蓄積ステップと、前記事後分布推論手段によって、前記観測量の集合データから複数のクラスタをノンパラメトリックベイズモデルにより生成するときの複数のパラメータの事後分布の値をＥＭアルゴリズムによりそれぞれ推定し、当該推定値を前記記憶手段に順次格納および更新する事後分布推定ステップと、前記出力制御手段によって、事前に設定された終了条件が成立したときに前記記憶手段に記憶されている前記複数のパラメータの事後分布の最新の推定値を出力する推定値出力ステップと、を含んで実行することを特徴とする。
【００２２】
かかる手順によれば、クラスタリング計算方法において、クラスタリング計算装置は、まず、会話の録音データから抽出された特徴量を読み込み、ノンパラメトリックベイズモデルに適応できるようにベクトルの観測量に変換し、蓄積する。そして、クラスタリング計算装置は、生成した観測量の集合データを用いて、ノンパラメトリックベイズモデルに従った演算により、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化する。そして、クラスタリング計算装置は、この最適化において、ＥＭアルゴリズムにより事後分布の推定および更新を繰り返し、蓄積記憶する。そして、クラスタリング計算装置は、終了条件が成立したときに蓄積記憶されている最新の推定値を出力する。
【００２３】
また、本発明に係るクラスタリング計算方法は、前記事後分布推論手段が、前記事後分布推定ステップにおいて、前記ＥＭアルゴリズムのＥステップの処理として、ｄＨＤＰモデルにおいて予め定められた事前分布、観測分布および最大クラスタ数と、ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＭステップまでに推定された隠れ変数の事後分布の推定値とを読み込んで、推定対象の時刻に過去分を含めた演算対象時刻ごと、かつ、クラスタごと、かつ、前記演算対象時刻別の全データごとに、前記クラスタと、混合比と、前記クラスタ分布の時間変化の程度を表す重みとに関するパラメータの事後分布の値を演算する段階と、前記ＥＭアルゴリズムのＭステップの処理として、前記ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＥステップまでに推定されたパラメータの事後分布の推定値とを読み込んで、推定対象の時刻に過去分を含めた演算対象時刻ごと、かつ、前記演算対象時刻別の全データごとに、２種類の隠れ変数の事後分布の値を推定し、前記２種類の隠れ変数のうち、前記クラスタ分布の時間変化の程度を表す重みに関連付けられた第１隠れ変数の事後分布の値についてはクラスタごとに演算し、前記第１隠れ変数および前記混合比に関連付けられた第２隠れ変数の事後分布の値については、前記演算対象時刻から過去に遡及した時刻ごとに演算する段階とを含み、前記Ｅステップの処理と前記Ｍステップの処理とを交互に予め定められた回数だけ繰り返し実行することを特徴とする。
【００２４】
かかる手順によれば、クラスタリング計算方法において、クラスタリング計算装置は、ノンパラメトリックベイズモデルのうちｄＨＤＰモデルに従った演算により、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化するためにＥＭアルゴリズムを用いる。ここで、ｄＨＤＰでは、時刻に依存して分布が少しずつ変化するというモデルとなっているため、話者の入れ替わり（turn-taking）によって観測量（サンプル）の分布が時間変化するダイアライゼーションのタスクにとって恰好の確率モデルである。したがって、話者の入れ替わりがある会話録音データから、未知の話者数と各話者を特徴づけるパラメータを正確に推定することができる。
【００２５】
また、本発明に係るクラスタリング計算方法は、前記事後分布推論手段が、前記Ｅステップにおいて、観測量の集合データと、隠れ変数の事後分布の推定値とについては、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込み、前記Ｍステップにおいて、観測量の集合データと、パラメータの事後分布の推定値とについては、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込むことを特徴とする。
【００２６】
かかる手順によれば、クラスタリング計算方法において、ＥステップおよびＭステップでは、当該推定対象の時刻から事前に定められた過去の一時点まで遡った観測量の集合データと、その演算に必要な推定値とを読みこむ。したがって、事前に過去の一時点を定めることなく、時間ステップの進展とともに推定すべき変数の個数が増加していくだけである場合と比べて、処理負荷の低減や処理の高速化を実現できる。
【００２７】
また、本発明に係るクラスタリング計算方法は、前記事後分布推論手段が、前記Ｅステップにおいて、推定対象の時刻から事前に設定された設定値分だけ遡及した過去の時刻までを演算対象時刻として、前記クラスタと、前記混合比と、前記重みとに関するパラメータの事後分布の値を演算し、前記Ｍステップにおいて、推定対象の時刻から事前に設定された設定値分まで遡及した過去の時刻までを演算対象時刻として、前記第１隠れ変数および第２隠れ変数の事後分布の値を演算することを特徴とする。
【００２８】
かかる手順によれば、クラスタリング計算方法において、ＥステップおよびＭステップでは、当該推定対象の時刻から事前に定められた設定値分だけ過去に遡った過去の時刻までを演算対象時刻として、推定値を求める演算を行う。したがって、事前に設定値を定めることなく、時間ステップの進展とともに推定すべき変数の個数が増加していくだけである場合と比べて、処理負荷の低減や処理の高速化を実現できる。
【００２９】
また、本発明に係るクラスタリング計算方法は、前記事後分布推論手段が、前記Ｅステップにおいて、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、事前に設定された再推定すべきクラスタの判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合にだけ、前記混合比に関するパラメータの事後分布の値を演算し、前記Ｍステップにおいて、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、前記判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合にだけ、前記第１隠れ変数の事後分布の値を演算することを特徴とする。
【００３０】
かかる手順によれば、クラスタリング計算方法において、ＥステップおよびＭステップでは、クラスタごとの演算のためにクラスタを更新したときに、必要な場合にだけそのクラスタにおける推定値の再推定処理を行う。したがって、ｄＨＤＰモデルにおいて予め設定される最大クラスタ数の個数のクラスタに関して、推定値の再推定処理を毎回実行する場合に比べて、処理負荷の低減や処理の高速化を実現できる。
【００３１】
また、本発明に係るクラスタリング計算プログラムは、前記いずれかのクラスタリング計算装置を構成する各手段としてコンピュータを機能させるためのプログラムである。このように構成されることにより、このプログラムをインストールされたコンピュータは、このプログラムに基づいた各機能を実現することができる。
【００３２】
また、本発明に係るコンピュータ読み取り可能な記録媒体は、前記クラスタリング計算プログラムが記録されたことを特徴とする。このように構成されることにより、この記録媒体を装着されたコンピュータは、この記録媒体に記録されたプログラムに基づいた各機能を実現することができる。
【発明の効果】
【００３３】
本発明によれば、ダイアライゼーションにおける話者クラスタリングの問題に対して、ノンパラメトリックベイズモデルを採用して、確率的なクラスタリングを用いることで、従来のようなパラメータの設定や探索によらずに容易に話者数を推定できる。また、本発明によれば、ノンパラメトリックベイズモデルを採用したので、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化できる。その結果、ダイアライゼーションにおいて、話者数と各話者を特徴づけるパラメータを正確に推定できる。
【図面の簡単な説明】
【００３４】
【図１】本発明の実施形態に係るクラスタリング計算装置を含むダイアライゼーションシステムの概要を示す構成図である。
【図２】本発明の実施形態に係るクラスタリング計算方法の全体処理の流れを示すフローチャートである。
【図３】図２に示す変分事後分布推論手順を示すフローチャートである。
【図４】図３に示すＥステップの計算手順の一例を示すフローチャートである。
【図５】図３に示すＭステップの計算手順の一例を示すフローチャートである。
【図６】本発明の実施形態に係るクラスタリング計算方法におけるＥステップの計算手順を示すフローチャートである。
【図７】本発明の実施形態に係るクラスタリング計算方法におけるＭステップの計算手順を示すフローチャートである。
【図８】本発明の実施形態に係るクラスタリング計算装置の構成の一例を示す機能ブロック図である。
【図９】本発明の実施形態に係るクラスタリング計算装置のクラスタリング性能を評価するために用いた人工音声データの時間平均パワー分布を示すグラフである。
【図１０】本発明の実施形態に係るクラスタリング計算装置のクラスタリング性能を評価するために用いた実音声データの時間平均パワー分布を示すグラフであって、（ａ）はＣＰ１、（ｂ）はＣＰ２、（ｃ）はＤＣ、（ｄ）はＣＮをそれぞれ示している。
【図１１】本発明の実施形態に係るクラスタリング計算装置によって人工音声データをクラスタリングした結果を示すグラフであって、（ａ）はＤＰＭ、（ｂ）はｄＨＤＰをそれぞれ示している。
【図１２】本発明の実施形態に係るクラスタリング計算装置によってＣＰ１の実音声データをクラスタリングした結果を示すグラフであって、（ａ）はＤＰＭ、（ｂ）はｄＨＤＰをそれぞれ示している。
【図１３】本発明の実施形態に係るクラスタリング計算装置によってＤＣの実音声データをクラスタリングした結果を示すグラフであって、（ａ）はＤＰＭ、（ｂ）はｄＨＤＰをそれぞれ示している。
【発明を実施するための形態】
【００３５】
図面を参照して本発明のクラスタリング計算装置およびクラスタリング計算方法を実施するための形態（以下「実施形態」という）について詳細に説明する。以下では、推論原理の概略、ダイアライゼーションシステムの概略、クラスタリング計算方法の概要、計算アルゴリズム、クラスタリング計算装置について順次説明する。
【００３６】
［推論原理の概略］
本実施形態では、ノンパラメトリックベイズモデルの一種として、例えば、ｄＨＤＰを用いることとする。ここでは、ｄＨＤＰを簡単に説明する。ｄＨＤＰの数学的なモデルを式（１）〜式（６）に示す。なお、ＤＰＭ等の他のノンパラメトリックベイズモデルを用いてもよいことはもちろんである。ＤＰＭの詳細は、例えば、「上田修功、山田武士、“ノンパラメトリックベイズモデル”，応用数理，Vol. 17， No. 3， pp. 196-214, 2007.」に記載されている。
【００３７】
【数１】
【００３８】
ここで、〜は確率分布からのサンプリングを表す。また、ＤＰ（・）はDirichlet Process（無限次元分布）を表し、γ，α０，ａ０，ｂ０は事前に設定するハイパーパラメータである。
【００３９】
ｄＨＤＰでは、まず、式（１）で離散無限個のパラメータ分布（クラスタ）Ｇ０を生成する。ダイアライゼーションにおいて、この分布Ｇ０は、データ全体を見たときの各発話者の発話割合と構成に相当する。ＤＰＭに代表されるノンパラメトリックベイズモデルでは、推定された解が自動的に少数のパラメータ（話者）からなる分布に集約される。
【００４０】
式（２）のＨtは、時刻ｔでの話者分布の時間変化（分布変化）を表す。
式（３）のｗtは、時刻ｔでの話者分布の時間変化の割合（程度）を表す。
式（５）のθt iは、時刻ｔでのクラスタのパラメータを表し、式（６）のｘt iは、時刻ｔでのサンプル分布を表す。なお、ｉは、時刻ｔでのｉ番目のデータを示す。
【００４１】
式（４）におけるＧtは、時刻ｔにおけるパラメータ分布である。ダイアライゼーションにおいて、分布Ｇtは、時刻ｔに実際に発話した発話者の構成を表現している。式（４）は、このＧtを、時刻ｔ−１での分布Ｇt-1と、時刻ｔでの分布変化を表すＨtとをｗtで重ね合わせた分布として表現したものである。よって、ｄＨＤＰでは、各時刻でのサンプル分布（ｘt iの分布）に変動が許される。一方で、時刻に不変な分布であるＧ０も推定するため、全体を通した話者クラスタも学習されている。
【００４２】
さらに、ｄＨＤＰでは、変化の割合であるｗtも動的に学習されるので、話者が交代するときには、劇的に分布が代わり、そうでないときには、ほとんど分布が変化しないというように、変化の割合が一定でないデータのモデル化にも対応できる。このように、ｄＨＤＰでは、時刻に依存して分布が少しずつ変化するというモデルとなっているため、話者の入れ替わり（turn-taking）によって観測量（サンプル）の分布が時間変化するダイアライゼーションのタスクにとって恰好の確率モデルとなっている。
【００４３】
［ダイアライゼーションシステムの概略］
図１は、本発明の実施形態に係るクラスタリング計算装置を含むダイアライゼーションシステムの概要を示す構成図である。本実施形態のダイアライゼーションは、前記した第２の方法（話者位置に関する情報を利用する方法）によるものとして説明する。予め、未知数話者による会話を録音し、ダイアライゼーションシステム１への入力とする。ここでは、室内１０１に３人の話者ＨA、ＨB、ＨCが図１のような定位置で会話したものとする。音声データ（会話の録音データ）１０２は時系列データとなる。
【００４４】
ダイアライゼーションシステム１は、１つの大型コンピュータ、または、複数のコンピュータからなる。ここでは、ダイアライゼーションシステム１は、３つのコンピュータ、すなわち、特徴量抽出部２と、クラスタリング計算装置３と、識別部４とを備える。
【００４５】
特徴量抽出部２は、ノイズ除去のような前処理を行うと共に、ダイアライゼーションに適した各種特徴量を抽出するものである。特徴量抽出部２は、例えば、マイクロホンアレーから取得した録音データから、ＤＯＡ（direction of arrival：音声到来方向）情報を抽出し、クラスタリング計算装置３に出力する。ＤＯＡ情報（音声の到達角度）は、マイクに対してどの方向からどの程度の強さの音声信号が観測されたかを推定した量になっている。
【００４６】
クラスタリング計算装置３は、ＤＯＡ情報に基づいてクラスタリング処理によって話者数と各話者を特徴づけるパラメータとを推定する。すなわち、クラスタリング計算装置３は、抽出した音声特徴量をクラスタリングし、クラスタ数と各クラスタのパラメータを推定する。クラスタリング計算装置３は、確率的なクラスタリングモデルを適用することで、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化する。
【００４７】
識別部４は、クラスタリング計算装置３で得られたクラスタリング結果（クラスタリング推定値）により、各時刻における話者の発話状態を識別する。この識別部４は、クラスタリング推定値を解析して、クラスタ数と、その位置とを、ユーザが識別できるような画面表示で提示する。
【００４８】
図１に例示したダイアライゼーション結果表示１０３のグラフでは、横軸が時刻（秒）、縦軸が方向（話者の位置）をそれぞれ示している。この例では、３人の話者ＨA、ＨB、ＨCに対応して、３つの方向に対して矩形波形が表示されている。各矩形波形の山となっている部分が各話者の発話（発言）を表している。まず話者ＨCが発言を終えると、次に話者ＨAが発言し、その発言途中に話者ＨCが再び発言し始めると、話者ＨAが黙り、話者ＨBが発言を開始する、というように、話者の入れ替わり（turn-taking）が生じていることがわかる。なお、２人の話者が同時に発話するタイミングも生じていることがわかる。
【００４９】
［クラスタリング計算方法の概要］
ここでは、推論原理のダイアライゼーションへの導入と、クラスタリングモデルとについて説明する。
＜推論原理のダイアライゼーションへの導入＞
ここでは、クラスタリング計算装置３において、ノンパラメトリックベイズモデル（ｄＨＤＰ）を用いるため、その前段の特徴量抽出部２で抽出する特徴量を定式化する。時刻ｔにおいて、角度ｄ（例えば、ｄ＝−１８０，−１７９，…，０，…，１８０）方向から聞こえてきた音声パワー（ＤＯＡ情報）をｆt dとする。すなわち、各時刻ｔにおける音声パワーベクトルは
である。なお、ｄ＝−１８０の方向と、ｄ＝１８０の方向とは同じものである。
【００５０】
このパワーベクトルを、本実施形態で用いるクラスタリングモデルに適合するように、１次元ベクトルｘt iの集合へと変換する。ここでは、各ｆt dについて、閾値定数τとパラメータμとを用いて、値ｇ（ｄ）をもつｎt d個のサンプル集合を生成することとした。ここで、関数ｇ（・）は、実装の都合に合わせて選択された適切なスケール関数を示す。例えば、関数ｇ（・）は、［−１８０：１８０］→［−０．５：０．５］に変換する関数を用いることができる。また、ｎt dは、式（７）で定義する。なお、１次元に限らず、例えば、３次元ベクトルｘt iの集合へ変換するように構成してもよい。
【００５１】
【数２】
【００５２】
式（７）の関数ｈ（・）は、パワーの値に応じて何らかの正整数を返す関数であり、定数でもよい。本実施形態では、例えば、ｈ＝１を用いる。以上の量子化変換を全ての角度ｄに関して行い、時刻ｔの観測量をサンプル集合に集約する。
【００５３】
【数３】
【００５４】
本実施形態で扱うクラスタリング問題は、サンプル集合Ｘt＝｛ｘ１，…，ｘt｝のクラスタリングと捉えることができる。人間の発話音声パワーは背景ノイズよりもはるかに大きいものとすれば、話者の位置を反映した値のサンプルが大量に観測されるため、主要なクラスタとして話者の位置が推定できるはずである。そこで、本実施形態のクラスタリング計算方法では、クラスタ数Ｋを推定しながら、同時に各サンプルｘt iに対してどのクラスタに属するかを表す変数ｚt iの確率分布を求めること、さらにＫ個のクラスタに対応するパラメータΘ＝｛θk｝を求めることをそれぞれ行う。
【００５５】
＜クラスタリングモデル＞
本実施形態では、計算量とアルゴリズムの簡易さとを考慮して、ｄＨＤＰ近似モデルを用いる。ｄＨＤＰ近似モデルの生成モデルは、以下の式（９）〜式（１５）のようになる。なお、ｄＨＤＰ近似モデルについては、「I. Pruteanu-Malinici, L. Ren, J. Paisley, E. Wang and L. Carin, “Dynamic Hierarchical Dirichlet Process for Modeling Topics in Time-Stamped Documents”, IEEE Transactions on Pattern Analysis and Machine Intelligence, submitted, 2008.」に記載されている。
【００５６】
【数４】
【００５７】
ｄＨＤＰ近似モデルでは、まず、最大クラスタ数Ｋを固定しておく。最大クラスタ数Ｋは、推定しようとしている話者数よりも充分大きい値（例えば、数十〜１００）であればよい。実際には、「有効な」クラスタ数Ｋeff（＜Ｋ）なる個数のクラスタが、推定しようとしている話者に対応することとなる。「有効な」クラスタ数Ｋeffについては、学習の結果、話者に対応しないそれ以外のクラスタの重み（混合比）が自動的にほぼ０になることを利用して判断することができる。すなわち「有効な」クラスタ数Ｋeffは、推定の過程で自動的に決定される。
【００５８】
≪式（９）、式（１０）≫
ｄＨＤＰ近似モデルでは、式（９）で、Ｋ個のクラスタに対応するパラメータをサンプリングする。式（１０）では、有限次元のDirichlet分布よりinnovation measure Ｈt（後記する式（１６ｃ）参照）のサンプリングを行う。より具体的には混合比πtをサンプリングする。これは、前記した式（４）より、式（１６ａ）およびその変形式である式（１６ｂ）が導かれるため、時刻ｔにおける話者分布Ｇtは、時刻ｔまでに生成したＨ1:lの重ね合わせだけで表現できるからである。ここで、「Ｈ1:l」は、Ｈ1〜Ｈlを表す。
【００５９】
【数５】
【００６０】
≪式（１１）、式（１２）≫
ｄＨＤＰ近似モデルでは、続いて、式（１１）により、話者分布Ｇtの時間変化の程度を表すｗtをサンプリングして、このｗtを用いて、式（１２）で定義されるｖtl（ｌ＝１，…，ｔ）を計算する。ここで、時刻ｌは、時刻ｔとそれよりも過去の時刻とを表す。
【００６１】
≪式（１３）≫
式（１３）に示す隠れ変数ｄt iは、ｔ次元の｛０，１｝ベクトルである（ｔは時刻、要素の値は０か１のみ）。隠れ変数ｄt iは、時刻ｔのときのｔ次元の要素のうち、ｌ次元目の要素ｄt i lが１の値をもつものである（ｌ≦ｔ）。隠れ変数ｄt iのｌ次元目の要素ｄt i lは、時刻ｌでの分布変化Ｈlに対応する要素である。時刻ｌでの分布変化Ｈlは、時刻ｔのｉ番目のサンプルｘt iをサンプルするための分布変化である。
【００６２】
≪式（１４）≫
式（１４）に示す隠れ変数ｚt iは、同様に、Ｋ次元の｛０，１｝ベクトルである。
隠れ変数ｚt iは、実際にサンプルｘt iをサンプリングするクラスタ（パラメータ）ｋに対応するｋ次元目の要素ｚt i kにのみ１の値を持つ。
【００６３】
≪式（１５）≫
式（１５）に示す観測量ｘt iは、与えられたクラスタ番号（ｋ）に対応するパラメータθから生成される。式（１５）は、前記した式（６）、式（８）と同様であって、別の表式である。
【００６４】
≪観測分布と事前分布≫
ｄＨＤＰ近似モデルでは、式（１５）に示す観測分布Ｆと、パラメータの事前分布Ｈ（式（９）参照）とを事前に定めておく必要がある。本実施形態では、一例として、観測分布Ｆが正規分布であり、かつ、パラメータの事前分布Ｈが共役事前分布であるようなＮｏｒｍａｌ−Ｇａｍｍａ分布を用いることにした。なお、Ｎｏｒｍａｌ−Ｇａｍｍａ分布については、参考文献１「C. M. ビショップ、“パターン認識と機械学習”、シュプリンガー・ジャパン、2007.」に記載されている。
【００６５】
本実施形態でＮｏｒｍａｌ−Ｇａｍｍａ分布を採用した理由は、これらの分布を利用したｄＨＤＰの近似モデル解法は、従来発表されていない上、多くの分野に応用可能なモデルであり、最も実用性が高いと考えたからである。ただし、目的や実際のデータに則して、これらを別の分布にすることも可能である。
【００６６】
≪ｄＨＤＰ近似モデルのダイアライゼーションシステムへの適用≫
図１のダイアライゼーションシステム１では、式（９）〜式（１５）に示すｄＨＤＰ近似モデルにおいて、ダイアライゼーションの話者数（＝クラスタ数）と各サンプルのクラスタリングを行い、各クラスタのパラメータ（＝話者の位置）を推定することとした。
【００６７】
このうち、各サンプルのクラスタリングは、式（１４）に示す隠れ変数ｚt iの分布ｐ（ｚt i）を求めることに等しい。隠れ変数ｚt iは、各サンプルｘt iに対してどのクラスタに属するかを表す変数なので、このクラスタリング結果が得られれば、各クラスタへ帰属するサンプル数（あるいはその期待値）が計算できる。また、この各クラスタへ帰属するサンプル数（あるいはその期待値）を用いて各クラスタｋの混合比（後記するβk＾）を計算可能である。さらに、この混合比（後記するβk＾）により、「有効な」クラスタ数Ｋeff（＝話者数Ｋeff）を決定することができる。また、クラスタリング結果が求まれば、各クラスタのパラメータ｛θk｝も容易に求めることができる。なお、本実施形態では、クラスタリング計算装置３によって、「有効な」クラスタ数Ｋeff（＝話者数Ｋeff）の決定まで行うこととするが、識別部４がこの処理を行ってもよい。つまり、クラスタリング計算装置３によって、クラスタリング結果を得て、各クラスタのパラメータ｛θk｝を得て、識別部４によって、混合比（後記するβk＾）を求め、「有効な」クラスタ数Ｋeff（＝話者数Ｋeff）を決定するようにしてもよい。
【００６８】
［計算アルゴリズム］
計算アルゴリズムについては、以下、１）ｄＨＤＰモデルのオンライン推定法、２）変分事後分布推論プロセス、３）観測モデルおよび事前分布の特定、４）推定結果とクラスタ数の決定方法、５）ｄＨＤＰの高速化方法に分けて詳細に説明する。
【００６９】
＜１）ｄＨＤＰモデルのオンライン推定法＞
ここでは、具体的な推論アルゴリズムを示す。図２は、本発明の実施形態に係るクラスタリング計算方法の全体処理の流れを示すフローチャートである。図２は、クラスタリングの全体の推論プロセスを示したものである。ダイアライゼーションは、一般的にオンライン計算のプロセスである。しかしながら、ｄＨＤＰモデルのオンライン推定方法は、従来研究されていない。本実施形態では、ｄＨＤＰモデルのオンライン推定法を開発した。図２の符号２０１で示す推論プロセスは、そのオンライン推定法を示したものである。
【００７０】
オンライン推定法は、時刻Ｔ（１≦Ｔ≦Ｔtotal）ごとに行うものである。ここで、刻々と進展する推定対象のある時点の時刻をＴ、それらのうち最終時刻をＴtotalとした。なお、以下では、過去を含めた演算対象時刻をｔとする。過去にどこまで遡及するかについては予め定めておく。例えば、推定対象時刻Ｔ＝５のとき、演算対象時刻ｔ＝１，２，３，４，５としたり、ｔ＝３，４，５としたりすることができる。この例では、ｔ＝１まで考慮することとする。
【００７１】
この推論プロセスでは、処理を開始すると、まず、Ｔを初期化する。すなわち、Ｔ＝１とする（ステップＳ１）。そして、時刻１からＴまでの隠れ変数およびパラメータ（隠れ変数・パラメータ）の推定値およびハイパーパラメータを入力する（ステップＳ２）。また、時刻１からＴ−１までの観測量（サンプル）を入力する（ステップＳ３）。また、時刻Ｔの音声パワーｆTを入力する（ステップＳ４）。この時刻Ｔの音声パワーｆTから、時刻Ｔの観測量｛ｘT i｝を生成して入力する（ステップＳ５）。また、時刻Ｔに対応する未知数（未知の隠れ変数・パラメータ）を初期化する（ステップＳ６）。
【００７２】
そして、隠れ変数・パラメータの変分事後分布を推定する（ステップＳ７）。推定後、Ｔをインクリメントする。すなわち、Ｔ←Ｔ＋１とする（ステップＳ８）。そして、入力が終了したか否かを判別する（ステップＳ９）。入力が終了していない場合（ステップＳ９：Ｎｏ）、ステップＳ２に戻る。一方、例えば、最終時刻Ｔtotalを超えたときに、入力が終了したと判定し（ステップＳ９：Ｙｅｓ）、推定結果を出力する（ステップＳ１０）。なお、ステップＳ２〜ステップＳ４の処理順序は、任意であり、並列処理してもよい。
また、ステップＳ６の処理は、ステップＳ７の前に行うのであれば処理順序は問わない。
【００７３】
オンライン推定法を繰り返すことで、時刻ステップの進展と共に観測サンプル｛ｘt i｝が蓄積され、その都度変数を再推定する。時刻Ｔ−１の時点での推定結果は、次の時刻Ｔでの推定の初期値として利用される。
【００７４】
推論の目的は、全観測データが与えられた時に全未知変数（｛ｚt i｝，｛ｄt i｝，ｗ，｛πt｝，｛θk｝）の推定値を求めることである。これは、確率モデルの観点からは、全変数の事後分布を求めることに相当する。本実施形態では、ｄＨＤＰ近似モデルに対して変分ベイズ法による事後分布推定法（変分事後分布推定法）を示す。実際に求めたい事後分布は、式（１７）であるが、変分法では、式（１８）のように全変数が独立であると仮定した分布（変分事後分布）を推定する。
【００７５】
【数６】
【００７６】
式（１８）に示すｑ（・）は変分事後分布を表す。ここで、時刻Ｔまでの観測量が得られたとする。このとき、現在の時刻Ｔからの過去を含む時刻ｔ（≦Ｔ）に関する変数の変分事後分布推定値ｑ＊は、式（１９）〜式（２３）のように表される。これら各変数の分布は、元の分布（式（９），式（１０），式（１１），式（１３），式（１４））からデータの情報が加わった分修正された形になる。ただし、式（１９）は、観測モデルＦやパラメータの事前分布Ｈに依存する。したがって、本実施形態では、式（１９）に基づいて後記する式（３６）を用いる。
【００７７】
【数７】
【数８】
【００７８】
ただし、式（１９）〜式（２３）を同時に最適化することは困難なので、変分法においてはＥＭアルゴリズム（Expectation−Maximization algorithm）という様式に従った繰り返し計算による推定で各変数を個別に最適化する。ＥＭアルゴリズムとは、複数の変数を同時最大化する計算手法であり、Ｅステップ（Expectation step）とＭステップ（Maximization step）からなる計算ステップを相互に繰り返し計算することで全体最適化を行うものである。なお、ＥＭアルゴリズムについては、前記した参考文献１に記載されている。
【００７９】
＜２）変分事後分布推論プロセス＞
図３は、図２に示す変分事後分布推論手順を示すフローチャートであって、ＥＭアルゴリズムを含んだ推論プロセスを示す。図３に示す推論プロセスでは、まず、時刻１からＴまでのサンプル｛ｘ1:T｝を入力する（ステップＳ２１）。また、時刻１からＴまでの隠れ変数・パラメータの推定値およびハイパーパラメータを入力する（ステップＳ２２）。なお、ステップＳ２１，Ｓ２２の処理順序は任意である。
【００８０】
そして、ＥＭステップの繰り返し回数を表す識別子ｊを初期化する。すなわち、ｊ＝１とする（ステップＳ２３）。そして、Ｅステップを計算する（ステップＳ２４）。Ｅステップでは、時刻１からＴまでのパラメータの推定値を更新する。続いて、Ｍステップを計算する（ステップＳ２５）Ｍステップでは、時刻１からＴまでの隠れ変数の推定値を更新する。
【００８１】
そして、ＥＭステップの繰り返し回数ｊをインクリメントする。すなわち、ｊ←ｊ＋１とする（ステップＳ２６）。そして、現在の繰り返し回数ｊが事前に設定されたしきい値（ｊmax）を超えたか否かを判別する。すなわち、ｊ＞ｊmaxとなったか否かを判別する（ステップＳ２７）。ｊ≦ｊmaxである場合（ステップＳ２７：Ｎｏ）、ステップＳ２４に戻る。一方、ｊ＞ｊmaxとなった場合（ステップＳ２７：Ｙｅｓ）、時刻Ｔが事前に設定された設定値（適切な正整数ｔupdateの倍数）を超えたか否かを判別する。すなわち、Ｔが「ｔupdateの倍数」となったか否かを判別する（ステップＳ２８）。
【００８２】
Ｔが「ｔupdateの倍数」となった場合（ステップＳ２８：Ｙｅｓ）、ハイパーパラメータを更新し（ステップＳ２９）、推定結果を保存する（ステップＳ３０）。Ｔが「ｔupdateの倍数」ではない場合（ステップＳ２８：Ｎｏ）、ステップＳ２９をスキップして、推定結果を保存する（ステップＳ３０）。
【００８３】
このようにＥステップ、Ｍステップをｊmax回に渡って相互に繰り返し計算することで変分事後分布の式（１９）〜式（２３）を得ることができる。なお、これらは、Ｍステップにおいて演算される。
【００８４】
本実施形態では、図３に示すようにステップＳ２８〜Ｓ３０の処理を行うことで、時間ステップＴを適切な正整数ｔupdate回インクリメントするごとにハイパーパラメータを推定することとした。通常、ハイパーパラメータは事前に与える定数であるが、ステップＳ２８〜Ｓ３０の処理を実行することで、ハイパーパラメータが固定値である場合よりも精度を上げることができる。ハイパーパラメータの更新方法については、事後分布からのサンプリング方法など、様々な公知の手法を用いることができる。また、ハイパーパラメータ自身もＥＭアルゴリズムで推定することも可能である。
【００８５】
≪Ｅステップ≫
Ｅステップの具体的な計算式は、以下の式（２４）〜式（２８）で表される。ここで、ψ（・）はプサイ関数（あるいはディガンマ関数）である。また、式（２７）は、観測モデルＦに依存する。したがって、本実施形態では、式（２７）に基づいて後記する式（３７）を用いる。
【００８６】
【数９】
【００８７】
Ｅステップでは、未知パラメータ（θ、π、ｗ）に関わる推定値を再計算する。ここでＥx［ｆ（ｘ）］は、変分分布上の期待値を示す。式（２８）にその定義を示す。
【００８８】
【数１０】
【００８９】
図４は、図３に示すＥステップの計算手順の一例を示すフローチャートである。Ｅステップでは、時刻Ｔまでに関係する全変数について再計算を行う。具体的には、Ｅステップでは、まず、時刻１からＴまでのサンプル｛ｘ1:T｝を入力する（ステップＳ３１）。また、時刻１からＴまでの隠れ変数・パラメータの推定値とハイパーパラメータ、最新のＭステップの演算結果を入力する（ステップＳ３２）。なお、ステップＳ３１，Ｓ３２の処理順序は任意である。
【００９０】
そして、過去を含む演算対象の時刻ｔを初期化する。すなわち、ｔ＝１とする（ステップＳ３３）。そして、この演算対象の時刻ｔに対して、式（２４）、式（２５）を計算する（ステップＳ３４）。次いで、クラスタの識別子ｋを初期化する。すなわち、ｋ＝１とする（ステップＳ３５）。そして、この演算対象の時刻ｔおよびクラスタｋに対して、式（２６）を計算する（ステップＳ３６）。さらに、時刻ｔにおけるデータの識別子ｉを初期化する。すなわち、ｉ＝１とする（ステップＳ３７）。そして、時刻ｔ、クラスタｋのｉ番目のデータに対して、式（２７）を計算する（ステップＳ３８）。
【００９１】
そして、時刻ｔにおけるデータの識別子ｉをインクリメントする。すなわち、ｉ←ｉ＋１とする（ステップＳ３９）。続いて、ｉ＞ｎtとなったか否かを判別する（ステップＳ４０）。なお、ｎtは、式（８）に示した個数である。ｉ≦ｎtである場合（ステップＳ４０：Ｎｏ）、ステップＳ３８に戻る。一方、ｉ＞ｎtとなった場合（ステップＳ４０：Ｙｅｓ）、次のクラスタに更新する。すなわち、ｋ←ｋ＋１とする（ステップＳ４１）。そして、すべてのクラスタについて演算したか否かを判別する。すなわち、ｋ＞Ｋとなったか否かを判別する（ステップＳ４２）。ｋ≦Ｋである場合（ステップＳ４２：Ｎｏ）、ステップＳ３６に戻る。一方、ｋ＞Ｋとなった場合（ステップＳ４２：Ｙｅｓ）、演算対象の時刻ｔを更新する。すなわち、ｔ←ｔ＋１とする（ステップＳ４１）。そして、演算対象の時刻ｔが推定対象の時刻Ｔとなったか否かを判別する。すなわち、ｔ＞Ｔとなったか否かを判別する（ステップＳ４４）。ｔ≦Ｔである場合（ステップＳ４４：Ｎｏ）、ステップＳ３４に戻る。一方、ｔ＞Ｔとなった場合（ステップＳ４４：Ｙｅｓ）、推定対象の時刻ＴにおけるＥステップの推定結果を保存する（ステップＳ４５）。
【００９２】
≪Ｍステップ≫
Ｍテップの具体的な計算式は、以下の式（２９）〜式（３２）で表される。Ｍステップでは、隠れ変数（ｚt，ｄt）に関わる推定値を再計算する。
【００９３】
【数１１】
【００９４】
図５は、図３に示すＭステップの計算手順の一例を示すフローチャートである。Ｍステップでも、Ｅステップと同様に、時刻Ｔまでに関係する全変数について再計算を行う。具体的には、Ｍステップでは、まず、時刻１からＴまでのサンプル｛ｘ1:T｝を入力する（ステップＳ５１）。また、時刻１からＴまでの隠れ変数・パラメータの推定値とハイパーパラメータ、最新のＥステップの演算結果を入力する（ステップＳ５２）。なお、ステップＳ５１，Ｓ５２の処理順序は任意である。
【００９５】
以下、Ｅステップと同様な処理については、説明を適宜省略する。Ｍステップでは、まず、ｔ＝１（ステップＳ５３）、ｉ＝１（ステップＳ５４）、ｋ＝１（ステップＳ５５）とした上で、前記した式（２９）、式（３０）を計算する（ステップＳ５６）。その後、ｋ←ｋ＋１として（ステップＳ５７）、ｋ＞ＫとなるまでステップＳ５６に戻る処理を繰り返す。ｋ＞Ｋとなった場合（ステップＳ５８：Ｙｅｓ）、時刻の識別子ｌを初期化する。すなわち、ｌ＝１とする（ステップＳ５９）。そして、時刻ｔ、時刻ｌ、クラスタｋのｉ番目のデータに対して、式（３１）、式（３２）を計算する（ステップＳ６０）。
【００９６】
そして、時刻ｌをインクリメントする。すなわち、ｌ←ｌ＋１とする（ステップＳ６１）。過去から時刻ｔまでの時刻を示す時刻ｌについてのすべての演算が終わるまでステップＳ６０に戻る処理を繰り返す。ｌ＞ｔとなった場合（ステップＳ６２：Ｙｅｓ）、ｉ←ｉ＋１とする（ステップＳ６３）。そして、ｉについてのすべての演算が終わるまでステップＳ５５に戻る処理を繰り返す。そして、ｉ＞ｎtとなった場合（ステップＳ６４：Ｙｅｓ）、ｔ←ｔ＋１とする（ステップＳ６５）。さらに、ｔについてのすべての演算が終わるまでステップＳ５４に戻る処理を繰り返す。そして、ｔ＞Ｔとなった場合（ステップＳ６６：Ｙｅｓ）、推定結果を保存する（ステップＳ６７）。
【００９７】
＜３）観測モデルおよび事前分布の特定＞
図３〜図５および式（１９）〜式（２７）を参照して説明した変分事後分布推論プロセスでは、観測モデルおよび事前分布を一般化した説明を行った。つまり、前記した式（１９）および式（２７）は、観測モデルＦ＝ｐ（ｘt i|θk）およびパラメータθkの事前分布Ｈに依存する。本実施形態では、観測モデルＦを正規分布と仮定し、事前分布ＨとしてＮｏｒｍａｌ−Ｇａｍｍａ分布を仮定する。Ｎｏｒｍａｌ−Ｇａｍｍａ分布は、式（３３）〜式（３５）で表されるモデルである。
【００９８】
【数１２】
【００９９】
このモデルの場合において、前記した式（１９）および式（２７）は、式（３６）、式（３７）でそれぞれ表される。また、式（３６）、式（３７）中のハイパーパラメータは、式（３８）〜式（４１）で表される。さらに、式（３８）〜式（４１）中の変数は、式（４２ａ）〜式（４２ｃ）で表される。式（３８）〜式（４２）の演算は、Ｅステップの中で完了する。
【０１００】
【数１３】
【０１０１】
【数１４】
【０１０２】
【数１５】
【０１０３】
＜４）推定結果とクラスタ数の決定方法＞
≪クラスタ数の決定方法≫
ＥＭアルゴリズムでは、Ｋ個のクラスタを常に保持しているが、推定が進むと少数のクラスタのみに大きな混合比が付き、他のクラスタの大きさ（cluster size）はほぼ０となる。時刻ｔにおいてクラスタｋへ振り分けられるサンプル数の期待値は、式（４３）の定義により計算できる。
【０１０４】
【数１６】
【０１０５】
本実施形態では、「有効な」クラスタ数Ｋeffを、式（４３）で定義した‖ｚt,k‖の時刻ｔに関する総和の比率で決定する。例えば、各クラスタｋの混合比は、式（４４ａ）で推定できる。このルールでは、式（４４ｂ）の条件が成立するならば、そのクラスタｋが「有効な」クラスタであると判断することとする。
【０１０６】
【数１７】
【０１０７】
このルールにおいて、そのようなクラスタｋの総数を、「有効な」クラスタ数Ｋeffとすることとする。このようにすれば、「有効な」クラスタ数が最大クラスタ数Ｋより小さくなること、すなわち、Ｋeff＜Ｋが保証される。
【０１０８】
≪保存すべき推定結果≫
保存すべき推定結果としては、第１に、ＥＭアルゴリズムで推定された各変数によって定められる変分事後分布である。また、第２に、クラスタリング結果を用いて求められる、「有効な」クラスタ数Ｋeff、式（４３）に示す‖ｚt,k‖、式（４４ａ）に示すβk＾などである。ここで、＾は、文字βの上に付される記号を意味する。
特に、保存すべき推定結果として２番目に挙げたクラスタリング結果を用いて求められる推定量“「有効な」クラスタ数Ｋeff、式（４３）に示す‖ｚt,k‖、式（４４ａ）に示すβk＾”は、図１における識別部４で利用することが可能な重要な推定量である。なお、式（４３）、式（４４ａ）、式（４４ｂ）の演算は、Ｍステップの中で行う。
【０１０９】
＜５）ｄＨＤＰの高速化方法＞
本実施形態のクラスタリング計算方法におけるオンライン推定方法では、図４や図５の処理の流れからも明らかなように、時間ステップの進展とともに推定すべき変数の数が増えていくことが特徴である。そのため、実時間性などを考慮して、計算省力化の技法を開発した。本実施形態では、ｄＨＤＰを用いたオンライン話者クラスタリングのための計算高速化法として、大別して３種類の高速化法（高速化１、高速化２、高速化３）を用いることとした。
【０１１０】
≪高速化１：データの忘却≫
前記した式（１６ａ）が意味することは、時刻ｔにおける話者分布Ｇtを計算するためには、時刻１≦ｌ≦ｔの分布変化Ｈｌが必要であるということである。従って、時刻１からの情報を保持し続ける必要があるので、時間ステップｔが進むとともに推論の計算量が増大する。ここで、次の仮定を導入する。すなわち、時刻ｌ＜ｔで話者の入れ替わりが発生したとする。すると、この時点で話者分布Ｈが大きく変化することが予想される。これは、前記した式（４）において、ｗl≒１となってＧl-1の影響がほぼ排除されることを意味する。このことと、前記した式（１６ａ）および式（１６ｂ）とを照らし合わせると、ｖt1≒…≒ｖl(t-1)＝０となる。従って、実際には、Ｇtの推論には、時刻ｌの分布変化Ｈl以降の分布のみしか関与しないことがわかる。このことから、分布選択を表す変数｛ｄt i｝の事後分布に相当するｓt i l（式（３０）参照）も大部分が０となることが期待される。
【０１１１】
ＥＭアルゴリズム中でｓが存在する部分は、ｓと、定数あるいは別変数と、の掛け算となっているため、ｓt i l＝０となる部分については、計算する必要がない。従って、ｓに対応する変数にはアクセスしなくてよい。この考察より、ある適切な時間ステップ幅Ｗ１より過去の変数あるいは定数については、ＥＭ更新式中でアクセスしない（忘却する）ことによって計算時間を削減できる。
【０１１２】
この時間ステップ幅Ｗ１の決定方法としては、事前に適切な定数を決めておく方法が最も簡便である。あるいは、この考察からの示唆より、次のような方法も考えることができる。すなわち、各時刻でのｓt iの推定結果より、ｓt i 1＝…＝ｓt i l＝０となるｌを見つけて、Ｗ１=ｔ−ｌとする方法がある。そのようなｌの見付け方としては、適切な閾値ｔｈ（＜１．０）に対して、式（４５）の関係を満たす最大のｌを探せばよい。この場合には、データを忘却する時間幅Ｗ１が推論結果に従って動的に変化することになる。
【０１１３】
【数１８】
【０１１４】
≪高速化２：推定時間幅の制限≫
図２に示したオンライン推定法では、毎時刻Ｔにおいて全変数のＥＭアルゴリズム推定を行うものとして説明した。このことは、早い時刻ステップに関する変数については、何度もＥＭ再推定を行っているということを意味する。ＥＭアルゴリズムは、局所最適解の計算アルゴリズムであるため、計算を繰り返すと、ある１つの解に収束する。よって、早い時刻ステップに関わる変数については再推定をしなくとも、値が収束している可能性が高い。
【０１１５】
この考察より、ＥＭアルゴリズムで再推定する変数を、適切な時間幅Ｗ２を使うことで、Ｔ−Ｗ２≦ｔ≦Ｔの範囲の変数に制限することが可能であることがわかる。このことで、推定する変数の個数が、時間ステップＴに対する線形増加とならないように、その個数を定数に保つことが可能となる。時間幅Ｗ２の決定方法としては、事前に適切な定数を決めておく方法が最も簡便である。
【０１１６】
≪高速化３：クラスタ変数再推定の制限≫
ｄＨＤＰは、ノンパラメトリックベイズモデルなので、最大クラスタ数（Ｋ）の個数のクラスタを予めモデルに準備していたとしても、実際には、「有効な」クラスタ数Ｋeff（＜Ｋ）の個数のクラスタだけが実質的なクラスタとなり、他のクラスタは、混合比がほぼ０となるために消去されたような形式となる。この消去されたような形式のクラスタには、有効な情報がないため、そのようなクラスタを用いてパラメータや混合比の推定計算を行っても無駄なだけである。
【０１１７】
このことから、前記した式（２６）、式（２７）、式（２９）などクラスタごとの推定を行う変数（ｋ）を用いて、実際にクラスタごとの変数再推定を行う演算処理では、その変数再推定の回数を確率的に（確率論に基づいて）低減することが考えられる。その方法は、例えば３種類挙げることができる。
【０１１８】
第１の低減法は、単純な方法として、全クラスタについて確率ｃ（≦１．０）で毎回ランダムに、再推定するか否かを決定する方法である。
第２の低減法は、各クラスタｋの混合比βk＾（式（４４）参照）に応じてＥＭアルゴリズムによる再推定を行う回数を増減する方法である。
第３の低減法は、前記した式（４４ｂ）の条件を満たす場合、つまり、対象のクラスタが「有効な」クラスタである場合には毎回再推定するが、それ以外の場合には、クラスタｋを確率的に更新する方法である。
これらの低減法において、各クラスタｋを更新する確率をｐupdate（ｋ）とおくと、第１〜第３の低減法は、それぞれ、式（４６）〜式（４８）のように表現できる。特に、第３の低減法は、推定精度をほとんど犠牲にせずに計算量をＫeff／Ｋ程度に削減できる効果がある。
【０１１９】
【数１９】
【０１２０】
前記した３つの高速化方法（高速化１、高速化２、高速化３）をすべて使う場合のＥＭアルゴリズムの一例を図６および図７にそれぞれ示す。なお、これら３つの方法は、それぞれ独立に利用可能である。
【０１２１】
図６は、本発明の実施形態に係るクラスタリング計算方法におけるＥステップの計算手順を示すフローチャートである。なお、図４のフローチャートと比較して異なる部分を太線及び破線で示し、説明を適宜省略して異なる処理について説明する。図６において、高速化１は、符号３０１で示すように、ステップＳ３１Ａ，Ｓ３２Ａからなる。ここでは、時刻Ｔ−Ｗ1からＴまでのサンプル｛ｘT-W1:T｝を入力する（ステップＳ３１Ａ）。また、時刻Ｔ−Ｗ1からＴまでの隠れ変数・パラメータの推定値とハイパーパラメータ、最新のＭステップの演算結果を入力する（ステップＳ３２Ａ）。次いで、高速化２は、符号３０２で示すように、ステップＳ３３Ａからなる。ここでは、ｔの初期値を１の代わりに、ｔ＝Ｔ−Ｗ2とする（ステップＳ３３Ａ）。
【０１２２】
さらに、高速化３は、ステップＳ３５に続いて行う処理であり、符号３０３で示すように、ステップＳ７１，Ｓ７２からなる。ここでは、［０，１］の一様乱数ｕを生成し（ステップＳ７１）、ｐupdate（ｋ）＞ｕであるか否かを判別することとした（ステップＳ７２）。ｐupdate（ｋ）＞ｕである場合（ステップＳ７２：Ｙｅｓ）、前記したステップＳ３６〜ステップＳ４５を実行する。つまり、前記した式（２６）および式（２７）の計算を行う。ただし、ステップＳ４２において、ｋ≦Ｋである場合（ステップＳ４２：Ｎｏ）、ステップＳ７１に戻る。一方、ステップＳ７２において、ｐupdate（ｋ）≦ｕである場合（ステップＳ７２：Ｎｏ）、前記したステップＳ３６〜ステップＳ４０をスキップして、ステップＳ４１に進む。つまり、前記した式（２６）および式（２７）の計算を行わない。
【０１２３】
図７は、本発明の実施形態に係るクラスタリング計算方法におけるＭステップの計算手順を示すフローチャートである。なお、図５のフローチャートと比較して異なる部分を太線及び破線で示し、説明を適宜省略して異なる処理について説明する。図７において、高速化１は、符号４０１で示すように、ステップＳ５１Ａ，Ｓ５２Ａからなる。ここでは、時刻Ｔ−Ｗ1からＴまでのサンプル｛ｘT-W1:T｝を入力する（ステップＳ５１Ａ）。また、時刻Ｔ−Ｗ1からＴまでの隠れ変数・パラメータの推定値とハイパーパラメータ、最新のＥステップの演算結果を入力する（ステップＳ５２Ａ）。次いで、高速化２は、符号４０２で示すように、ステップＳ５３Ａからなる。ここでは、ｔの初期値を１の代わりに、ｔ＝Ｔ−Ｗ2とする（ステップＳ５３Ａ）。
【０１２４】
さらに、高速化３は、ステップＳ５５に続いて行う処理であり、符号４０３で示すように、ステップＳ８１，Ｓ８２からなる。ここでは、［０，１］の一様乱数ｕを生成し（ステップＳ８１）、ｐupdate（ｋ）＞ｕであるか否かを判別することとした（ステップＳ８２）。ｐupdate（ｋ）＞ｕである場合（ステップＳ８２：Ｙｅｓ）、前記したステップＳ５６〜ステップＳ６７を実行する。つまり、前記した式（２９）〜式（３２）の計算を行う。ただし、ステップＳ５８において、ｋ≦Ｋである場合（ステップＳ５８：Ｎｏ）、ステップＳ８１に戻る。一方、ステップＳ８２において、ｐupdate（ｋ）≦ｕである場合（ステップＳ８２：Ｎｏ）、前記したステップＳ５６をスキップして、ステップＳ５７に進む。つまり、前記した式（２９）および式（３０）の計算を行わない。
【０１２５】
［クラスタリング計算装置］
図８は、本発明の実施形態に係るクラスタリング計算装置の構成の一例を示す機能ブロック図である。このクラスタリング計算装置３は、例えば、図２の推論プロセス等を実現するものであり、ＣＰＵ等の演算装置と、メモリ、ハードディスク等の記憶装置（記憶手段）と、マウスやキーボード等の外部から情報の入力を検出する入力装置と、外部との各種情報の送受信を行うインタフェース装置と、ＬＣＤ（Liquid Crystal Display）等の表示装置を備えたコンピュータと、このコンピュータにインストールされたプログラムとから構成される。
【０１２６】
クラスタリング計算装置３は、ハードウェア装置とソフトウェアとが協働することによって、前記したハードウェア資源がプログラムによって制御されることにより実現され、図８に示すように、記憶手段１０と、以下ＣＰＵの機能として、変分事後分布推論部３０と、パワーベクトル読込部２１と、入力制御部２２と、割当部２３と、未知数初期化部２４と、パワーベクトル書込部２５と、観測量生成部２６と、終了判定部２７と、出力制御部２８とを備えることとした。
【０１２７】
＜記憶手段＞
記憶手段１０は、ＲＯＭ、ＲＡＭ、ＨＤＤ等からなる。記憶手段１０は、プログラム記憶領域、設定データ記憶領域、演算データ一時記憶領域、推定結果記憶領域等に区分され、コマンド、データ、プログラム等各種情報を記憶する。例えば、推定結果記憶領域には、データとして、図９に示すように、未知数の事後分布推定値１１、ハイパーパラメータ１２、Ｅステップにおける推定値（まとめてＥステップ１３と表記した）、Ｍステップにおける推定値（まとめてＭステップ１４と表記した）、パワーベクトル１５、観測量（サンプル）１６が保存される。
【０１２８】
＜変分事後分布推論部＞
メインの変分事後分布推論部３０は、例えば、図３、図６、図７の推論プロセスを実現するものであり、ここでは、Ｅステップ用計算部３１と、Ｍステップ用計算部３２と、ＥＭ収束判定部３３と、パラメータ更新条件判定部３４と、ハイパーパラメータ更新部３５とを備えることとした。詳細は後記する。
【０１２９】
＜パワーベクトル読込部＞
パワーベクトル読込部２１は、推定対象の時刻Ｔにしたがって、パワーベクトルｆTを順次読み込み、入力制御部２２に渡す。
【０１３０】
＜入力制御部＞
入力制御部２２は、パワーベクトルｆTを取得すると、パワーベクトル書込部２５にそのまま送ると共に、割当部２３、未知数初期化部２４、観測量生成部２６および終了判定部２７に対しては、それぞれ必要なコマンドを出力する。
【０１３１】
＜割当部＞
割当部２３は、入力制御部２２からのコマンドにしたがって、そのときのＴ（＝１，２，…）や、角度ｄ（１８０，−１７９，…）に対応したｉに応じて、前記した式（９）〜（１４）の各未知数（パラメータ・隠れ変数）を順次生成し（初期値は例えば０）、記憶手段１０において、未知数の事後分布推定値１１として保存する。すなわち、割当部２３は、推定すべき未知数にＴとｉとを割り当てる。この意味で、図９においては、Ｔとｉの割当部２３と表記した。なお、記憶手段１０には、ｄＨＤＰ近似モデルで用いる各ハイパーパラメータ１２の初期値が予め格納されている。
【０１３２】
＜未知数初期化部＞
未知数初期化部２４は、入力制御部２２からのコマンドにしたがって、そのときのＴや角度ｄに対応したｉに応じて、Ｅステップ演算に用いるための、前記した式（２４）〜（２６）、（３７）の左辺パラメータを順次生成し（初期値は例えば０）、記憶手段１０において、Ｅステップ１３の初期値として保存する。また、未知数初期化部２４は、Ｍステップ演算に用いるための、式（２９）〜（３２）の左辺パラメータを順次生成し（初期値は例えば０）、記憶手段１０において、Ｍステップ１４の初期値として保存する。さらに、未知数初期化部２４は、ＥＭ推定値の演算に用いるための、式（３６），（２０）〜（２３）の左辺パラメータを順次生成し（初期値は例えば０）、記憶手段１０において、未知数の事後分布推定値１１に上書き保存する（更新する）。
【０１３３】
＜パワーベクトル書込部＞
パワーベクトル書込部２５は、入力制御部２２から取得したパワーベクトルｆTを、記憶手段１０において、パワーベクトル１５として順次保存する。
【０１３４】
＜観測量生成部＞
観測量生成部２６は、入力制御部２２からのコマンドにしたがって、そのときのＴに応じて、記憶手段１０から読み出したパワーベクトル１５のデータの角度ｄ毎のパワー値を、所定のルールに則って変換し、識別子ｉ（個数ｎｔ）のデータに変換することで、観測量（サンプル）を生成し、記憶手段１０に、観測量１６として順次保存する。なお、本実施形態では、所定のルールとしては、前記した式（７）、式（８）を用いる。
【０１３５】
＜終了判定部＞
終了判定部２７は、入力制御部２２からの入力信号（コマンド）が所定期間途絶えたときに、パワーベクトルｆTの入力が終了したと判定し、出力制御部２８に通知する。なお、本実施形態では、最終時刻Ｔtotalになったときに、入力が終了したと判定する。
【０１３６】
＜出力制御部＞
出力制御部２８は、終了通知を受け取ると、記憶手段１０から、最終的な推定値として、未知数の事後分布推定値１１を取得し、識別部４に出力する。
【０１３７】
［変分事後分布推論部の詳細］
≪Ｅステップ用計算部≫
Ｅステップ用計算部３１は、処理時点のＴに応じて、記憶手段１０から、サンプル、ハイパーパラメータ、Ｍステップの演算結果等を、過去を含めた演算対象時刻ｔ（ｔ≦Ｔ）に関して読み出す。そして、Ｅステップ用計算部３１は、処理時点のＴに応じて、過去を含むすべての演算対象時刻ｔ（ｔ≦Ｔ）に関して、式（２４）〜（２６）、（３７）を演算し、その演算結果を、記憶手段１０にＥステップ１３として保存する。
【０１３８】
≪Ｍステップ用計算部≫
Ｍステップ用計算部３２は、処理時点のＴに応じて、記憶手段１０から、サンプル、ハイパーパラメータ、Ｅステップの演算結果等を、過去を含めた演算対象時刻ｔ（ｔ≦Ｔ）に関して読み出す。そして、Ｍステップ用計算部３２は、処理時点のＴに応じて、過去を含むすべての演算対象時刻ｔ（ｔ≦Ｔ）に関して、式（２９）〜（３２）を演算し、その演算結果を、記憶手段１０にＭステップ１４として保存する。そして、Ｍステップ用計算部３２は、記憶手段１０に保存されたこれらの演算結果を合わせて用いて、処理時点のＴに応じて、過去を含むすべての演算対象時刻ｔ（１≦ｔ≦Ｔ）に関して、式（２０）〜（２３）、（３６）を演算し、その演算結果を、記憶手段１０に、未知数の事後分布推定値１１として上書き保存する（更新する）。
【０１３９】
≪ＥＭ収束判定部≫
ＥＭ収束判定部３３は、ＥステップとＭステップを合わせた１セットの処理の繰り返し回数ｊのしきい値（ｊmax）に達したか否かを判定することで、ＥＭアルゴリズムが収束したか否かを判別する。しきい値（ｊmax）は、予め設定される。ＥＭアルゴリズムが収束していない場合、ＥＭ収束判定部３３は、ＥステップとＭステップとを繰り返す制御を行う。ＥＭアルゴリズムが収束した場合、ＥＭ収束判定部３３は、ＥステップとＭステップとの処理を停止する制御を行う。本実施形態では、ＥＭ収束判定部３３は、収束の如何に関わらず、その時点のＴの値をパラメータ更新条件判定部３４に通知することとした。
【０１４０】
≪パラメータ更新条件判定部≫
パラメータ更新条件判定部３４は、受け取った処理時点のＴが、事前に設定された設定値（適切な正整数ｔupdateの倍数）と同じであるか否かを判別する。同じある場合には、パラメータ更新条件判定部３４は、ハイパーパラメータ更新部３５に通知する。
【０１４１】
≪ハイパーパラメータ更新部≫
ハイパーパラメータ更新部３５は、通知を受けると、記憶手段１０に保存されているハイパーパラメータ１２を任意の方法で更新する。これにより、Ｅステップ用計算部３１およびＭステップ用計算部３２は、処理時点のＴが、予め設定された適切な正整数の倍数と同じである次のタイミングから、更新されたハイパーパラメータを用いることができる。なお、パラメータ更新条件判定部３４およびハイパーパラメータ更新部３５は、前記したステップＳ２８〜Ｓ３０の処理に対応している。ただし、これらの構成は、必ずしも備える必要はない。
【０１４２】
なお、クラスタリング計算装置３は、一般的なコンピュータを、クラスタリング計算装置３を構成する前記した各手段として機能させるプログラム（クラスタリング計算プログラム）により動作させることで実現することができる。このプログラムは、通信回線を介して提供することも可能であるし、ＣＤ−ＲＯＭ等の記録媒体に書き込んで配布することも可能である。このプログラムをインストールされたコンピュータは、ＣＰＵが、ＲＯＭ等に格納されたこのプログラムをＲＡＭに展開することにより、クラスタリング計算装置３と同等の効果を奏することができる。
【０１４３】
本実施形態によれば、複数話者の会話からなる録音データから、その会話に参加した話者数、各話者の位置、さらに各話者の発話行動のタイミングを推定するダイアライゼーションにおける話者クラスタリングの問題に対して、確率的なクラスタリングを用いたので、従来のようなパラメータの設定や探索によらずに容易に話者数を推定できる。
また、本実施形態によれば、ｄＨＤＰ近似モデルを採用したことにより、時間ごとに発話に参加する話者が変化する状況も適切にモデリングできる。その結果、より精度の高い話者クラスタリングが実現できる。
さらに、本実施形態によれば、ｄＨＤＰのオンライン推定法とその高速化法を用いることで、高速に推論を行うことが可能である。なお、精度と時間のトレードオフによって実時間推論も可能である。
【０１４４】
以上、本発明の実施形態について説明したが、本発明はこれに限定されるものではなく、その趣旨を変えない範囲で実施することができる。例えば、クラスタリング計算装置３は、ｄＨＤＰ近似モデルにおいて、Ｎｏｒｍａｌ−Ｇａｍｍａ分布を仮定したが、別の分布を仮定してもよい。この場合、前記した式（３６）、式（３７）および関連する関係式と同様な計算式を、仮定された別の分布で改めて導出すればよい。
【０１４５】
また、本実施形態では、クラスタリング計算装置３は、ｄＨＤＰ近似モデルにより推論するものとしたが、例えば、ｄＨＤＰモデルや他のノンパラメトリックベイズモデルであってもよい。他のノンパラメトリックベイズモデルの場合、確率分布の連続的な時間変化をモデル化した分布であることが好ましい。
【０１４６】
また、本実施形態では、話者の位置に基づいたクラスタリングを行うために、特徴量抽出部２では、一例として、マイクロホンアレーを用いてＤＯＡ情報（音声の到達角度）を抽出するものとしたが、ダイアライゼーションに適していれば、それ以外の様々な特徴量を抽出することが可能である。例えば、話者ごとの声質の特徴に基づいたクラスタリングをする場合にはＭＦＣＣ特徴量（Mel Frequency Cepstrum Coefficient）などを抽出することができる。そして、このような特徴量であっても前記したアルゴリズムに適用することが可能である。
【実施例】
【０１４７】
本発明の効果を確認するために本実施形態に係るクラスタリング計算装置の性能を測定した。まず、第１段階の実験として、設計されたとおりのクラスタ数と適切なパラメータとを持つクラスタリングが、ｄＨＤＰによって実現できるか否かを検証する実験（クラスタリング検証実験）を行った。次に、第２段階の実験として、得られたクラスタ結果を用いたダイアライゼーション精度を評価する実験（ＤＥＲによるダイアライゼーション精度の評価）を行った。前記２種類の実験では、人工音声データと実音声データとを用いてクラスタリング計算装置の性能を確認した。
【０１４８】
＜実験データ＞
≪人工音声データ≫
人工音声データは、３人の発話者が交代で発話・非発話を切り替える状況をシミュレートしたデータである。この人工音声データは、６４［ｍｓｅｃ］の時間ステップで計算されたＤＯＡデータ（音声の到達角度のデータ）と、ＶＡＤ（音声区間検出器）による音声・非音声判定結果とからなる。また、人工音声データは、ノイズがほとんど重畳しないデータである。各実験では、ＶＡＤによる判定結果を用いて非音声区間を閾値処理で除外して４２２ステップの連続シーケンスの人工音声データを作成した。また、時間ステップごとのサンプル分布をある程度安定させるため、４２２ステップを数ステップごとに重複なくまとめて１つの長いメタステップとした。各実験では、５ステップのデータをまとめて１メタステップとし、メタステップのステップ数をｔ，Ｔに対応させた。よって、ステップ数はＴtotal＝ｃｅｉｌ（４２２／５）＝８５である。ここで、ｃｅｉｌは切り上げを示す。この人工音声データには、複数の話者が同時に喋る区間も存在している。
【０１４９】
人工音声データについての全時刻での平均パワーベクトルの分布を図９に示す。図９において、横軸は、ＤＯＡデータ、すなわち、音声到来方向角度［ｄｅｇ］を示し、縦軸は、時間平均パワー［ｄＢ］を示す。図９に示すように、各話者に対応するパワー分布の３つのピークを観察することができる。
【０１５０】
≪実音声データ≫
実音声データは、実際の複数話者の会話の様子を録音したデータである。実音声データとして、非特許文献１に記載されている４データを利用した。４データの詳細を表１に示す。表１において、ＣＰはcrossword puzzleデータ、ＤＣはdiscussionデータ、ＣＮはconversationデータを表す。
【０１５１】
【表１】
【０１５２】
実音声データは、どのデータも３００秒の音声データである。これらの実音声データについての平均パワーベクトルの分布を図１０にそれぞれ示す。図１０の横軸および縦軸は、図９のグラフと同様である。ただし、時間平均パワーのオーダーは低い。
【０１５３】
図１０（ａ）に示すＣＰ１と、図１０（ｂ）に示すＣＰ２に関しては、表１に記載の発話者数「４」に対応するパワー分布のピークの個数も「４個」なので、比較的良好なクラスタリング結果が期待される。一方、図１０（ｃ）に示すＤＣと、図１０（ｄ）に示すＣＮに関しては、表１に記載の発話者数「３」に対応するパワー分布の３つのピークが明確には観察できないため、正しいクラスタ数の推定は困難を伴うことが予想される。
【０１５４】
［クラスタリング検証実験］
まず、第１段階として、クラスタリングの性能を確認した。ここでは、ＤＰＭ（参考例１）とｄＨＤＰ（実施例１）とによって、オンラインで話者クラスタの推定を行った。このクラスタリング検証実験では、最終時刻Ｔtotalにおいて、最終的な混合比が、チャンスレベル（１／Ｋ）を超えたクラスタのみを有効なクラスタとしてカウントすることで、最終的なクラスタリング結果を求めた。そして、得られたクラスタリング結果により、ＤＰＭ（参考例１）とｄＨＤＰ（実施例１）とを比較し、話者の分布と話者数に対応したクラスタリング結果が得られているか否かを判定した。なお、チャンスレベルは偶然の一致が起こる確率である。
【０１５５】
≪人工音声データの場合≫
図９に示す人工音声データに対して、オンラインで、ＤＰＭ（参考例１）とｄＨＤＰ（実施例１）とを適用した結果を図１１に示す。図１１において、横軸は、ＤＯＡデータを、［−１８０：１８０］→［−０．５：０．５］に変換する関数を用いて無次元量に正規化した角度、すなわち、正規化角度を示す。また、縦軸は、確率密度関数（probabilistic density function：p. d. f）の値（無次元数）を示す。
【０１５６】
図１１（ｂ）に示すｄＨＤＰ（実施例１）の結果と、図９のグラフとの比較の結果、ｄＨＤＰでは、正しいクラスタ数とパラメータを得ることができたことが分かる。一方、図１１（ａ）に示すＤＰＭ（参考例１）の場合、クラスタ数が「１」となり、図９に示す人工音声データに対しては不適切な結果が得られたことを確認した。
【０１５７】
≪実音声データ（ＣＰ１、ＣＰ２）の場合≫
図１０（ａ）に示す実音声データ（ＣＰ１）に対して、オンラインで、ＤＰＭ（参考例２）とｄＨＤＰ（実施例２）とを適用した結果を図１２に示す。図１２の各軸は、図１１のグラフと同様である。図１２（ｂ）に示すｄＨＤＰ（実施例２）の結果と、図１０（ａ）のグラフとの比較の結果、ｄＨＤＰでは、正しいクラスタ数とパラメータを得ることができたことが分かる。一方、図１２（ａ）に示すＤＰＭ（参考例２）の場合、多数のクラスタに分かれてしまい、図１０（ａ）に示すＣＰ１に対しては不適切な結果が得られたことを確認した。なお、同様の傾向がＣＰ２データに対しても確認された。このときのＤＰＭ（参考例３）とｄＨＤＰ（実施例３）との結果の記述は省略した。
【０１５８】
≪実音声データ（ＤＣ）の場合≫
図１０（ｃ）に示す実音声データ（ＤＣ）に対して、オンラインで、ＤＰＭ（参考例４）とｄＨＤＰ（実施例４）とを適用した結果を図１３に示す。図１３の各軸は、図１１のグラフと同様である。図１３（ｂ）に示すｄＨＤＰ（実施例４）の結果と、図１０（ｃ）のグラフや表１における発話者数との比較の結果、ｄＨＤＰでは、話者数である３クラスタを得ることはできなかった。しかしながら、各クラスタのうち、サイズ数（図中のcluster size）の上位３つである、「cluster 4」、「cluster 6」、「cluster 14」の正規化角度の位置は、ＤＣデータにおける話者の位置に対応することができた。ここで、サイズ数（図中のcluster size）は、前記した式（４４）の右辺の分母で示される数値で定義した。なお、サイズ数の下位２つは、ノイズクラスタとなった。
【０１５９】
一方、図１３（ａ）に示すＤＰＭ（参考例４）の場合、より多くのクラスタに分かれてしまい、サイズ数上位３つのクラスタも話者の位置にそれぞれ対応することができなかった。この点からも、ｄＨＤＰ（実施例４）は、ＤＰＭ（参考例４）によるクラスタリングに比してより正確なクラスタリングを実現できたと考えられる。なお、ＤＰＭ（参考例４）の場合、各クラスタのうち、サイズ数の上位３つである、「cluster 2」、「cluster 4」、「cluster 5」の正規化角度の位置は、ＤＣデータにおける話者の位置のうちの２つにしか対応できなかった。また、第４〜６位の位置でもあと１つの話者位置に対応できなかった。また、ＣＮデータについては、ＤＣデータと同様な傾向を有すると考えられる。
【０１６０】
［ＤＥＲによるダイアライゼーション精度の評価］
第１段階のクラスタリング検証実験に続いて、第２段階では、ダイアライゼーションのためのクラスタリングとしての性能を評価するため、ＤＥＲ（diarization error ratio）による評価を試みた。ＤＥＲとはＮＩＳＴが提案した話者識別能力の指標である。具体的には、ＤＥＲは、全音声区間長に対して、以下の（１）〜（３）の３種類の誤識別区間がどれだけあったかを百分率で示したものである。ＤＥＲ値が少ないほど良いダイアライゼーションができた、という評価になる。
【０１６１】
（１）false alarm speaker time：誰も話していないのに、誰かが話したと誤検出した区間長
（２）missed speaker time：誰かが話しているのに、誰も話していないと判断した区間長
（３）speaker error time：誰かが話していることは正しく検出したが、話者を誤った区間長
【０１６２】
なお、ＤＥＲについての詳細は、下記ＵＲＬに記載されている。
「NIST Speech group,”Spring2007(RT-07) Rich Transcription Meeting Recognition Evaluation Plan”,[online],[平成２１年１月21日検索]、インターネット＜URL:http://www.nist.gov/speech/tests/rt/2007/index.html＞」
【０１６３】
第１段階のｄＨＤＰによるクラスタリングは、ダイアライゼーションのサブ問題でしかなく、そのままでは話者識別はできない。しかしながら、ｄＨＤＰによるクラスタリングでは、各時刻におけるサンプル（＝方向付きの音声パワーデータ）をクラスタリングしているため、各フレーム（時間ステップ）において各クラスタにアサインされたサンプル数を数えれば、所定の閾値を用いることで話者ごとの発声または非発声を決定することが可能である。そこで、前記した式（４３）で定義した‖ｚt,k‖を用いて、各時刻における話者ｋの発話または非発話を、式（４９）および式（５０）に示すルールで決定した。なお、所定の閾値として式（４９）に示したτＤＥＲは、適切な値に設定した。表１に示した各実音声データに対して、ｄＨＤＰ（実施例５、実施例６）について算出されたＤＥＲの結果を表２に示す。このＤＥＲ値を非特許文献１における結果（比較例）と比較した。
【０１６４】
【数２０】
【０１６５】
【表２】
【０１６６】
表２において、比較例は、既存手法である非特許文献１の中で報告された最良の結果を示す。本発明のｄＨＤＰクラスタリングの結果（実施例５、実施例６）のうち、実施例５（naive）は、式（４９）および式（５０）に示すルールのみでＤＥＲを算出したときの値を示す。また、実施例６（heuristic）は、式（４９）および式（５０）に示すルール（識別則）に加え、１フレーム内（１時間ステップ内）での同時発話人数に上限の仮定をおくことで、非発話区間のfalse alarmを低減した実施例である。実施例６で採用した方法は、非特許文献１で最良の方法が得られた識別則と共通点を持っている。すなわち、実施例６は、式（４９）に示したサンプル数閾値τＤＥＲと、同時発話人数の上限の仮定とに基づいて探索を行った結果である。そのため、実施例６は実施例５よりも良好な結果となった。さらに、実施例６は、同時発話人数の上限の仮定をした比較例と比較しても、明らかに、ＤＥＲ値として良好な値を示すことがわかる。
【０１６７】
以上の検証実験および評価実験をまとめると、ｄＨＤＰで解決されるのはサブ問題（第１段階）であるクラスタリングであるが、第２段階で、クラスタリングの過程で計算されるサンプルアサインメントｚを用いることで、ダイアライゼーションの精度を示すＤＥＲ指標の向上を実現可能である。すなわち、本発明によれば、ｄＨＤＰでクラスタリングの問題を解決したことで、自然と良好なダイアライゼーションが可能となったと結論できる。
【符号の説明】
【０１６８】
１ ダイアライゼーションシステム
２ 特徴量抽出部（特徴量抽出手段）
３ クラスタリング計算装置
４ 識別部（識別手段）
１０ 記憶手段（推定値記憶手段、観測量記憶手段）
２１ パワーベクトル読込部（読込手段）
２２ 入力制御部
２３ 割当部
２４ 未知数初期化部
２５ パワーベクトル書込部
２６ 観測量生成部（観測量生成手段）
２７ 終了判定部
２８ 出力制御部（出力制御手段）
３０ 変分事後分布推論部（事後分布推論手段）
３１ Ｅステップ用計算部（Ｅステップ用計算手段）
３２ Ｍステップ用計算部（Ｍステップ用計算手段）
３３ ＥＭ収束判定部（収束判定手段）
３４ ハイパーパラメータ更新条件判定部
３５ ハイパーパラメータ更新部
【技術分野】
【０００１】
本発明は、話者数が未知である会話の録音データから話者数を推定する技術に関する。
【背景技術】
【０００２】
従来、複数話者の会話からなる録音データから、その会話に参加した話者数とさらに各話者の発話したタイミングを推定する問題が知られている。この問題やこの問題を解決する技術は、ダイアライゼーション（diarization）と呼ばれている。ダイアライゼーションは、簡単には「いつ、誰が発話したか」を自動推定する技術といえる。この技術の応用としては、会議データへのアノテーション（annotation）やそれを用いた検索、音声強調など幅広い用途が期待されている。
【０００３】
既存のダイアライゼーションの解法は、大きく２つに分けることができる。第１の方法は、話者固有の声質を推定して話者を区別する方法であり、第２の方法は、話者の位置を推定する方法である。このうち、第１の方法（話者声質を利用する方法）では、各話者の声の特徴を抽出することで現在発話している話者を識別する。この第１の方法では、話者が移動したとしても話者を識別できるという利点があるが、複数の話者が同時に発話した場合には話者の識別が困難となるという欠点がある。
【０００４】
第２の方法（話者位置に関する情報を利用する方法）は、話者の位置を推定することで話者数とその位置を推定する方法である（例えば、非特許文献１〜３参照）。非特許文献１〜３の手法では、マイクロホンアレーを用いて各話者の位置を推定することで話者の識別を行う。従って、非特許文献１〜３の手法は、話者が移動すると同一話者として識別することができないという欠点があるが、複数話者が同時に発話する場合でも各話者の発話行動を識別することができるという利点がある。
【０００５】
ダイアライゼーションでは、一般に、録音された音声データ中の発話者数は未知であり、推定しなければならない。また、第１の方法（話者声質を利用する方法）によるもの、第２の方法（話者位置に関する情報を利用する方法）によるものに関わらず、各話者を識別するためには、各話者を特徴づける量（パラメータ）を推定する必要がある。これは、いわゆるクラスタリングの問題に相当する。クラスタリングの問題とは、観測されたデータを適切な数のクラスタに分割し、各クラスタのパラメータを推定する問題である。ダイアライゼーションを単純にクラスタリングの問題に置き換えようとしても、そのクラスタ数とデータの分割は未知なので、パラメータと共に学習する必要がある。ダイアライゼーションの問題では、クラスタ数は「話者数」を表し、各クラスタのパラメータは「話者の特徴量」を表し、結果得られるデータの分割結果から「各話者の発話タイミング」が示唆されることになる。
【０００６】
非特許文献１ではleader-followerアルゴリズムと呼ばれるクラスタリング手法を用いている。これは、逐次的なクラスタリング手法であり、少しずつサンプルを入力していきながらクラスタリングを行うものである。このクラスタリング手法では、クラスタリングの際に、新しく入力されたサンプルが既存のどのクラスタからも一定距離以上離れていた場合に、そのサンプルを中心として新しいクラスタを生成する。なお、leader-followerアルゴリズムについては、「R. O. Duda, P. E. Hart and D. G. Stork, “Pattern Classification”, John Wiley & Sons, 2001.」に記載されている。
【０００７】
また、非特許文献３では、ＢＩＣ基準と呼ばれる評価値を最大化するようなクラスタ数を選択する方法を提案している。この方法では、クラスタ数をＫに設定してから実際にクラスタリングを行って評価値を計算し、次はクラスタ数をＫ＋１に設定して評価値を計算する、ということを繰り返えして評価値を最大化するクラスタ数を探索する。
【０００８】
ところで、クラスタ数未知のデータに対するクラスタリングモデルとして、ノンパラメトリックベイズモデルが、近年、多くの場面で利用されるようになってきた。例えば、ノンパラメトリックベイズモデルの１種であるDirichlet Process Mixture（ＤＰＭ）は、クラスタ数と各サンプルのクラスタリングを確率的に同時に最適化することができる。従って、ＤＰＭは、既存のクラスタリングモデルのように、クラスタ数の最適化を簡単に実現できる点が大きな特徴である。このＤＰＭにおいて、確率分布の連続的な時間変化をモデル化した拡張モデルとしては、dynamic Hierarchical Dirichlet Process（ｄＨＤＰ）と呼ばれるモデルが知られている（非特許文献４参照）。
【先行技術文献】
【非特許文献】
【０００９】
【非特許文献１】S. Araki, M. Fujimoto, K. Ishizuka, H. Sawada and S. Makino, ”A DOA Based Speaker Diarization System for Real Meetings”, Proceedings of the Joint Workshop on Hands-Free Speech Communication and Microphone Arrays, p.29-32, 2008.
【非特許文献２】荒木章子、藤本雅清、石塚健太郎、澤田宏、牧野昭二、“音声区間検出と方向情報を用いた会議音声話者識別システムとその評価”、日本音響学会講演論文集、春季、vol. 1-10-1, p. 1-4, 2008.
【非特許文献３】J. M. Pardo, X. Anguera and C. Wooters, “Speaker Diarization for Multi-Microphone Meetings Using Only Between-Channel Differences”, Proceedings of the Third Joint Workshop on Multimodal Interaction and Related machine Learning Algorithms”, p. 257-264, 2008.
【非特許文献４】L. Ren, D. B. Dunson and L. Carin, “The Dynamic Hierarchical Dirichlet process”, Proceedings of International Conference on Machine Learning, p. 824-831, 2008.
【発明の概要】
【発明が解決しようとする課題】
【００１０】
従来のダイアライゼーションの研究において、例えば、非特許文献１に記載の手法は、計算が簡単で高速に実行できるが、新しいクラスタを生成するための距離閾値の設定が必要である。この閾値の設定によって、最終的に得られるクラスタリングとクラスタ数が決定される。一方で、クラスタ数やクラスタリング結果の推定値を真の値に近づけるという意味で、この閾値を最適化することは困難である。
【００１１】
また、例えば、非特許文献３に記載の方法は、実際には不適切なクラスタ数のもとでもクラスタリングを行う必要があり、処理の過程において、計算量や時間の観点からは大きな無駄が発生する。
【００１２】
したがって、従来のダイアライゼーションの研究では、話者数と話者の特徴を推定する部分、つまりクラスタリングの問題に改良の余地があった。また、ダイアライゼーションに対して、処理負荷の低減や処理の高速化が要望されている。さらに、従来のダイアライゼーションの研究では、ｄＨＤＰやＤＰＭなどのノンパラメトリックベイズモデルは、クラスタリングモデルとして採用されておらず、その適用方法が知られていなかった。
【００１３】
そこで、本発明では、前記した問題を解決し、ダイアライゼーションにおいて、話者数と各話者を特徴づけるパラメータを正確に推定するクラスタリング技術を提供することを目的とする。
【課題を解決するための手段】
【００１４】
前記目的を達成するために、本願発明者らは、ダイアライゼーションにおいて、話者数と各話者を特徴づけるパラメータを推定するクラスタリングにおいて種々検討を行った。その結果、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化できるモデルとしてノンパラメトリックベイズモデルを採用したときに、話者数を正確に推定できることを見出した。
【００１５】
そこで、本発明に係るクラスタリング計算装置は、話者数が未知である会話の録音データから前記会話の話者数を推定するために、各話者を特徴付ける特徴量を抽出する特徴量抽出手段と、前記抽出された特徴量から前記各話者に対応する複数のクラスタを生成するときの複数の未知パラメータの値をそれぞれ推定するクラスタリング計算装置と、前記推定された複数のパラメータ値により前記会話の各話者を識別する識別手段とを有したダイアライゼーションシステムの前記クラスタリング計算装置であって、前記抽出された特徴量を読み込む読込手段と、前記読み込んだ特徴量をノンパラメトリックベイズモデルに対応したベクトルの観測量に変換する観測量生成手段と、前記変換された観測量の集合データを蓄積記憶する観測量記憶手段と、前記観測量の集合データから複数のクラスタをノンパラメトリックベイズモデルにより生成するときの複数のパラメータの事後分布の値をＥＭアルゴリズムによりそれぞれ推定および更新する事後分布推論手段と、前記推定および更新された複数のパラメータの事後分布の値を蓄積記憶する推定値記憶手段と、事前に設定された終了条件が成立したときに前記推定値記憶手段に記憶されている前記複数のパラメータの事後分布の最新の推定値を出力する出力制御手段と、を備えることを特徴とする。
【００１６】
かかる構成によれば、クラスタリング計算装置は、会話の録音データから抽出された特徴量をノンパラメトリックベイズモデルに適応できるようにベクトルの観測量に変換することで生成した観測量の集合データを用いて、ノンパラメトリックベイズモデルに従った演算により、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化する。そして、クラスタリング計算装置は、この最適化において、ＥＭアルゴリズムにより事後分布の推定および更新を繰り返し、蓄積記憶する。ここで、ＥＭアルゴリズムは、局所最適解の計算アルゴリズムであるため、計算を繰り返すことで、ある１つの解に収束させることができる。また、ノンパラメトリックベイズモデルでは、予め準備した最大クラスタ数の個数のクラスタが最適化処理の過程で、有効なクラスタ数の個数に収束し、他のクラスタは、クラスタの混合比がほぼ０となる。このクラスタリング結果が得られれば、各クラスタへ帰属するサンプル数（観測量の個数）が計算できる。また、この各クラスタへ帰属するサンプル数を用いて各クラスタの混合比を計算可能である。さらに、この混合比により、有効なクラスタ数を決定することができる。ここで、クラスタ数は、会話の録音データ中の話者に対応しているので、話者数が決定できるようになる。また、これらクラスタリング結果と共に、それを用いた各種変数の推定値も同様に蓄積記憶することができる。そして、クラスタリング計算装置は、終了条件が成立したときに蓄積記憶されている最新の推定値を出力する。
【００１７】
また、本発明に係るクラスタリング計算装置は、前記事後分布推論手段が、前記ＥＭアルゴリズムのＥステップの処理として、ｄＨＤＰ（dynamic Hierarchical Dirichlet Process）モデルにおいて予め定められた事前分布、観測分布および最大クラスタ数と、ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＭステップまでに推定された隠れ変数の事後分布の推定値とを読み込んで、過去分を含めた演算対象時刻ごと、かつ、クラスタごと、かつ、演算対象時刻別の全データごとに、前記クラスタと、混合比と、前記クラスタ分布の時間変化の程度を表す重みとに関するパラメータの事後分布の値を演算するＥステップ用計算手段と、前記ＥＭアルゴリズムのＭステップの処理として、前記ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＥステップまでに推定されたパラメータの事後分布の推定値とを読み込んで、過去分を含めた演算対象時刻ごと、かつ、演算対象時刻別の全データごとに、２種類の隠れ変数の事後分布の値を推定し、前記２種類の隠れ変数のうち、前記クラスタ分布の時間変化の程度を表す重みに関連付けられた第１隠れ変数の事後分布の値についてはクラスタごとに演算し、前記第１隠れ変数および前記混合比に関連付けられた第２隠れ変数の事後分布の値については、演算対象時刻から過去に遡及した時刻ごとに演算するＭステップ用計算手段と、前記Ｅステップの処理と前記Ｍステップの処理とを交互に予め定められた回数だけ繰り返し実行させる制御を行う収束判定手段と、を備えることを特徴とする。
【００１８】
かかる構成によれば、クラスタリング計算装置は、ノンパラメトリックベイズモデルのうちｄＨＤＰモデルに従った演算により、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化するためにＥＭアルゴリズムを用いる。また、ｄＨＤＰモデルでは、ある時刻ｔにおけるパラメータ分布を、前の時刻ｔ−１でのパラメータ分布と、当該時刻ｔでのパラメータ分布の変化量とを、時刻ｔにおける変化量の変化の割合（重み）で重ね合わせた分布として表現している。すなわち、ｄＨＤＰでは、時刻に依存して分布が少しずつ変化するというモデルとなっている。クラスタリング計算装置のアルゴリズムは、ある時刻ｔにおけるパラメータ分布を、時刻ｔに実際に発話した発話者を構成しているものとして構築した。そのため、ｄＨＤＰモデルは、話者の入れ替わり（turn-taking）によって観測量（サンプル）の分布が時間変化するダイアライゼーションのタスクにとって好都合の確率モデルとなった。したがって、話者の入れ替わりがある会話録音データから、未知の話者数と各話者を特徴づけるパラメータを正確に推定することができる。
【００１９】
また、本発明に係るクラスタリング計算装置は、前記Ｅステップ用計算手段が、観測量の集合データと、隠れ変数の事後分布の推定値とについて、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込み、推定対象の時刻から事前に設定された設定値分だけ遡及した過去の時刻までを演算対象時刻として、前記クラスタと、前記混合比と、前記重みとに関するパラメータの事後分布の値を演算し、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、事前に設定された再推定すべきクラスタの判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合に、前記混合比に関するパラメータの事後分布の値を演算し、前記Ｍステップ用計算手段が、観測量の集合データと、パラメータの事後分布の推定値とについて、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込み、推定対象の時刻から事前に設定された設定値分まで遡及した過去の時刻までを演算対象時刻として、前記第１隠れ変数および第２隠れ変数の事後分布の値を演算し、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、前記判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合に、前記第１隠れ変数の事後分布の値を演算することを特徴とする。
【００２０】
かかる構成によれば、クラスタリング計算装置は、ＥステップおよびＭステップにおいて、推定対象の時刻から事前に定められた過去の一時点まで遡った観測量の集合データと、その演算に必要な推定値とを読み込み、当該推定対象の時刻から事前に定められた設定値分だけ過去に遡った過去の時刻までを演算対象時刻として、推定値を求める演算を行う。したがって、事前に過去の一時点や設定値を定めることなく、時間ステップの進展とともに推定すべき変数の個数が増加していくだけである場合と比べて、処理負荷の低減や処理の高速化を実現できる。また、クラスタリング計算装置は、ＥステップおよびＭステップにおいて、クラスタごとの演算のためにクラスタを更新したときに、必要な場合にだけそのクラスタにおける推定値の再推定処理を行う。したがって、ｄＨＤＰモデルにおいて予め設定される最大クラスタ数の個数のクラスタに関して、推定値の再推定処理を毎回実行する場合に比べて、処理負荷の低減や処理の高速化を実現できる。
【００２１】
また、課題を解決するため、本発明に係るクラスタリング計算方法は、話者数が未知である会話の録音データから前記会話の話者数を推定するダイアライゼーションシステムにおいて、記憶手段と、読込手段と、観測量生成手段と、事後分布推論手段と、出力制御手段とを備えて、前記録音データから抽出された各話者を特徴付ける特徴量から前記各話者に対応する複数のクラスタを生成するときの複数の未知パラメータの値をそれぞれ推定するクラスタリング計算装置のクラスタリング計算方法であって、前記読込手段によって、前記抽出された特徴量を読み込む特徴量読込ステップと、前記観測量生成手段によって、前記読み込んだ特徴量をノンパラメトリックベイズモデルに対応したベクトルの観測量に変換し、前記変換された観測量の集合データを記憶手段に順次蓄積する観測量蓄積ステップと、前記事後分布推論手段によって、前記観測量の集合データから複数のクラスタをノンパラメトリックベイズモデルにより生成するときの複数のパラメータの事後分布の値をＥＭアルゴリズムによりそれぞれ推定し、当該推定値を前記記憶手段に順次格納および更新する事後分布推定ステップと、前記出力制御手段によって、事前に設定された終了条件が成立したときに前記記憶手段に記憶されている前記複数のパラメータの事後分布の最新の推定値を出力する推定値出力ステップと、を含んで実行することを特徴とする。
【００２２】
かかる手順によれば、クラスタリング計算方法において、クラスタリング計算装置は、まず、会話の録音データから抽出された特徴量を読み込み、ノンパラメトリックベイズモデルに適応できるようにベクトルの観測量に変換し、蓄積する。そして、クラスタリング計算装置は、生成した観測量の集合データを用いて、ノンパラメトリックベイズモデルに従った演算により、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化する。そして、クラスタリング計算装置は、この最適化において、ＥＭアルゴリズムにより事後分布の推定および更新を繰り返し、蓄積記憶する。そして、クラスタリング計算装置は、終了条件が成立したときに蓄積記憶されている最新の推定値を出力する。
【００２３】
また、本発明に係るクラスタリング計算方法は、前記事後分布推論手段が、前記事後分布推定ステップにおいて、前記ＥＭアルゴリズムのＥステップの処理として、ｄＨＤＰモデルにおいて予め定められた事前分布、観測分布および最大クラスタ数と、ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＭステップまでに推定された隠れ変数の事後分布の推定値とを読み込んで、推定対象の時刻に過去分を含めた演算対象時刻ごと、かつ、クラスタごと、かつ、前記演算対象時刻別の全データごとに、前記クラスタと、混合比と、前記クラスタ分布の時間変化の程度を表す重みとに関するパラメータの事後分布の値を演算する段階と、前記ＥＭアルゴリズムのＭステップの処理として、前記ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＥステップまでに推定されたパラメータの事後分布の推定値とを読み込んで、推定対象の時刻に過去分を含めた演算対象時刻ごと、かつ、前記演算対象時刻別の全データごとに、２種類の隠れ変数の事後分布の値を推定し、前記２種類の隠れ変数のうち、前記クラスタ分布の時間変化の程度を表す重みに関連付けられた第１隠れ変数の事後分布の値についてはクラスタごとに演算し、前記第１隠れ変数および前記混合比に関連付けられた第２隠れ変数の事後分布の値については、前記演算対象時刻から過去に遡及した時刻ごとに演算する段階とを含み、前記Ｅステップの処理と前記Ｍステップの処理とを交互に予め定められた回数だけ繰り返し実行することを特徴とする。
【００２４】
かかる手順によれば、クラスタリング計算方法において、クラスタリング計算装置は、ノンパラメトリックベイズモデルのうちｄＨＤＰモデルに従った演算により、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化するためにＥＭアルゴリズムを用いる。ここで、ｄＨＤＰでは、時刻に依存して分布が少しずつ変化するというモデルとなっているため、話者の入れ替わり（turn-taking）によって観測量（サンプル）の分布が時間変化するダイアライゼーションのタスクにとって恰好の確率モデルである。したがって、話者の入れ替わりがある会話録音データから、未知の話者数と各話者を特徴づけるパラメータを正確に推定することができる。
【００２５】
また、本発明に係るクラスタリング計算方法は、前記事後分布推論手段が、前記Ｅステップにおいて、観測量の集合データと、隠れ変数の事後分布の推定値とについては、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込み、前記Ｍステップにおいて、観測量の集合データと、パラメータの事後分布の推定値とについては、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込むことを特徴とする。
【００２６】
かかる手順によれば、クラスタリング計算方法において、ＥステップおよびＭステップでは、当該推定対象の時刻から事前に定められた過去の一時点まで遡った観測量の集合データと、その演算に必要な推定値とを読みこむ。したがって、事前に過去の一時点を定めることなく、時間ステップの進展とともに推定すべき変数の個数が増加していくだけである場合と比べて、処理負荷の低減や処理の高速化を実現できる。
【００２７】
また、本発明に係るクラスタリング計算方法は、前記事後分布推論手段が、前記Ｅステップにおいて、推定対象の時刻から事前に設定された設定値分だけ遡及した過去の時刻までを演算対象時刻として、前記クラスタと、前記混合比と、前記重みとに関するパラメータの事後分布の値を演算し、前記Ｍステップにおいて、推定対象の時刻から事前に設定された設定値分まで遡及した過去の時刻までを演算対象時刻として、前記第１隠れ変数および第２隠れ変数の事後分布の値を演算することを特徴とする。
【００２８】
かかる手順によれば、クラスタリング計算方法において、ＥステップおよびＭステップでは、当該推定対象の時刻から事前に定められた設定値分だけ過去に遡った過去の時刻までを演算対象時刻として、推定値を求める演算を行う。したがって、事前に設定値を定めることなく、時間ステップの進展とともに推定すべき変数の個数が増加していくだけである場合と比べて、処理負荷の低減や処理の高速化を実現できる。
【００２９】
また、本発明に係るクラスタリング計算方法は、前記事後分布推論手段が、前記Ｅステップにおいて、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、事前に設定された再推定すべきクラスタの判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合にだけ、前記混合比に関するパラメータの事後分布の値を演算し、前記Ｍステップにおいて、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、前記判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合にだけ、前記第１隠れ変数の事後分布の値を演算することを特徴とする。
【００３０】
かかる手順によれば、クラスタリング計算方法において、ＥステップおよびＭステップでは、クラスタごとの演算のためにクラスタを更新したときに、必要な場合にだけそのクラスタにおける推定値の再推定処理を行う。したがって、ｄＨＤＰモデルにおいて予め設定される最大クラスタ数の個数のクラスタに関して、推定値の再推定処理を毎回実行する場合に比べて、処理負荷の低減や処理の高速化を実現できる。
【００３１】
また、本発明に係るクラスタリング計算プログラムは、前記いずれかのクラスタリング計算装置を構成する各手段としてコンピュータを機能させるためのプログラムである。このように構成されることにより、このプログラムをインストールされたコンピュータは、このプログラムに基づいた各機能を実現することができる。
【００３２】
また、本発明に係るコンピュータ読み取り可能な記録媒体は、前記クラスタリング計算プログラムが記録されたことを特徴とする。このように構成されることにより、この記録媒体を装着されたコンピュータは、この記録媒体に記録されたプログラムに基づいた各機能を実現することができる。
【発明の効果】
【００３３】
本発明によれば、ダイアライゼーションにおける話者クラスタリングの問題に対して、ノンパラメトリックベイズモデルを採用して、確率的なクラスタリングを用いることで、従来のようなパラメータの設定や探索によらずに容易に話者数を推定できる。また、本発明によれば、ノンパラメトリックベイズモデルを採用したので、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化できる。その結果、ダイアライゼーションにおいて、話者数と各話者を特徴づけるパラメータを正確に推定できる。
【図面の簡単な説明】
【００３４】
【図１】本発明の実施形態に係るクラスタリング計算装置を含むダイアライゼーションシステムの概要を示す構成図である。
【図２】本発明の実施形態に係るクラスタリング計算方法の全体処理の流れを示すフローチャートである。
【図３】図２に示す変分事後分布推論手順を示すフローチャートである。
【図４】図３に示すＥステップの計算手順の一例を示すフローチャートである。
【図５】図３に示すＭステップの計算手順の一例を示すフローチャートである。
【図６】本発明の実施形態に係るクラスタリング計算方法におけるＥステップの計算手順を示すフローチャートである。
【図７】本発明の実施形態に係るクラスタリング計算方法におけるＭステップの計算手順を示すフローチャートである。
【図８】本発明の実施形態に係るクラスタリング計算装置の構成の一例を示す機能ブロック図である。
【図９】本発明の実施形態に係るクラスタリング計算装置のクラスタリング性能を評価するために用いた人工音声データの時間平均パワー分布を示すグラフである。
【図１０】本発明の実施形態に係るクラスタリング計算装置のクラスタリング性能を評価するために用いた実音声データの時間平均パワー分布を示すグラフであって、（ａ）はＣＰ１、（ｂ）はＣＰ２、（ｃ）はＤＣ、（ｄ）はＣＮをそれぞれ示している。
【図１１】本発明の実施形態に係るクラスタリング計算装置によって人工音声データをクラスタリングした結果を示すグラフであって、（ａ）はＤＰＭ、（ｂ）はｄＨＤＰをそれぞれ示している。
【図１２】本発明の実施形態に係るクラスタリング計算装置によってＣＰ１の実音声データをクラスタリングした結果を示すグラフであって、（ａ）はＤＰＭ、（ｂ）はｄＨＤＰをそれぞれ示している。
【図１３】本発明の実施形態に係るクラスタリング計算装置によってＤＣの実音声データをクラスタリングした結果を示すグラフであって、（ａ）はＤＰＭ、（ｂ）はｄＨＤＰをそれぞれ示している。
【発明を実施するための形態】
【００３５】
図面を参照して本発明のクラスタリング計算装置およびクラスタリング計算方法を実施するための形態（以下「実施形態」という）について詳細に説明する。以下では、推論原理の概略、ダイアライゼーションシステムの概略、クラスタリング計算方法の概要、計算アルゴリズム、クラスタリング計算装置について順次説明する。
【００３６】
［推論原理の概略］
本実施形態では、ノンパラメトリックベイズモデルの一種として、例えば、ｄＨＤＰを用いることとする。ここでは、ｄＨＤＰを簡単に説明する。ｄＨＤＰの数学的なモデルを式（１）〜式（６）に示す。なお、ＤＰＭ等の他のノンパラメトリックベイズモデルを用いてもよいことはもちろんである。ＤＰＭの詳細は、例えば、「上田修功、山田武士、“ノンパラメトリックベイズモデル”，応用数理，Vol. 17， No. 3， pp. 196-214, 2007.」に記載されている。
【００３７】
【数１】
【００３８】
ここで、〜は確率分布からのサンプリングを表す。また、ＤＰ（・）はDirichlet Process（無限次元分布）を表し、γ，α０，ａ０，ｂ０は事前に設定するハイパーパラメータである。
【００３９】
ｄＨＤＰでは、まず、式（１）で離散無限個のパラメータ分布（クラスタ）Ｇ０を生成する。ダイアライゼーションにおいて、この分布Ｇ０は、データ全体を見たときの各発話者の発話割合と構成に相当する。ＤＰＭに代表されるノンパラメトリックベイズモデルでは、推定された解が自動的に少数のパラメータ（話者）からなる分布に集約される。
【００４０】
式（２）のＨtは、時刻ｔでの話者分布の時間変化（分布変化）を表す。
式（３）のｗtは、時刻ｔでの話者分布の時間変化の割合（程度）を表す。
式（５）のθt iは、時刻ｔでのクラスタのパラメータを表し、式（６）のｘt iは、時刻ｔでのサンプル分布を表す。なお、ｉは、時刻ｔでのｉ番目のデータを示す。
【００４１】
式（４）におけるＧtは、時刻ｔにおけるパラメータ分布である。ダイアライゼーションにおいて、分布Ｇtは、時刻ｔに実際に発話した発話者の構成を表現している。式（４）は、このＧtを、時刻ｔ−１での分布Ｇt-1と、時刻ｔでの分布変化を表すＨtとをｗtで重ね合わせた分布として表現したものである。よって、ｄＨＤＰでは、各時刻でのサンプル分布（ｘt iの分布）に変動が許される。一方で、時刻に不変な分布であるＧ０も推定するため、全体を通した話者クラスタも学習されている。
【００４２】
さらに、ｄＨＤＰでは、変化の割合であるｗtも動的に学習されるので、話者が交代するときには、劇的に分布が代わり、そうでないときには、ほとんど分布が変化しないというように、変化の割合が一定でないデータのモデル化にも対応できる。このように、ｄＨＤＰでは、時刻に依存して分布が少しずつ変化するというモデルとなっているため、話者の入れ替わり（turn-taking）によって観測量（サンプル）の分布が時間変化するダイアライゼーションのタスクにとって恰好の確率モデルとなっている。
【００４３】
［ダイアライゼーションシステムの概略］
図１は、本発明の実施形態に係るクラスタリング計算装置を含むダイアライゼーションシステムの概要を示す構成図である。本実施形態のダイアライゼーションは、前記した第２の方法（話者位置に関する情報を利用する方法）によるものとして説明する。予め、未知数話者による会話を録音し、ダイアライゼーションシステム１への入力とする。ここでは、室内１０１に３人の話者ＨA、ＨB、ＨCが図１のような定位置で会話したものとする。音声データ（会話の録音データ）１０２は時系列データとなる。
【００４４】
ダイアライゼーションシステム１は、１つの大型コンピュータ、または、複数のコンピュータからなる。ここでは、ダイアライゼーションシステム１は、３つのコンピュータ、すなわち、特徴量抽出部２と、クラスタリング計算装置３と、識別部４とを備える。
【００４５】
特徴量抽出部２は、ノイズ除去のような前処理を行うと共に、ダイアライゼーションに適した各種特徴量を抽出するものである。特徴量抽出部２は、例えば、マイクロホンアレーから取得した録音データから、ＤＯＡ（direction of arrival：音声到来方向）情報を抽出し、クラスタリング計算装置３に出力する。ＤＯＡ情報（音声の到達角度）は、マイクに対してどの方向からどの程度の強さの音声信号が観測されたかを推定した量になっている。
【００４６】
クラスタリング計算装置３は、ＤＯＡ情報に基づいてクラスタリング処理によって話者数と各話者を特徴づけるパラメータとを推定する。すなわち、クラスタリング計算装置３は、抽出した音声特徴量をクラスタリングし、クラスタ数と各クラスタのパラメータを推定する。クラスタリング計算装置３は、確率的なクラスタリングモデルを適用することで、クラスタ数やデータの分割、そしてパラメータを事後確率最大化の意味で同時に最適化する。
【００４７】
識別部４は、クラスタリング計算装置３で得られたクラスタリング結果（クラスタリング推定値）により、各時刻における話者の発話状態を識別する。この識別部４は、クラスタリング推定値を解析して、クラスタ数と、その位置とを、ユーザが識別できるような画面表示で提示する。
【００４８】
図１に例示したダイアライゼーション結果表示１０３のグラフでは、横軸が時刻（秒）、縦軸が方向（話者の位置）をそれぞれ示している。この例では、３人の話者ＨA、ＨB、ＨCに対応して、３つの方向に対して矩形波形が表示されている。各矩形波形の山となっている部分が各話者の発話（発言）を表している。まず話者ＨCが発言を終えると、次に話者ＨAが発言し、その発言途中に話者ＨCが再び発言し始めると、話者ＨAが黙り、話者ＨBが発言を開始する、というように、話者の入れ替わり（turn-taking）が生じていることがわかる。なお、２人の話者が同時に発話するタイミングも生じていることがわかる。
【００４９】
［クラスタリング計算方法の概要］
ここでは、推論原理のダイアライゼーションへの導入と、クラスタリングモデルとについて説明する。
＜推論原理のダイアライゼーションへの導入＞
ここでは、クラスタリング計算装置３において、ノンパラメトリックベイズモデル（ｄＨＤＰ）を用いるため、その前段の特徴量抽出部２で抽出する特徴量を定式化する。時刻ｔにおいて、角度ｄ（例えば、ｄ＝−１８０，−１７９，…，０，…，１８０）方向から聞こえてきた音声パワー（ＤＯＡ情報）をｆt dとする。すなわち、各時刻ｔにおける音声パワーベクトルは
である。なお、ｄ＝−１８０の方向と、ｄ＝１８０の方向とは同じものである。
【００５０】
このパワーベクトルを、本実施形態で用いるクラスタリングモデルに適合するように、１次元ベクトルｘt iの集合へと変換する。ここでは、各ｆt dについて、閾値定数τとパラメータμとを用いて、値ｇ（ｄ）をもつｎt d個のサンプル集合を生成することとした。ここで、関数ｇ（・）は、実装の都合に合わせて選択された適切なスケール関数を示す。例えば、関数ｇ（・）は、［−１８０：１８０］→［−０．５：０．５］に変換する関数を用いることができる。また、ｎt dは、式（７）で定義する。なお、１次元に限らず、例えば、３次元ベクトルｘt iの集合へ変換するように構成してもよい。
【００５１】
【数２】
【００５２】
式（７）の関数ｈ（・）は、パワーの値に応じて何らかの正整数を返す関数であり、定数でもよい。本実施形態では、例えば、ｈ＝１を用いる。以上の量子化変換を全ての角度ｄに関して行い、時刻ｔの観測量をサンプル集合に集約する。
【００５３】
【数３】
【００５４】
本実施形態で扱うクラスタリング問題は、サンプル集合Ｘt＝｛ｘ１，…，ｘt｝のクラスタリングと捉えることができる。人間の発話音声パワーは背景ノイズよりもはるかに大きいものとすれば、話者の位置を反映した値のサンプルが大量に観測されるため、主要なクラスタとして話者の位置が推定できるはずである。そこで、本実施形態のクラスタリング計算方法では、クラスタ数Ｋを推定しながら、同時に各サンプルｘt iに対してどのクラスタに属するかを表す変数ｚt iの確率分布を求めること、さらにＫ個のクラスタに対応するパラメータΘ＝｛θk｝を求めることをそれぞれ行う。
【００５５】
＜クラスタリングモデル＞
本実施形態では、計算量とアルゴリズムの簡易さとを考慮して、ｄＨＤＰ近似モデルを用いる。ｄＨＤＰ近似モデルの生成モデルは、以下の式（９）〜式（１５）のようになる。なお、ｄＨＤＰ近似モデルについては、「I. Pruteanu-Malinici, L. Ren, J. Paisley, E. Wang and L. Carin, “Dynamic Hierarchical Dirichlet Process for Modeling Topics in Time-Stamped Documents”, IEEE Transactions on Pattern Analysis and Machine Intelligence, submitted, 2008.」に記載されている。
【００５６】
【数４】
【００５７】
ｄＨＤＰ近似モデルでは、まず、最大クラスタ数Ｋを固定しておく。最大クラスタ数Ｋは、推定しようとしている話者数よりも充分大きい値（例えば、数十〜１００）であればよい。実際には、「有効な」クラスタ数Ｋeff（＜Ｋ）なる個数のクラスタが、推定しようとしている話者に対応することとなる。「有効な」クラスタ数Ｋeffについては、学習の結果、話者に対応しないそれ以外のクラスタの重み（混合比）が自動的にほぼ０になることを利用して判断することができる。すなわち「有効な」クラスタ数Ｋeffは、推定の過程で自動的に決定される。
【００５８】
≪式（９）、式（１０）≫
ｄＨＤＰ近似モデルでは、式（９）で、Ｋ個のクラスタに対応するパラメータをサンプリングする。式（１０）では、有限次元のDirichlet分布よりinnovation measure Ｈt（後記する式（１６ｃ）参照）のサンプリングを行う。より具体的には混合比πtをサンプリングする。これは、前記した式（４）より、式（１６ａ）およびその変形式である式（１６ｂ）が導かれるため、時刻ｔにおける話者分布Ｇtは、時刻ｔまでに生成したＨ1:lの重ね合わせだけで表現できるからである。ここで、「Ｈ1:l」は、Ｈ1〜Ｈlを表す。
【００５９】
【数５】
【００６０】
≪式（１１）、式（１２）≫
ｄＨＤＰ近似モデルでは、続いて、式（１１）により、話者分布Ｇtの時間変化の程度を表すｗtをサンプリングして、このｗtを用いて、式（１２）で定義されるｖtl（ｌ＝１，…，ｔ）を計算する。ここで、時刻ｌは、時刻ｔとそれよりも過去の時刻とを表す。
【００６１】
≪式（１３）≫
式（１３）に示す隠れ変数ｄt iは、ｔ次元の｛０，１｝ベクトルである（ｔは時刻、要素の値は０か１のみ）。隠れ変数ｄt iは、時刻ｔのときのｔ次元の要素のうち、ｌ次元目の要素ｄt i lが１の値をもつものである（ｌ≦ｔ）。隠れ変数ｄt iのｌ次元目の要素ｄt i lは、時刻ｌでの分布変化Ｈlに対応する要素である。時刻ｌでの分布変化Ｈlは、時刻ｔのｉ番目のサンプルｘt iをサンプルするための分布変化である。
【００６２】
≪式（１４）≫
式（１４）に示す隠れ変数ｚt iは、同様に、Ｋ次元の｛０，１｝ベクトルである。
隠れ変数ｚt iは、実際にサンプルｘt iをサンプリングするクラスタ（パラメータ）ｋに対応するｋ次元目の要素ｚt i kにのみ１の値を持つ。
【００６３】
≪式（１５）≫
式（１５）に示す観測量ｘt iは、与えられたクラスタ番号（ｋ）に対応するパラメータθから生成される。式（１５）は、前記した式（６）、式（８）と同様であって、別の表式である。
【００６４】
≪観測分布と事前分布≫
ｄＨＤＰ近似モデルでは、式（１５）に示す観測分布Ｆと、パラメータの事前分布Ｈ（式（９）参照）とを事前に定めておく必要がある。本実施形態では、一例として、観測分布Ｆが正規分布であり、かつ、パラメータの事前分布Ｈが共役事前分布であるようなＮｏｒｍａｌ−Ｇａｍｍａ分布を用いることにした。なお、Ｎｏｒｍａｌ−Ｇａｍｍａ分布については、参考文献１「C. M. ビショップ、“パターン認識と機械学習”、シュプリンガー・ジャパン、2007.」に記載されている。
【００６５】
本実施形態でＮｏｒｍａｌ−Ｇａｍｍａ分布を採用した理由は、これらの分布を利用したｄＨＤＰの近似モデル解法は、従来発表されていない上、多くの分野に応用可能なモデルであり、最も実用性が高いと考えたからである。ただし、目的や実際のデータに則して、これらを別の分布にすることも可能である。
【００６６】
≪ｄＨＤＰ近似モデルのダイアライゼーションシステムへの適用≫
図１のダイアライゼーションシステム１では、式（９）〜式（１５）に示すｄＨＤＰ近似モデルにおいて、ダイアライゼーションの話者数（＝クラスタ数）と各サンプルのクラスタリングを行い、各クラスタのパラメータ（＝話者の位置）を推定することとした。
【００６７】
このうち、各サンプルのクラスタリングは、式（１４）に示す隠れ変数ｚt iの分布ｐ（ｚt i）を求めることに等しい。隠れ変数ｚt iは、各サンプルｘt iに対してどのクラスタに属するかを表す変数なので、このクラスタリング結果が得られれば、各クラスタへ帰属するサンプル数（あるいはその期待値）が計算できる。また、この各クラスタへ帰属するサンプル数（あるいはその期待値）を用いて各クラスタｋの混合比（後記するβk＾）を計算可能である。さらに、この混合比（後記するβk＾）により、「有効な」クラスタ数Ｋeff（＝話者数Ｋeff）を決定することができる。また、クラスタリング結果が求まれば、各クラスタのパラメータ｛θk｝も容易に求めることができる。なお、本実施形態では、クラスタリング計算装置３によって、「有効な」クラスタ数Ｋeff（＝話者数Ｋeff）の決定まで行うこととするが、識別部４がこの処理を行ってもよい。つまり、クラスタリング計算装置３によって、クラスタリング結果を得て、各クラスタのパラメータ｛θk｝を得て、識別部４によって、混合比（後記するβk＾）を求め、「有効な」クラスタ数Ｋeff（＝話者数Ｋeff）を決定するようにしてもよい。
【００６８】
［計算アルゴリズム］
計算アルゴリズムについては、以下、１）ｄＨＤＰモデルのオンライン推定法、２）変分事後分布推論プロセス、３）観測モデルおよび事前分布の特定、４）推定結果とクラスタ数の決定方法、５）ｄＨＤＰの高速化方法に分けて詳細に説明する。
【００６９】
＜１）ｄＨＤＰモデルのオンライン推定法＞
ここでは、具体的な推論アルゴリズムを示す。図２は、本発明の実施形態に係るクラスタリング計算方法の全体処理の流れを示すフローチャートである。図２は、クラスタリングの全体の推論プロセスを示したものである。ダイアライゼーションは、一般的にオンライン計算のプロセスである。しかしながら、ｄＨＤＰモデルのオンライン推定方法は、従来研究されていない。本実施形態では、ｄＨＤＰモデルのオンライン推定法を開発した。図２の符号２０１で示す推論プロセスは、そのオンライン推定法を示したものである。
【００７０】
オンライン推定法は、時刻Ｔ（１≦Ｔ≦Ｔtotal）ごとに行うものである。ここで、刻々と進展する推定対象のある時点の時刻をＴ、それらのうち最終時刻をＴtotalとした。なお、以下では、過去を含めた演算対象時刻をｔとする。過去にどこまで遡及するかについては予め定めておく。例えば、推定対象時刻Ｔ＝５のとき、演算対象時刻ｔ＝１，２，３，４，５としたり、ｔ＝３，４，５としたりすることができる。この例では、ｔ＝１まで考慮することとする。
【００７１】
この推論プロセスでは、処理を開始すると、まず、Ｔを初期化する。すなわち、Ｔ＝１とする（ステップＳ１）。そして、時刻１からＴまでの隠れ変数およびパラメータ（隠れ変数・パラメータ）の推定値およびハイパーパラメータを入力する（ステップＳ２）。また、時刻１からＴ−１までの観測量（サンプル）を入力する（ステップＳ３）。また、時刻Ｔの音声パワーｆTを入力する（ステップＳ４）。この時刻Ｔの音声パワーｆTから、時刻Ｔの観測量｛ｘT i｝を生成して入力する（ステップＳ５）。また、時刻Ｔに対応する未知数（未知の隠れ変数・パラメータ）を初期化する（ステップＳ６）。
【００７２】
そして、隠れ変数・パラメータの変分事後分布を推定する（ステップＳ７）。推定後、Ｔをインクリメントする。すなわち、Ｔ←Ｔ＋１とする（ステップＳ８）。そして、入力が終了したか否かを判別する（ステップＳ９）。入力が終了していない場合（ステップＳ９：Ｎｏ）、ステップＳ２に戻る。一方、例えば、最終時刻Ｔtotalを超えたときに、入力が終了したと判定し（ステップＳ９：Ｙｅｓ）、推定結果を出力する（ステップＳ１０）。なお、ステップＳ２〜ステップＳ４の処理順序は、任意であり、並列処理してもよい。
また、ステップＳ６の処理は、ステップＳ７の前に行うのであれば処理順序は問わない。
【００７３】
オンライン推定法を繰り返すことで、時刻ステップの進展と共に観測サンプル｛ｘt i｝が蓄積され、その都度変数を再推定する。時刻Ｔ−１の時点での推定結果は、次の時刻Ｔでの推定の初期値として利用される。
【００７４】
推論の目的は、全観測データが与えられた時に全未知変数（｛ｚt i｝，｛ｄt i｝，ｗ，｛πt｝，｛θk｝）の推定値を求めることである。これは、確率モデルの観点からは、全変数の事後分布を求めることに相当する。本実施形態では、ｄＨＤＰ近似モデルに対して変分ベイズ法による事後分布推定法（変分事後分布推定法）を示す。実際に求めたい事後分布は、式（１７）であるが、変分法では、式（１８）のように全変数が独立であると仮定した分布（変分事後分布）を推定する。
【００７５】
【数６】
【００７６】
式（１８）に示すｑ（・）は変分事後分布を表す。ここで、時刻Ｔまでの観測量が得られたとする。このとき、現在の時刻Ｔからの過去を含む時刻ｔ（≦Ｔ）に関する変数の変分事後分布推定値ｑ＊は、式（１９）〜式（２３）のように表される。これら各変数の分布は、元の分布（式（９），式（１０），式（１１），式（１３），式（１４））からデータの情報が加わった分修正された形になる。ただし、式（１９）は、観測モデルＦやパラメータの事前分布Ｈに依存する。したがって、本実施形態では、式（１９）に基づいて後記する式（３６）を用いる。
【００７７】
【数７】
【数８】
【００７８】
ただし、式（１９）〜式（２３）を同時に最適化することは困難なので、変分法においてはＥＭアルゴリズム（Expectation−Maximization algorithm）という様式に従った繰り返し計算による推定で各変数を個別に最適化する。ＥＭアルゴリズムとは、複数の変数を同時最大化する計算手法であり、Ｅステップ（Expectation step）とＭステップ（Maximization step）からなる計算ステップを相互に繰り返し計算することで全体最適化を行うものである。なお、ＥＭアルゴリズムについては、前記した参考文献１に記載されている。
【００７９】
＜２）変分事後分布推論プロセス＞
図３は、図２に示す変分事後分布推論手順を示すフローチャートであって、ＥＭアルゴリズムを含んだ推論プロセスを示す。図３に示す推論プロセスでは、まず、時刻１からＴまでのサンプル｛ｘ1:T｝を入力する（ステップＳ２１）。また、時刻１からＴまでの隠れ変数・パラメータの推定値およびハイパーパラメータを入力する（ステップＳ２２）。なお、ステップＳ２１，Ｓ２２の処理順序は任意である。
【００８０】
そして、ＥＭステップの繰り返し回数を表す識別子ｊを初期化する。すなわち、ｊ＝１とする（ステップＳ２３）。そして、Ｅステップを計算する（ステップＳ２４）。Ｅステップでは、時刻１からＴまでのパラメータの推定値を更新する。続いて、Ｍステップを計算する（ステップＳ２５）Ｍステップでは、時刻１からＴまでの隠れ変数の推定値を更新する。
【００８１】
そして、ＥＭステップの繰り返し回数ｊをインクリメントする。すなわち、ｊ←ｊ＋１とする（ステップＳ２６）。そして、現在の繰り返し回数ｊが事前に設定されたしきい値（ｊmax）を超えたか否かを判別する。すなわち、ｊ＞ｊmaxとなったか否かを判別する（ステップＳ２７）。ｊ≦ｊmaxである場合（ステップＳ２７：Ｎｏ）、ステップＳ２４に戻る。一方、ｊ＞ｊmaxとなった場合（ステップＳ２７：Ｙｅｓ）、時刻Ｔが事前に設定された設定値（適切な正整数ｔupdateの倍数）を超えたか否かを判別する。すなわち、Ｔが「ｔupdateの倍数」となったか否かを判別する（ステップＳ２８）。
【００８２】
Ｔが「ｔupdateの倍数」となった場合（ステップＳ２８：Ｙｅｓ）、ハイパーパラメータを更新し（ステップＳ２９）、推定結果を保存する（ステップＳ３０）。Ｔが「ｔupdateの倍数」ではない場合（ステップＳ２８：Ｎｏ）、ステップＳ２９をスキップして、推定結果を保存する（ステップＳ３０）。
【００８３】
このようにＥステップ、Ｍステップをｊmax回に渡って相互に繰り返し計算することで変分事後分布の式（１９）〜式（２３）を得ることができる。なお、これらは、Ｍステップにおいて演算される。
【００８４】
本実施形態では、図３に示すようにステップＳ２８〜Ｓ３０の処理を行うことで、時間ステップＴを適切な正整数ｔupdate回インクリメントするごとにハイパーパラメータを推定することとした。通常、ハイパーパラメータは事前に与える定数であるが、ステップＳ２８〜Ｓ３０の処理を実行することで、ハイパーパラメータが固定値である場合よりも精度を上げることができる。ハイパーパラメータの更新方法については、事後分布からのサンプリング方法など、様々な公知の手法を用いることができる。また、ハイパーパラメータ自身もＥＭアルゴリズムで推定することも可能である。
【００８５】
≪Ｅステップ≫
Ｅステップの具体的な計算式は、以下の式（２４）〜式（２８）で表される。ここで、ψ（・）はプサイ関数（あるいはディガンマ関数）である。また、式（２７）は、観測モデルＦに依存する。したがって、本実施形態では、式（２７）に基づいて後記する式（３７）を用いる。
【００８６】
【数９】
【００８７】
Ｅステップでは、未知パラメータ（θ、π、ｗ）に関わる推定値を再計算する。ここでＥx［ｆ（ｘ）］は、変分分布上の期待値を示す。式（２８）にその定義を示す。
【００８８】
【数１０】
【００８９】
図４は、図３に示すＥステップの計算手順の一例を示すフローチャートである。Ｅステップでは、時刻Ｔまでに関係する全変数について再計算を行う。具体的には、Ｅステップでは、まず、時刻１からＴまでのサンプル｛ｘ1:T｝を入力する（ステップＳ３１）。また、時刻１からＴまでの隠れ変数・パラメータの推定値とハイパーパラメータ、最新のＭステップの演算結果を入力する（ステップＳ３２）。なお、ステップＳ３１，Ｓ３２の処理順序は任意である。
【００９０】
そして、過去を含む演算対象の時刻ｔを初期化する。すなわち、ｔ＝１とする（ステップＳ３３）。そして、この演算対象の時刻ｔに対して、式（２４）、式（２５）を計算する（ステップＳ３４）。次いで、クラスタの識別子ｋを初期化する。すなわち、ｋ＝１とする（ステップＳ３５）。そして、この演算対象の時刻ｔおよびクラスタｋに対して、式（２６）を計算する（ステップＳ３６）。さらに、時刻ｔにおけるデータの識別子ｉを初期化する。すなわち、ｉ＝１とする（ステップＳ３７）。そして、時刻ｔ、クラスタｋのｉ番目のデータに対して、式（２７）を計算する（ステップＳ３８）。
【００９１】
そして、時刻ｔにおけるデータの識別子ｉをインクリメントする。すなわち、ｉ←ｉ＋１とする（ステップＳ３９）。続いて、ｉ＞ｎtとなったか否かを判別する（ステップＳ４０）。なお、ｎtは、式（８）に示した個数である。ｉ≦ｎtである場合（ステップＳ４０：Ｎｏ）、ステップＳ３８に戻る。一方、ｉ＞ｎtとなった場合（ステップＳ４０：Ｙｅｓ）、次のクラスタに更新する。すなわち、ｋ←ｋ＋１とする（ステップＳ４１）。そして、すべてのクラスタについて演算したか否かを判別する。すなわち、ｋ＞Ｋとなったか否かを判別する（ステップＳ４２）。ｋ≦Ｋである場合（ステップＳ４２：Ｎｏ）、ステップＳ３６に戻る。一方、ｋ＞Ｋとなった場合（ステップＳ４２：Ｙｅｓ）、演算対象の時刻ｔを更新する。すなわち、ｔ←ｔ＋１とする（ステップＳ４１）。そして、演算対象の時刻ｔが推定対象の時刻Ｔとなったか否かを判別する。すなわち、ｔ＞Ｔとなったか否かを判別する（ステップＳ４４）。ｔ≦Ｔである場合（ステップＳ４４：Ｎｏ）、ステップＳ３４に戻る。一方、ｔ＞Ｔとなった場合（ステップＳ４４：Ｙｅｓ）、推定対象の時刻ＴにおけるＥステップの推定結果を保存する（ステップＳ４５）。
【００９２】
≪Ｍステップ≫
Ｍテップの具体的な計算式は、以下の式（２９）〜式（３２）で表される。Ｍステップでは、隠れ変数（ｚt，ｄt）に関わる推定値を再計算する。
【００９３】
【数１１】
【００９４】
図５は、図３に示すＭステップの計算手順の一例を示すフローチャートである。Ｍステップでも、Ｅステップと同様に、時刻Ｔまでに関係する全変数について再計算を行う。具体的には、Ｍステップでは、まず、時刻１からＴまでのサンプル｛ｘ1:T｝を入力する（ステップＳ５１）。また、時刻１からＴまでの隠れ変数・パラメータの推定値とハイパーパラメータ、最新のＥステップの演算結果を入力する（ステップＳ５２）。なお、ステップＳ５１，Ｓ５２の処理順序は任意である。
【００９５】
以下、Ｅステップと同様な処理については、説明を適宜省略する。Ｍステップでは、まず、ｔ＝１（ステップＳ５３）、ｉ＝１（ステップＳ５４）、ｋ＝１（ステップＳ５５）とした上で、前記した式（２９）、式（３０）を計算する（ステップＳ５６）。その後、ｋ←ｋ＋１として（ステップＳ５７）、ｋ＞ＫとなるまでステップＳ５６に戻る処理を繰り返す。ｋ＞Ｋとなった場合（ステップＳ５８：Ｙｅｓ）、時刻の識別子ｌを初期化する。すなわち、ｌ＝１とする（ステップＳ５９）。そして、時刻ｔ、時刻ｌ、クラスタｋのｉ番目のデータに対して、式（３１）、式（３２）を計算する（ステップＳ６０）。
【００９６】
そして、時刻ｌをインクリメントする。すなわち、ｌ←ｌ＋１とする（ステップＳ６１）。過去から時刻ｔまでの時刻を示す時刻ｌについてのすべての演算が終わるまでステップＳ６０に戻る処理を繰り返す。ｌ＞ｔとなった場合（ステップＳ６２：Ｙｅｓ）、ｉ←ｉ＋１とする（ステップＳ６３）。そして、ｉについてのすべての演算が終わるまでステップＳ５５に戻る処理を繰り返す。そして、ｉ＞ｎtとなった場合（ステップＳ６４：Ｙｅｓ）、ｔ←ｔ＋１とする（ステップＳ６５）。さらに、ｔについてのすべての演算が終わるまでステップＳ５４に戻る処理を繰り返す。そして、ｔ＞Ｔとなった場合（ステップＳ６６：Ｙｅｓ）、推定結果を保存する（ステップＳ６７）。
【００９７】
＜３）観測モデルおよび事前分布の特定＞
図３〜図５および式（１９）〜式（２７）を参照して説明した変分事後分布推論プロセスでは、観測モデルおよび事前分布を一般化した説明を行った。つまり、前記した式（１９）および式（２７）は、観測モデルＦ＝ｐ（ｘt i|θk）およびパラメータθkの事前分布Ｈに依存する。本実施形態では、観測モデルＦを正規分布と仮定し、事前分布ＨとしてＮｏｒｍａｌ−Ｇａｍｍａ分布を仮定する。Ｎｏｒｍａｌ−Ｇａｍｍａ分布は、式（３３）〜式（３５）で表されるモデルである。
【００９８】
【数１２】
【００９９】
このモデルの場合において、前記した式（１９）および式（２７）は、式（３６）、式（３７）でそれぞれ表される。また、式（３６）、式（３７）中のハイパーパラメータは、式（３８）〜式（４１）で表される。さらに、式（３８）〜式（４１）中の変数は、式（４２ａ）〜式（４２ｃ）で表される。式（３８）〜式（４２）の演算は、Ｅステップの中で完了する。
【０１００】
【数１３】
【０１０１】
【数１４】
【０１０２】
【数１５】
【０１０３】
＜４）推定結果とクラスタ数の決定方法＞
≪クラスタ数の決定方法≫
ＥＭアルゴリズムでは、Ｋ個のクラスタを常に保持しているが、推定が進むと少数のクラスタのみに大きな混合比が付き、他のクラスタの大きさ（cluster size）はほぼ０となる。時刻ｔにおいてクラスタｋへ振り分けられるサンプル数の期待値は、式（４３）の定義により計算できる。
【０１０４】
【数１６】
【０１０５】
本実施形態では、「有効な」クラスタ数Ｋeffを、式（４３）で定義した‖ｚt,k‖の時刻ｔに関する総和の比率で決定する。例えば、各クラスタｋの混合比は、式（４４ａ）で推定できる。このルールでは、式（４４ｂ）の条件が成立するならば、そのクラスタｋが「有効な」クラスタであると判断することとする。
【０１０６】
【数１７】
【０１０７】
このルールにおいて、そのようなクラスタｋの総数を、「有効な」クラスタ数Ｋeffとすることとする。このようにすれば、「有効な」クラスタ数が最大クラスタ数Ｋより小さくなること、すなわち、Ｋeff＜Ｋが保証される。
【０１０８】
≪保存すべき推定結果≫
保存すべき推定結果としては、第１に、ＥＭアルゴリズムで推定された各変数によって定められる変分事後分布である。また、第２に、クラスタリング結果を用いて求められる、「有効な」クラスタ数Ｋeff、式（４３）に示す‖ｚt,k‖、式（４４ａ）に示すβk＾などである。ここで、＾は、文字βの上に付される記号を意味する。
特に、保存すべき推定結果として２番目に挙げたクラスタリング結果を用いて求められる推定量“「有効な」クラスタ数Ｋeff、式（４３）に示す‖ｚt,k‖、式（４４ａ）に示すβk＾”は、図１における識別部４で利用することが可能な重要な推定量である。なお、式（４３）、式（４４ａ）、式（４４ｂ）の演算は、Ｍステップの中で行う。
【０１０９】
＜５）ｄＨＤＰの高速化方法＞
本実施形態のクラスタリング計算方法におけるオンライン推定方法では、図４や図５の処理の流れからも明らかなように、時間ステップの進展とともに推定すべき変数の数が増えていくことが特徴である。そのため、実時間性などを考慮して、計算省力化の技法を開発した。本実施形態では、ｄＨＤＰを用いたオンライン話者クラスタリングのための計算高速化法として、大別して３種類の高速化法（高速化１、高速化２、高速化３）を用いることとした。
【０１１０】
≪高速化１：データの忘却≫
前記した式（１６ａ）が意味することは、時刻ｔにおける話者分布Ｇtを計算するためには、時刻１≦ｌ≦ｔの分布変化Ｈｌが必要であるということである。従って、時刻１からの情報を保持し続ける必要があるので、時間ステップｔが進むとともに推論の計算量が増大する。ここで、次の仮定を導入する。すなわち、時刻ｌ＜ｔで話者の入れ替わりが発生したとする。すると、この時点で話者分布Ｈが大きく変化することが予想される。これは、前記した式（４）において、ｗl≒１となってＧl-1の影響がほぼ排除されることを意味する。このことと、前記した式（１６ａ）および式（１６ｂ）とを照らし合わせると、ｖt1≒…≒ｖl(t-1)＝０となる。従って、実際には、Ｇtの推論には、時刻ｌの分布変化Ｈl以降の分布のみしか関与しないことがわかる。このことから、分布選択を表す変数｛ｄt i｝の事後分布に相当するｓt i l（式（３０）参照）も大部分が０となることが期待される。
【０１１１】
ＥＭアルゴリズム中でｓが存在する部分は、ｓと、定数あるいは別変数と、の掛け算となっているため、ｓt i l＝０となる部分については、計算する必要がない。従って、ｓに対応する変数にはアクセスしなくてよい。この考察より、ある適切な時間ステップ幅Ｗ１より過去の変数あるいは定数については、ＥＭ更新式中でアクセスしない（忘却する）ことによって計算時間を削減できる。
【０１１２】
この時間ステップ幅Ｗ１の決定方法としては、事前に適切な定数を決めておく方法が最も簡便である。あるいは、この考察からの示唆より、次のような方法も考えることができる。すなわち、各時刻でのｓt iの推定結果より、ｓt i 1＝…＝ｓt i l＝０となるｌを見つけて、Ｗ１=ｔ−ｌとする方法がある。そのようなｌの見付け方としては、適切な閾値ｔｈ（＜１．０）に対して、式（４５）の関係を満たす最大のｌを探せばよい。この場合には、データを忘却する時間幅Ｗ１が推論結果に従って動的に変化することになる。
【０１１３】
【数１８】
【０１１４】
≪高速化２：推定時間幅の制限≫
図２に示したオンライン推定法では、毎時刻Ｔにおいて全変数のＥＭアルゴリズム推定を行うものとして説明した。このことは、早い時刻ステップに関する変数については、何度もＥＭ再推定を行っているということを意味する。ＥＭアルゴリズムは、局所最適解の計算アルゴリズムであるため、計算を繰り返すと、ある１つの解に収束する。よって、早い時刻ステップに関わる変数については再推定をしなくとも、値が収束している可能性が高い。
【０１１５】
この考察より、ＥＭアルゴリズムで再推定する変数を、適切な時間幅Ｗ２を使うことで、Ｔ−Ｗ２≦ｔ≦Ｔの範囲の変数に制限することが可能であることがわかる。このことで、推定する変数の個数が、時間ステップＴに対する線形増加とならないように、その個数を定数に保つことが可能となる。時間幅Ｗ２の決定方法としては、事前に適切な定数を決めておく方法が最も簡便である。
【０１１６】
≪高速化３：クラスタ変数再推定の制限≫
ｄＨＤＰは、ノンパラメトリックベイズモデルなので、最大クラスタ数（Ｋ）の個数のクラスタを予めモデルに準備していたとしても、実際には、「有効な」クラスタ数Ｋeff（＜Ｋ）の個数のクラスタだけが実質的なクラスタとなり、他のクラスタは、混合比がほぼ０となるために消去されたような形式となる。この消去されたような形式のクラスタには、有効な情報がないため、そのようなクラスタを用いてパラメータや混合比の推定計算を行っても無駄なだけである。
【０１１７】
このことから、前記した式（２６）、式（２７）、式（２９）などクラスタごとの推定を行う変数（ｋ）を用いて、実際にクラスタごとの変数再推定を行う演算処理では、その変数再推定の回数を確率的に（確率論に基づいて）低減することが考えられる。その方法は、例えば３種類挙げることができる。
【０１１８】
第１の低減法は、単純な方法として、全クラスタについて確率ｃ（≦１．０）で毎回ランダムに、再推定するか否かを決定する方法である。
第２の低減法は、各クラスタｋの混合比βk＾（式（４４）参照）に応じてＥＭアルゴリズムによる再推定を行う回数を増減する方法である。
第３の低減法は、前記した式（４４ｂ）の条件を満たす場合、つまり、対象のクラスタが「有効な」クラスタである場合には毎回再推定するが、それ以外の場合には、クラスタｋを確率的に更新する方法である。
これらの低減法において、各クラスタｋを更新する確率をｐupdate（ｋ）とおくと、第１〜第３の低減法は、それぞれ、式（４６）〜式（４８）のように表現できる。特に、第３の低減法は、推定精度をほとんど犠牲にせずに計算量をＫeff／Ｋ程度に削減できる効果がある。
【０１１９】
【数１９】
【０１２０】
前記した３つの高速化方法（高速化１、高速化２、高速化３）をすべて使う場合のＥＭアルゴリズムの一例を図６および図７にそれぞれ示す。なお、これら３つの方法は、それぞれ独立に利用可能である。
【０１２１】
図６は、本発明の実施形態に係るクラスタリング計算方法におけるＥステップの計算手順を示すフローチャートである。なお、図４のフローチャートと比較して異なる部分を太線及び破線で示し、説明を適宜省略して異なる処理について説明する。図６において、高速化１は、符号３０１で示すように、ステップＳ３１Ａ，Ｓ３２Ａからなる。ここでは、時刻Ｔ−Ｗ1からＴまでのサンプル｛ｘT-W1:T｝を入力する（ステップＳ３１Ａ）。また、時刻Ｔ−Ｗ1からＴまでの隠れ変数・パラメータの推定値とハイパーパラメータ、最新のＭステップの演算結果を入力する（ステップＳ３２Ａ）。次いで、高速化２は、符号３０２で示すように、ステップＳ３３Ａからなる。ここでは、ｔの初期値を１の代わりに、ｔ＝Ｔ−Ｗ2とする（ステップＳ３３Ａ）。
【０１２２】
さらに、高速化３は、ステップＳ３５に続いて行う処理であり、符号３０３で示すように、ステップＳ７１，Ｓ７２からなる。ここでは、［０，１］の一様乱数ｕを生成し（ステップＳ７１）、ｐupdate（ｋ）＞ｕであるか否かを判別することとした（ステップＳ７２）。ｐupdate（ｋ）＞ｕである場合（ステップＳ７２：Ｙｅｓ）、前記したステップＳ３６〜ステップＳ４５を実行する。つまり、前記した式（２６）および式（２７）の計算を行う。ただし、ステップＳ４２において、ｋ≦Ｋである場合（ステップＳ４２：Ｎｏ）、ステップＳ７１に戻る。一方、ステップＳ７２において、ｐupdate（ｋ）≦ｕである場合（ステップＳ７２：Ｎｏ）、前記したステップＳ３６〜ステップＳ４０をスキップして、ステップＳ４１に進む。つまり、前記した式（２６）および式（２７）の計算を行わない。
【０１２３】
図７は、本発明の実施形態に係るクラスタリング計算方法におけるＭステップの計算手順を示すフローチャートである。なお、図５のフローチャートと比較して異なる部分を太線及び破線で示し、説明を適宜省略して異なる処理について説明する。図７において、高速化１は、符号４０１で示すように、ステップＳ５１Ａ，Ｓ５２Ａからなる。ここでは、時刻Ｔ−Ｗ1からＴまでのサンプル｛ｘT-W1:T｝を入力する（ステップＳ５１Ａ）。また、時刻Ｔ−Ｗ1からＴまでの隠れ変数・パラメータの推定値とハイパーパラメータ、最新のＥステップの演算結果を入力する（ステップＳ５２Ａ）。次いで、高速化２は、符号４０２で示すように、ステップＳ５３Ａからなる。ここでは、ｔの初期値を１の代わりに、ｔ＝Ｔ−Ｗ2とする（ステップＳ５３Ａ）。
【０１２４】
さらに、高速化３は、ステップＳ５５に続いて行う処理であり、符号４０３で示すように、ステップＳ８１，Ｓ８２からなる。ここでは、［０，１］の一様乱数ｕを生成し（ステップＳ８１）、ｐupdate（ｋ）＞ｕであるか否かを判別することとした（ステップＳ８２）。ｐupdate（ｋ）＞ｕである場合（ステップＳ８２：Ｙｅｓ）、前記したステップＳ５６〜ステップＳ６７を実行する。つまり、前記した式（２９）〜式（３２）の計算を行う。ただし、ステップＳ５８において、ｋ≦Ｋである場合（ステップＳ５８：Ｎｏ）、ステップＳ８１に戻る。一方、ステップＳ８２において、ｐupdate（ｋ）≦ｕである場合（ステップＳ８２：Ｎｏ）、前記したステップＳ５６をスキップして、ステップＳ５７に進む。つまり、前記した式（２９）および式（３０）の計算を行わない。
【０１２５】
［クラスタリング計算装置］
図８は、本発明の実施形態に係るクラスタリング計算装置の構成の一例を示す機能ブロック図である。このクラスタリング計算装置３は、例えば、図２の推論プロセス等を実現するものであり、ＣＰＵ等の演算装置と、メモリ、ハードディスク等の記憶装置（記憶手段）と、マウスやキーボード等の外部から情報の入力を検出する入力装置と、外部との各種情報の送受信を行うインタフェース装置と、ＬＣＤ（Liquid Crystal Display）等の表示装置を備えたコンピュータと、このコンピュータにインストールされたプログラムとから構成される。
【０１２６】
クラスタリング計算装置３は、ハードウェア装置とソフトウェアとが協働することによって、前記したハードウェア資源がプログラムによって制御されることにより実現され、図８に示すように、記憶手段１０と、以下ＣＰＵの機能として、変分事後分布推論部３０と、パワーベクトル読込部２１と、入力制御部２２と、割当部２３と、未知数初期化部２４と、パワーベクトル書込部２５と、観測量生成部２６と、終了判定部２７と、出力制御部２８とを備えることとした。
【０１２７】
＜記憶手段＞
記憶手段１０は、ＲＯＭ、ＲＡＭ、ＨＤＤ等からなる。記憶手段１０は、プログラム記憶領域、設定データ記憶領域、演算データ一時記憶領域、推定結果記憶領域等に区分され、コマンド、データ、プログラム等各種情報を記憶する。例えば、推定結果記憶領域には、データとして、図９に示すように、未知数の事後分布推定値１１、ハイパーパラメータ１２、Ｅステップにおける推定値（まとめてＥステップ１３と表記した）、Ｍステップにおける推定値（まとめてＭステップ１４と表記した）、パワーベクトル１５、観測量（サンプル）１６が保存される。
【０１２８】
＜変分事後分布推論部＞
メインの変分事後分布推論部３０は、例えば、図３、図６、図７の推論プロセスを実現するものであり、ここでは、Ｅステップ用計算部３１と、Ｍステップ用計算部３２と、ＥＭ収束判定部３３と、パラメータ更新条件判定部３４と、ハイパーパラメータ更新部３５とを備えることとした。詳細は後記する。
【０１２９】
＜パワーベクトル読込部＞
パワーベクトル読込部２１は、推定対象の時刻Ｔにしたがって、パワーベクトルｆTを順次読み込み、入力制御部２２に渡す。
【０１３０】
＜入力制御部＞
入力制御部２２は、パワーベクトルｆTを取得すると、パワーベクトル書込部２５にそのまま送ると共に、割当部２３、未知数初期化部２４、観測量生成部２６および終了判定部２７に対しては、それぞれ必要なコマンドを出力する。
【０１３１】
＜割当部＞
割当部２３は、入力制御部２２からのコマンドにしたがって、そのときのＴ（＝１，２，…）や、角度ｄ（１８０，−１７９，…）に対応したｉに応じて、前記した式（９）〜（１４）の各未知数（パラメータ・隠れ変数）を順次生成し（初期値は例えば０）、記憶手段１０において、未知数の事後分布推定値１１として保存する。すなわち、割当部２３は、推定すべき未知数にＴとｉとを割り当てる。この意味で、図９においては、Ｔとｉの割当部２３と表記した。なお、記憶手段１０には、ｄＨＤＰ近似モデルで用いる各ハイパーパラメータ１２の初期値が予め格納されている。
【０１３２】
＜未知数初期化部＞
未知数初期化部２４は、入力制御部２２からのコマンドにしたがって、そのときのＴや角度ｄに対応したｉに応じて、Ｅステップ演算に用いるための、前記した式（２４）〜（２６）、（３７）の左辺パラメータを順次生成し（初期値は例えば０）、記憶手段１０において、Ｅステップ１３の初期値として保存する。また、未知数初期化部２４は、Ｍステップ演算に用いるための、式（２９）〜（３２）の左辺パラメータを順次生成し（初期値は例えば０）、記憶手段１０において、Ｍステップ１４の初期値として保存する。さらに、未知数初期化部２４は、ＥＭ推定値の演算に用いるための、式（３６），（２０）〜（２３）の左辺パラメータを順次生成し（初期値は例えば０）、記憶手段１０において、未知数の事後分布推定値１１に上書き保存する（更新する）。
【０１３３】
＜パワーベクトル書込部＞
パワーベクトル書込部２５は、入力制御部２２から取得したパワーベクトルｆTを、記憶手段１０において、パワーベクトル１５として順次保存する。
【０１３４】
＜観測量生成部＞
観測量生成部２６は、入力制御部２２からのコマンドにしたがって、そのときのＴに応じて、記憶手段１０から読み出したパワーベクトル１５のデータの角度ｄ毎のパワー値を、所定のルールに則って変換し、識別子ｉ（個数ｎｔ）のデータに変換することで、観測量（サンプル）を生成し、記憶手段１０に、観測量１６として順次保存する。なお、本実施形態では、所定のルールとしては、前記した式（７）、式（８）を用いる。
【０１３５】
＜終了判定部＞
終了判定部２７は、入力制御部２２からの入力信号（コマンド）が所定期間途絶えたときに、パワーベクトルｆTの入力が終了したと判定し、出力制御部２８に通知する。なお、本実施形態では、最終時刻Ｔtotalになったときに、入力が終了したと判定する。
【０１３６】
＜出力制御部＞
出力制御部２８は、終了通知を受け取ると、記憶手段１０から、最終的な推定値として、未知数の事後分布推定値１１を取得し、識別部４に出力する。
【０１３７】
［変分事後分布推論部の詳細］
≪Ｅステップ用計算部≫
Ｅステップ用計算部３１は、処理時点のＴに応じて、記憶手段１０から、サンプル、ハイパーパラメータ、Ｍステップの演算結果等を、過去を含めた演算対象時刻ｔ（ｔ≦Ｔ）に関して読み出す。そして、Ｅステップ用計算部３１は、処理時点のＴに応じて、過去を含むすべての演算対象時刻ｔ（ｔ≦Ｔ）に関して、式（２４）〜（２６）、（３７）を演算し、その演算結果を、記憶手段１０にＥステップ１３として保存する。
【０１３８】
≪Ｍステップ用計算部≫
Ｍステップ用計算部３２は、処理時点のＴに応じて、記憶手段１０から、サンプル、ハイパーパラメータ、Ｅステップの演算結果等を、過去を含めた演算対象時刻ｔ（ｔ≦Ｔ）に関して読み出す。そして、Ｍステップ用計算部３２は、処理時点のＴに応じて、過去を含むすべての演算対象時刻ｔ（ｔ≦Ｔ）に関して、式（２９）〜（３２）を演算し、その演算結果を、記憶手段１０にＭステップ１４として保存する。そして、Ｍステップ用計算部３２は、記憶手段１０に保存されたこれらの演算結果を合わせて用いて、処理時点のＴに応じて、過去を含むすべての演算対象時刻ｔ（１≦ｔ≦Ｔ）に関して、式（２０）〜（２３）、（３６）を演算し、その演算結果を、記憶手段１０に、未知数の事後分布推定値１１として上書き保存する（更新する）。
【０１３９】
≪ＥＭ収束判定部≫
ＥＭ収束判定部３３は、ＥステップとＭステップを合わせた１セットの処理の繰り返し回数ｊのしきい値（ｊmax）に達したか否かを判定することで、ＥＭアルゴリズムが収束したか否かを判別する。しきい値（ｊmax）は、予め設定される。ＥＭアルゴリズムが収束していない場合、ＥＭ収束判定部３３は、ＥステップとＭステップとを繰り返す制御を行う。ＥＭアルゴリズムが収束した場合、ＥＭ収束判定部３３は、ＥステップとＭステップとの処理を停止する制御を行う。本実施形態では、ＥＭ収束判定部３３は、収束の如何に関わらず、その時点のＴの値をパラメータ更新条件判定部３４に通知することとした。
【０１４０】
≪パラメータ更新条件判定部≫
パラメータ更新条件判定部３４は、受け取った処理時点のＴが、事前に設定された設定値（適切な正整数ｔupdateの倍数）と同じであるか否かを判別する。同じある場合には、パラメータ更新条件判定部３４は、ハイパーパラメータ更新部３５に通知する。
【０１４１】
≪ハイパーパラメータ更新部≫
ハイパーパラメータ更新部３５は、通知を受けると、記憶手段１０に保存されているハイパーパラメータ１２を任意の方法で更新する。これにより、Ｅステップ用計算部３１およびＭステップ用計算部３２は、処理時点のＴが、予め設定された適切な正整数の倍数と同じである次のタイミングから、更新されたハイパーパラメータを用いることができる。なお、パラメータ更新条件判定部３４およびハイパーパラメータ更新部３５は、前記したステップＳ２８〜Ｓ３０の処理に対応している。ただし、これらの構成は、必ずしも備える必要はない。
【０１４２】
なお、クラスタリング計算装置３は、一般的なコンピュータを、クラスタリング計算装置３を構成する前記した各手段として機能させるプログラム（クラスタリング計算プログラム）により動作させることで実現することができる。このプログラムは、通信回線を介して提供することも可能であるし、ＣＤ−ＲＯＭ等の記録媒体に書き込んで配布することも可能である。このプログラムをインストールされたコンピュータは、ＣＰＵが、ＲＯＭ等に格納されたこのプログラムをＲＡＭに展開することにより、クラスタリング計算装置３と同等の効果を奏することができる。
【０１４３】
本実施形態によれば、複数話者の会話からなる録音データから、その会話に参加した話者数、各話者の位置、さらに各話者の発話行動のタイミングを推定するダイアライゼーションにおける話者クラスタリングの問題に対して、確率的なクラスタリングを用いたので、従来のようなパラメータの設定や探索によらずに容易に話者数を推定できる。
また、本実施形態によれば、ｄＨＤＰ近似モデルを採用したことにより、時間ごとに発話に参加する話者が変化する状況も適切にモデリングできる。その結果、より精度の高い話者クラスタリングが実現できる。
さらに、本実施形態によれば、ｄＨＤＰのオンライン推定法とその高速化法を用いることで、高速に推論を行うことが可能である。なお、精度と時間のトレードオフによって実時間推論も可能である。
【０１４４】
以上、本発明の実施形態について説明したが、本発明はこれに限定されるものではなく、その趣旨を変えない範囲で実施することができる。例えば、クラスタリング計算装置３は、ｄＨＤＰ近似モデルにおいて、Ｎｏｒｍａｌ−Ｇａｍｍａ分布を仮定したが、別の分布を仮定してもよい。この場合、前記した式（３６）、式（３７）および関連する関係式と同様な計算式を、仮定された別の分布で改めて導出すればよい。
【０１４５】
また、本実施形態では、クラスタリング計算装置３は、ｄＨＤＰ近似モデルにより推論するものとしたが、例えば、ｄＨＤＰモデルや他のノンパラメトリックベイズモデルであってもよい。他のノンパラメトリックベイズモデルの場合、確率分布の連続的な時間変化をモデル化した分布であることが好ましい。
【０１４６】
また、本実施形態では、話者の位置に基づいたクラスタリングを行うために、特徴量抽出部２では、一例として、マイクロホンアレーを用いてＤＯＡ情報（音声の到達角度）を抽出するものとしたが、ダイアライゼーションに適していれば、それ以外の様々な特徴量を抽出することが可能である。例えば、話者ごとの声質の特徴に基づいたクラスタリングをする場合にはＭＦＣＣ特徴量（Mel Frequency Cepstrum Coefficient）などを抽出することができる。そして、このような特徴量であっても前記したアルゴリズムに適用することが可能である。
【実施例】
【０１４７】
本発明の効果を確認するために本実施形態に係るクラスタリング計算装置の性能を測定した。まず、第１段階の実験として、設計されたとおりのクラスタ数と適切なパラメータとを持つクラスタリングが、ｄＨＤＰによって実現できるか否かを検証する実験（クラスタリング検証実験）を行った。次に、第２段階の実験として、得られたクラスタ結果を用いたダイアライゼーション精度を評価する実験（ＤＥＲによるダイアライゼーション精度の評価）を行った。前記２種類の実験では、人工音声データと実音声データとを用いてクラスタリング計算装置の性能を確認した。
【０１４８】
＜実験データ＞
≪人工音声データ≫
人工音声データは、３人の発話者が交代で発話・非発話を切り替える状況をシミュレートしたデータである。この人工音声データは、６４［ｍｓｅｃ］の時間ステップで計算されたＤＯＡデータ（音声の到達角度のデータ）と、ＶＡＤ（音声区間検出器）による音声・非音声判定結果とからなる。また、人工音声データは、ノイズがほとんど重畳しないデータである。各実験では、ＶＡＤによる判定結果を用いて非音声区間を閾値処理で除外して４２２ステップの連続シーケンスの人工音声データを作成した。また、時間ステップごとのサンプル分布をある程度安定させるため、４２２ステップを数ステップごとに重複なくまとめて１つの長いメタステップとした。各実験では、５ステップのデータをまとめて１メタステップとし、メタステップのステップ数をｔ，Ｔに対応させた。よって、ステップ数はＴtotal＝ｃｅｉｌ（４２２／５）＝８５である。ここで、ｃｅｉｌは切り上げを示す。この人工音声データには、複数の話者が同時に喋る区間も存在している。
【０１４９】
人工音声データについての全時刻での平均パワーベクトルの分布を図９に示す。図９において、横軸は、ＤＯＡデータ、すなわち、音声到来方向角度［ｄｅｇ］を示し、縦軸は、時間平均パワー［ｄＢ］を示す。図９に示すように、各話者に対応するパワー分布の３つのピークを観察することができる。
【０１５０】
≪実音声データ≫
実音声データは、実際の複数話者の会話の様子を録音したデータである。実音声データとして、非特許文献１に記載されている４データを利用した。４データの詳細を表１に示す。表１において、ＣＰはcrossword puzzleデータ、ＤＣはdiscussionデータ、ＣＮはconversationデータを表す。
【０１５１】
【表１】
【０１５２】
実音声データは、どのデータも３００秒の音声データである。これらの実音声データについての平均パワーベクトルの分布を図１０にそれぞれ示す。図１０の横軸および縦軸は、図９のグラフと同様である。ただし、時間平均パワーのオーダーは低い。
【０１５３】
図１０（ａ）に示すＣＰ１と、図１０（ｂ）に示すＣＰ２に関しては、表１に記載の発話者数「４」に対応するパワー分布のピークの個数も「４個」なので、比較的良好なクラスタリング結果が期待される。一方、図１０（ｃ）に示すＤＣと、図１０（ｄ）に示すＣＮに関しては、表１に記載の発話者数「３」に対応するパワー分布の３つのピークが明確には観察できないため、正しいクラスタ数の推定は困難を伴うことが予想される。
【０１５４】
［クラスタリング検証実験］
まず、第１段階として、クラスタリングの性能を確認した。ここでは、ＤＰＭ（参考例１）とｄＨＤＰ（実施例１）とによって、オンラインで話者クラスタの推定を行った。このクラスタリング検証実験では、最終時刻Ｔtotalにおいて、最終的な混合比が、チャンスレベル（１／Ｋ）を超えたクラスタのみを有効なクラスタとしてカウントすることで、最終的なクラスタリング結果を求めた。そして、得られたクラスタリング結果により、ＤＰＭ（参考例１）とｄＨＤＰ（実施例１）とを比較し、話者の分布と話者数に対応したクラスタリング結果が得られているか否かを判定した。なお、チャンスレベルは偶然の一致が起こる確率である。
【０１５５】
≪人工音声データの場合≫
図９に示す人工音声データに対して、オンラインで、ＤＰＭ（参考例１）とｄＨＤＰ（実施例１）とを適用した結果を図１１に示す。図１１において、横軸は、ＤＯＡデータを、［−１８０：１８０］→［−０．５：０．５］に変換する関数を用いて無次元量に正規化した角度、すなわち、正規化角度を示す。また、縦軸は、確率密度関数（probabilistic density function：p. d. f）の値（無次元数）を示す。
【０１５６】
図１１（ｂ）に示すｄＨＤＰ（実施例１）の結果と、図９のグラフとの比較の結果、ｄＨＤＰでは、正しいクラスタ数とパラメータを得ることができたことが分かる。一方、図１１（ａ）に示すＤＰＭ（参考例１）の場合、クラスタ数が「１」となり、図９に示す人工音声データに対しては不適切な結果が得られたことを確認した。
【０１５７】
≪実音声データ（ＣＰ１、ＣＰ２）の場合≫
図１０（ａ）に示す実音声データ（ＣＰ１）に対して、オンラインで、ＤＰＭ（参考例２）とｄＨＤＰ（実施例２）とを適用した結果を図１２に示す。図１２の各軸は、図１１のグラフと同様である。図１２（ｂ）に示すｄＨＤＰ（実施例２）の結果と、図１０（ａ）のグラフとの比較の結果、ｄＨＤＰでは、正しいクラスタ数とパラメータを得ることができたことが分かる。一方、図１２（ａ）に示すＤＰＭ（参考例２）の場合、多数のクラスタに分かれてしまい、図１０（ａ）に示すＣＰ１に対しては不適切な結果が得られたことを確認した。なお、同様の傾向がＣＰ２データに対しても確認された。このときのＤＰＭ（参考例３）とｄＨＤＰ（実施例３）との結果の記述は省略した。
【０１５８】
≪実音声データ（ＤＣ）の場合≫
図１０（ｃ）に示す実音声データ（ＤＣ）に対して、オンラインで、ＤＰＭ（参考例４）とｄＨＤＰ（実施例４）とを適用した結果を図１３に示す。図１３の各軸は、図１１のグラフと同様である。図１３（ｂ）に示すｄＨＤＰ（実施例４）の結果と、図１０（ｃ）のグラフや表１における発話者数との比較の結果、ｄＨＤＰでは、話者数である３クラスタを得ることはできなかった。しかしながら、各クラスタのうち、サイズ数（図中のcluster size）の上位３つである、「cluster 4」、「cluster 6」、「cluster 14」の正規化角度の位置は、ＤＣデータにおける話者の位置に対応することができた。ここで、サイズ数（図中のcluster size）は、前記した式（４４）の右辺の分母で示される数値で定義した。なお、サイズ数の下位２つは、ノイズクラスタとなった。
【０１５９】
一方、図１３（ａ）に示すＤＰＭ（参考例４）の場合、より多くのクラスタに分かれてしまい、サイズ数上位３つのクラスタも話者の位置にそれぞれ対応することができなかった。この点からも、ｄＨＤＰ（実施例４）は、ＤＰＭ（参考例４）によるクラスタリングに比してより正確なクラスタリングを実現できたと考えられる。なお、ＤＰＭ（参考例４）の場合、各クラスタのうち、サイズ数の上位３つである、「cluster 2」、「cluster 4」、「cluster 5」の正規化角度の位置は、ＤＣデータにおける話者の位置のうちの２つにしか対応できなかった。また、第４〜６位の位置でもあと１つの話者位置に対応できなかった。また、ＣＮデータについては、ＤＣデータと同様な傾向を有すると考えられる。
【０１６０】
［ＤＥＲによるダイアライゼーション精度の評価］
第１段階のクラスタリング検証実験に続いて、第２段階では、ダイアライゼーションのためのクラスタリングとしての性能を評価するため、ＤＥＲ（diarization error ratio）による評価を試みた。ＤＥＲとはＮＩＳＴが提案した話者識別能力の指標である。具体的には、ＤＥＲは、全音声区間長に対して、以下の（１）〜（３）の３種類の誤識別区間がどれだけあったかを百分率で示したものである。ＤＥＲ値が少ないほど良いダイアライゼーションができた、という評価になる。
【０１６１】
（１）false alarm speaker time：誰も話していないのに、誰かが話したと誤検出した区間長
（２）missed speaker time：誰かが話しているのに、誰も話していないと判断した区間長
（３）speaker error time：誰かが話していることは正しく検出したが、話者を誤った区間長
【０１６２】
なお、ＤＥＲについての詳細は、下記ＵＲＬに記載されている。
「NIST Speech group,”Spring2007(RT-07) Rich Transcription Meeting Recognition Evaluation Plan”,[online],[平成２１年１月21日検索]、インターネット＜URL:http://www.nist.gov/speech/tests/rt/2007/index.html＞」
【０１６３】
第１段階のｄＨＤＰによるクラスタリングは、ダイアライゼーションのサブ問題でしかなく、そのままでは話者識別はできない。しかしながら、ｄＨＤＰによるクラスタリングでは、各時刻におけるサンプル（＝方向付きの音声パワーデータ）をクラスタリングしているため、各フレーム（時間ステップ）において各クラスタにアサインされたサンプル数を数えれば、所定の閾値を用いることで話者ごとの発声または非発声を決定することが可能である。そこで、前記した式（４３）で定義した‖ｚt,k‖を用いて、各時刻における話者ｋの発話または非発話を、式（４９）および式（５０）に示すルールで決定した。なお、所定の閾値として式（４９）に示したτＤＥＲは、適切な値に設定した。表１に示した各実音声データに対して、ｄＨＤＰ（実施例５、実施例６）について算出されたＤＥＲの結果を表２に示す。このＤＥＲ値を非特許文献１における結果（比較例）と比較した。
【０１６４】
【数２０】
【０１６５】
【表２】
【０１６６】
表２において、比較例は、既存手法である非特許文献１の中で報告された最良の結果を示す。本発明のｄＨＤＰクラスタリングの結果（実施例５、実施例６）のうち、実施例５（naive）は、式（４９）および式（５０）に示すルールのみでＤＥＲを算出したときの値を示す。また、実施例６（heuristic）は、式（４９）および式（５０）に示すルール（識別則）に加え、１フレーム内（１時間ステップ内）での同時発話人数に上限の仮定をおくことで、非発話区間のfalse alarmを低減した実施例である。実施例６で採用した方法は、非特許文献１で最良の方法が得られた識別則と共通点を持っている。すなわち、実施例６は、式（４９）に示したサンプル数閾値τＤＥＲと、同時発話人数の上限の仮定とに基づいて探索を行った結果である。そのため、実施例６は実施例５よりも良好な結果となった。さらに、実施例６は、同時発話人数の上限の仮定をした比較例と比較しても、明らかに、ＤＥＲ値として良好な値を示すことがわかる。
【０１６７】
以上の検証実験および評価実験をまとめると、ｄＨＤＰで解決されるのはサブ問題（第１段階）であるクラスタリングであるが、第２段階で、クラスタリングの過程で計算されるサンプルアサインメントｚを用いることで、ダイアライゼーションの精度を示すＤＥＲ指標の向上を実現可能である。すなわち、本発明によれば、ｄＨＤＰでクラスタリングの問題を解決したことで、自然と良好なダイアライゼーションが可能となったと結論できる。
【符号の説明】
【０１６８】
１ ダイアライゼーションシステム
２ 特徴量抽出部（特徴量抽出手段）
３ クラスタリング計算装置
４ 識別部（識別手段）
１０ 記憶手段（推定値記憶手段、観測量記憶手段）
２１ パワーベクトル読込部（読込手段）
２２ 入力制御部
２３ 割当部
２４ 未知数初期化部
２５ パワーベクトル書込部
２６ 観測量生成部（観測量生成手段）
２７ 終了判定部
２８ 出力制御部（出力制御手段）
３０ 変分事後分布推論部（事後分布推論手段）
３１ Ｅステップ用計算部（Ｅステップ用計算手段）
３２ Ｍステップ用計算部（Ｍステップ用計算手段）
３３ ＥＭ収束判定部（収束判定手段）
３４ ハイパーパラメータ更新条件判定部
３５ ハイパーパラメータ更新部
【特許請求の範囲】
【請求項１】
話者数が未知である会話の録音データから前記会話の話者数を推定するために、各話者を特徴付ける特徴量を抽出する特徴量抽出手段と、前記抽出された特徴量から前記各話者に対応する複数のクラスタを生成するときの複数の未知パラメータの値をそれぞれ推定するクラスタリング計算装置と、前記推定された複数のパラメータ値により前記会話の各話者を識別する識別手段とを有したダイアライゼーションシステムの前記クラスタリング計算装置であって、
前記抽出された特徴量を読み込む読込手段と、
前記読み込んだ特徴量をノンパラメトリックベイズモデルに対応したベクトルの観測量に変換する観測量生成手段と、
前記変換された観測量の集合データを蓄積記憶する観測量記憶手段と、
前記観測量の集合データから複数のクラスタをノンパラメトリックベイズモデルにより生成するときの複数のパラメータの事後分布の値をＥＭアルゴリズムによりそれぞれ推定および更新する事後分布推論手段と、
前記推定および更新された複数のパラメータの事後分布の値を蓄積記憶する推定値記憶手段と、
事前に設定された終了条件が成立したときに前記推定値記憶手段に記憶されている前記複数のパラメータの事後分布の最新の推定値を出力する出力制御手段と、
を備えることを特徴とするクラスタリング計算装置。
【請求項２】
前記事後分布推論手段は、
前記ＥＭアルゴリズムのＥステップの処理として、ｄＨＤＰ（dynamic Hierarchical Dirichlet Process）モデルにおいて予め定められた事前分布、観測分布および最大クラスタ数と、ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＭステップまでに推定された隠れ変数の事後分布の推定値とを読み込んで、過去分を含めた演算対象時刻ごと、かつ、クラスタごと、かつ、演算対象時刻別の全データごとに、前記クラスタと、混合比と、前記クラスタ分布の時間変化の程度を表す重みとに関するパラメータの事後分布の値を演算するＥステップ用計算手段と、
前記ＥＭアルゴリズムのＭステップの処理として、前記ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＥステップまでに推定されたパラメータの事後分布の推定値とを読み込んで、過去分を含めた演算対象時刻ごと、かつ、演算対象時刻別の全データごとに、２種類の隠れ変数の事後分布の値を推定し、前記２種類の隠れ変数のうち、前記クラスタ分布の時間変化の程度を表す重みに関連付けられた第１隠れ変数の事後分布の値についてはクラスタごとに演算し、前記第１隠れ変数および前記混合比に関連付けられた第２隠れ変数の事後分布の値については、演算対象時刻から過去に遡及した時刻ごとに演算するＭステップ用計算手段と、
前記Ｅステップの処理と前記Ｍステップの処理とを交互に予め定められた回数だけ繰り返し実行させる制御を行う収束判定手段と、
を備えることを特徴とする請求項１に記載のクラスタリング計算装置。
【請求項３】
前記Ｅステップ用計算手段は、
観測量の集合データと、隠れ変数の事後分布の推定値とについて、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込み、推定対象の時刻から事前に設定された設定値分だけ遡及した過去の時刻までを演算対象時刻として、前記クラスタと、前記混合比と、前記重みとに関するパラメータの事後分布の値を演算し、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、事前に設定された再推定すべきクラスタの判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合に、前記混合比に関するパラメータの事後分布の値を演算し、
前記Ｍステップ用計算手段は、
観測量の集合データと、パラメータの事後分布の推定値とについて、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込み、推定対象の時刻から事前に設定された設定値分まで遡及した過去の時刻までを演算対象時刻として、前記第１隠れ変数および第２隠れ変数の事後分布の値を演算し、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、前記判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合に、前記第１隠れ変数の事後分布の値を演算することを特徴とする請求項２に記載のクラスタリング計算装置。
【請求項４】
話者数が未知である会話の録音データから前記会話の話者数を推定するダイアライゼーションシステムにおいて、記憶手段と、読込手段と、観測量生成手段と、事後分布推論手段と、出力制御手段とを備えて、前記録音データから抽出された各話者を特徴付ける特徴量から前記各話者に対応する複数のクラスタを生成するときの複数の未知パラメータの値をそれぞれ推定するクラスタリング計算装置のクラスタリング計算方法であって、
前記読込手段によって、前記抽出された特徴量を読み込む特徴量読込ステップと、
前記観測量生成手段によって、前記読み込んだ特徴量をノンパラメトリックベイズモデルに対応したベクトルの観測量に変換し、前記変換された観測量の集合データを記憶手段に順次蓄積する観測量蓄積ステップと、
前記事後分布推論手段によって、前記観測量の集合データから複数のクラスタをノンパラメトリックベイズモデルにより生成するときの複数のパラメータの事後分布の値をＥＭアルゴリズムによりそれぞれ推定し、当該推定値を前記記憶手段に順次格納および更新する事後分布推定ステップと、
前記出力制御手段によって、事前に設定された終了条件が成立したときに前記記憶手段に記憶されている前記複数のパラメータの事後分布の最新の推定値を出力する推定値出力ステップと、
を含んで実行することを特徴とするクラスタリング計算方法。
【請求項５】
前記事後分布推論手段は、
前記事後分布推定ステップにおいて、
前記ＥＭアルゴリズムのＥステップの処理として、ｄＨＤＰモデルにおいて予め定められた事前分布、観測分布および最大クラスタ数と、ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＭステップまでに推定された隠れ変数の事後分布の推定値とを読み込んで、推定対象の時刻に過去分を含めた演算対象時刻ごと、かつ、クラスタごと、かつ、前記演算対象時刻別の全データごとに、前記クラスタと、混合比と、前記クラスタ分布の時間変化の程度を表す重みとに関するパラメータの事後分布の値を演算する段階と、
前記ＥＭアルゴリズムのＭステップの処理として、前記ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＥステップまでに推定されたパラメータの事後分布の推定値とを読み込んで、推定対象の時刻に過去分を含めた演算対象時刻ごと、かつ、前記演算対象時刻別の全データごとに、２種類の隠れ変数の事後分布の値を推定し、前記２種類の隠れ変数のうち、前記クラスタ分布の時間変化の程度を表す重みに関連付けられた第１隠れ変数の事後分布の値についてはクラスタごとに演算し、前記第１隠れ変数および前記混合比に関連付けられた第２隠れ変数の事後分布の値については、前記演算対象時刻から過去に遡及した時刻ごとに演算する段階とを含み、
前記Ｅステップの処理と前記Ｍステップの処理とを交互に予め定められた回数だけ繰り返し実行することを特徴とする請求項４に記載のクラスタリング計算方法。
【請求項６】
前記事後分布推論手段は、
前記Ｅステップにおいて、観測量の集合データと、隠れ変数の事後分布の推定値とについては、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込み、
前記Ｍステップにおいて、観測量の集合データと、パラメータの事後分布の推定値とについては、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込むことを特徴とする請求項５に記載のクラスタリング計算方法。
【請求項７】
前記事後分布推論手段は、
前記Ｅステップにおいて、推定対象の時刻から事前に設定された設定値分だけ遡及した過去の時刻までを演算対象時刻として、前記クラスタと、前記混合比と、前記重みとに関するパラメータの事後分布の値を演算し、
前記Ｍステップにおいて、推定対象の時刻から事前に設定された設定値分まで遡及した過去の時刻までを演算対象時刻として、前記第１隠れ変数および第２隠れ変数の事後分布の値を演算することを特徴とする請求項５に記載のクラスタリング計算方法。
【請求項８】
前記事後分布推論手段は、
前記Ｅステップにおいて、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、事前に設定された再推定すべきクラスタの判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合にだけ、前記混合比に関するパラメータの事後分布の値を演算し、
前記Ｍステップにおいて、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、前記判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合にだけ、前記第１隠れ変数の事後分布の値を演算することを特徴とする請求項５に記載のクラスタリング計算方法。
【請求項９】
請求項１乃至請求項３のいずれか一項に記載のクラスタリング計算装置を構成する各手段としてコンピュータを機能させるためのクラスタリング計算プログラム。
【請求項１０】
請求項９に記載のクラスタリング計算プログラムが記録されたことを特徴とするコンピュータ読み取り可能な記録媒体。
【請求項１】
話者数が未知である会話の録音データから前記会話の話者数を推定するために、各話者を特徴付ける特徴量を抽出する特徴量抽出手段と、前記抽出された特徴量から前記各話者に対応する複数のクラスタを生成するときの複数の未知パラメータの値をそれぞれ推定するクラスタリング計算装置と、前記推定された複数のパラメータ値により前記会話の各話者を識別する識別手段とを有したダイアライゼーションシステムの前記クラスタリング計算装置であって、
前記抽出された特徴量を読み込む読込手段と、
前記読み込んだ特徴量をノンパラメトリックベイズモデルに対応したベクトルの観測量に変換する観測量生成手段と、
前記変換された観測量の集合データを蓄積記憶する観測量記憶手段と、
前記観測量の集合データから複数のクラスタをノンパラメトリックベイズモデルにより生成するときの複数のパラメータの事後分布の値をＥＭアルゴリズムによりそれぞれ推定および更新する事後分布推論手段と、
前記推定および更新された複数のパラメータの事後分布の値を蓄積記憶する推定値記憶手段と、
事前に設定された終了条件が成立したときに前記推定値記憶手段に記憶されている前記複数のパラメータの事後分布の最新の推定値を出力する出力制御手段と、
を備えることを特徴とするクラスタリング計算装置。
【請求項２】
前記事後分布推論手段は、
前記ＥＭアルゴリズムのＥステップの処理として、ｄＨＤＰ（dynamic Hierarchical Dirichlet Process）モデルにおいて予め定められた事前分布、観測分布および最大クラスタ数と、ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＭステップまでに推定された隠れ変数の事後分布の推定値とを読み込んで、過去分を含めた演算対象時刻ごと、かつ、クラスタごと、かつ、演算対象時刻別の全データごとに、前記クラスタと、混合比と、前記クラスタ分布の時間変化の程度を表す重みとに関するパラメータの事後分布の値を演算するＥステップ用計算手段と、
前記ＥＭアルゴリズムのＭステップの処理として、前記ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＥステップまでに推定されたパラメータの事後分布の推定値とを読み込んで、過去分を含めた演算対象時刻ごと、かつ、演算対象時刻別の全データごとに、２種類の隠れ変数の事後分布の値を推定し、前記２種類の隠れ変数のうち、前記クラスタ分布の時間変化の程度を表す重みに関連付けられた第１隠れ変数の事後分布の値についてはクラスタごとに演算し、前記第１隠れ変数および前記混合比に関連付けられた第２隠れ変数の事後分布の値については、演算対象時刻から過去に遡及した時刻ごとに演算するＭステップ用計算手段と、
前記Ｅステップの処理と前記Ｍステップの処理とを交互に予め定められた回数だけ繰り返し実行させる制御を行う収束判定手段と、
を備えることを特徴とする請求項１に記載のクラスタリング計算装置。
【請求項３】
前記Ｅステップ用計算手段は、
観測量の集合データと、隠れ変数の事後分布の推定値とについて、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込み、推定対象の時刻から事前に設定された設定値分だけ遡及した過去の時刻までを演算対象時刻として、前記クラスタと、前記混合比と、前記重みとに関するパラメータの事後分布の値を演算し、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、事前に設定された再推定すべきクラスタの判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合に、前記混合比に関するパラメータの事後分布の値を演算し、
前記Ｍステップ用計算手段は、
観測量の集合データと、パラメータの事後分布の推定値とについて、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込み、推定対象の時刻から事前に設定された設定値分まで遡及した過去の時刻までを演算対象時刻として、前記第１隠れ変数および第２隠れ変数の事後分布の値を演算し、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、前記判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合に、前記第１隠れ変数の事後分布の値を演算することを特徴とする請求項２に記載のクラスタリング計算装置。
【請求項４】
話者数が未知である会話の録音データから前記会話の話者数を推定するダイアライゼーションシステムにおいて、記憶手段と、読込手段と、観測量生成手段と、事後分布推論手段と、出力制御手段とを備えて、前記録音データから抽出された各話者を特徴付ける特徴量から前記各話者に対応する複数のクラスタを生成するときの複数の未知パラメータの値をそれぞれ推定するクラスタリング計算装置のクラスタリング計算方法であって、
前記読込手段によって、前記抽出された特徴量を読み込む特徴量読込ステップと、
前記観測量生成手段によって、前記読み込んだ特徴量をノンパラメトリックベイズモデルに対応したベクトルの観測量に変換し、前記変換された観測量の集合データを記憶手段に順次蓄積する観測量蓄積ステップと、
前記事後分布推論手段によって、前記観測量の集合データから複数のクラスタをノンパラメトリックベイズモデルにより生成するときの複数のパラメータの事後分布の値をＥＭアルゴリズムによりそれぞれ推定し、当該推定値を前記記憶手段に順次格納および更新する事後分布推定ステップと、
前記出力制御手段によって、事前に設定された終了条件が成立したときに前記記憶手段に記憶されている前記複数のパラメータの事後分布の最新の推定値を出力する推定値出力ステップと、
を含んで実行することを特徴とするクラスタリング計算方法。
【請求項５】
前記事後分布推論手段は、
前記事後分布推定ステップにおいて、
前記ＥＭアルゴリズムのＥステップの処理として、ｄＨＤＰモデルにおいて予め定められた事前分布、観測分布および最大クラスタ数と、ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＭステップまでに推定された隠れ変数の事後分布の推定値とを読み込んで、推定対象の時刻に過去分を含めた演算対象時刻ごと、かつ、クラスタごと、かつ、前記演算対象時刻別の全データごとに、前記クラスタと、混合比と、前記クラスタ分布の時間変化の程度を表す重みとに関するパラメータの事後分布の値を演算する段階と、
前記ＥＭアルゴリズムのＭステップの処理として、前記ハイパーパラメータの設定値と、過去から推定対象の時刻までに変換された観測量の集合データと、過去から最新のＥステップまでに推定されたパラメータの事後分布の推定値とを読み込んで、推定対象の時刻に過去分を含めた演算対象時刻ごと、かつ、前記演算対象時刻別の全データごとに、２種類の隠れ変数の事後分布の値を推定し、前記２種類の隠れ変数のうち、前記クラスタ分布の時間変化の程度を表す重みに関連付けられた第１隠れ変数の事後分布の値についてはクラスタごとに演算し、前記第１隠れ変数および前記混合比に関連付けられた第２隠れ変数の事後分布の値については、前記演算対象時刻から過去に遡及した時刻ごとに演算する段階とを含み、
前記Ｅステップの処理と前記Ｍステップの処理とを交互に予め定められた回数だけ繰り返し実行することを特徴とする請求項４に記載のクラスタリング計算方法。
【請求項６】
前記事後分布推論手段は、
前記Ｅステップにおいて、観測量の集合データと、隠れ変数の事後分布の推定値とについては、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込み、
前記Ｍステップにおいて、観測量の集合データと、パラメータの事後分布の推定値とについては、推定対象の時刻から事前に設定された定数または直前のＭステップの推定値を反映した変動数分だけ過去に遡及した時刻のデータを読み込むことを特徴とする請求項５に記載のクラスタリング計算方法。
【請求項７】
前記事後分布推論手段は、
前記Ｅステップにおいて、推定対象の時刻から事前に設定された設定値分だけ遡及した過去の時刻までを演算対象時刻として、前記クラスタと、前記混合比と、前記重みとに関するパラメータの事後分布の値を演算し、
前記Ｍステップにおいて、推定対象の時刻から事前に設定された設定値分まで遡及した過去の時刻までを演算対象時刻として、前記第１隠れ変数および第２隠れ変数の事後分布の値を演算することを特徴とする請求項５に記載のクラスタリング計算方法。
【請求項８】
前記事後分布推論手段は、
前記Ｅステップにおいて、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、事前に設定された再推定すべきクラスタの判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合にだけ、前記混合比に関するパラメータの事後分布の値を演算し、
前記Ｍステップにおいて、クラスタごとの演算のためにクラスタを更新する処理を実行するたびに、前記判断基準に基づいて、推定処理対象のクラスタを再推定すべきか否か判別し、再推定すべきクラスタである場合にだけ、前記第１隠れ変数の事後分布の値を演算することを特徴とする請求項５に記載のクラスタリング計算方法。
【請求項９】
請求項１乃至請求項３のいずれか一項に記載のクラスタリング計算装置を構成する各手段としてコンピュータを機能させるためのクラスタリング計算プログラム。
【請求項１０】
請求項９に記載のクラスタリング計算プログラムが記録されたことを特徴とするコンピュータ読み取り可能な記録媒体。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【公開番号】特開２０１０−１７５６１４（Ｐ２０１０−１７５６１４Ａ）
【公開日】平成２２年８月１２日（２０１０．８．１２）
【国際特許分類】
【出願番号】特願２００９−１５３３８（Ｐ２００９−１５３３８）
【出願日】平成２１年１月２７日（２００９．１．２７）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【Ｆターム（参考）】
【公開日】平成２２年８月１２日（２０１０．８．１２）
【国際特許分類】
【出願日】平成２１年１月２７日（２００９．１．２７）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【Ｆターム（参考）】
[ Back to top ]