以下、添付図面を参照しながら、本発明における隠れマルコフモデルの推定方法,推定装置および推定プログラムの好ましい各実施例について説明する。
先ず、本発明の具体的な実施例を提示する前に、当該実施例における独自のalpha-HMM再推定アルゴリズムを導入するまでの理論的な経緯について、以下の数式を参照しながら説明する。
一般に、隠れマルコフモデルによる学習アルゴリズムは、凸ダイバージェンスの最適化から始まって生成されている。まず、次の式に示すように、2つの確率密度p,qの間の凸ダイバージェンスについて考える。
ここで、YはK次元ユークリッド空間である。関数f(r)は、(0,∞)において凸である。関数f(r)の双対関数g(r)は、次式を満たす必要がある。
上記数1における「≧」という記号は、p=qを殆どどこでも保持する場合に限り保持する。f(1)の正規化が任意であるため、ここでは、Df(p・q)とDg(q・p)がpとqとの間の擬似距離として作用するように、f(1)を次のように選定する。
留意すべきなのは、次の2つの式が、カルバック・ライブラー・ダイバージェンス(Kullback-Leibler divergence)を発生させる、ということである。もしpとqが確率量の関数であれば、数1における積分は総和に置き換えられる。これは、以後説明する積分の全てにおいて同様である。
上記数1〜数5の関係は、対数より有能であるものの、対数と類似した関数であることを意味している。連続して2回微分可能な関数のクラスについて考えると、数2は以下の等式を与える。
また、r=1付近では以下の式がそれぞれ成り立つ。
ここで、o(1)は高次無限小の項を意味する。もし、乗算的に分離可能な凸関数を使用するのであれば、o(1)は不要であることに注目されたい。
以下の関数は、その一例である。
もし、上記の関数を、f(1)=0およびf”(1)=1であるようにシフトして正規化すれば、k=−c(1−c)となる。前記数7のコアとなる項は以下の式を有している。
上式において、次の式は興味深い関数であり、これはパラメータcによって凸性が調整できる単調関数である。
これは、以下の式の関係が保たれることを指し示すのに重要である。
したがって、L(c)(r)はパラメータcを有する法則化された対数と見なされる。つまり、c-対数と呼ぶことができる。しかしながら、パラメータcの増加はその数値を減少させる。そこで、ここでは以下の関係によって、法則化された対数をパラメータ化する。
このように、α-対数(alpha-log)は以下の式に示すような凸ダイバージェンスDf(p・q)から始まることによって得られる。
ここでの対数は次式のようになる。
注目すべきは、次の式に示す凸関数が、r=p/qとして上記数1に適用されるとαダイバージェンスが発生し、これはalpha-EMアルゴリズムにおいて重要な役割を果たす。
αダイバージェンスは、次の式であらわせる。
上式において、p=qが殆どどこでも適用される場合に限り、等式のゼロが達成される。
次に、alpha-HMM再推定アルゴリズムの原点となるalpha-EMアルゴリズムの期待値と最大化の処理ステップについて説明する。ここではPY|ψ(y|ψ)を、ψによってパラメータ化されたYにおける観測データyの確率密度または確率量とする。集合Yは、不完全データの集合と見なされる。x∈X(太字)は、未知または消失したデータを含む理想的観測の完全データあるいは増補データとする。そして、不完全データのpdf(probability density function:確率密度関数)あるいはpmf(probability mass function:確率量関数)は次式に示される。
ここでの積分範囲は、次式のようになり、pdfの積分はpmfの総和となる。
そして、条件付きのpdfまたはpmfは次式の通りとなる。
pdfまたはpmfの不完全なデータのalpha-log尤度比は、次式の通りである。
ここで、次の式におけるαダイバージェンスの計算は、alpha-EMアルゴリズムの基本的な関係を示している。
その式は、以下の通りである。
したがって、最大化のための目的関数は以下の式のような定量化(quantity)によって表される。
上記説明によって、alpha-EMアルゴリズムがalpha-HMM再推定アルゴリズムの開始点と位置付けられることがわかる。
alpha-EMアルゴリズムでは、次のステップで処理が行なわれる。
初期設定:最初のサイクル用のパラメータφを設定する。
E-ステップ:φが与えられたときのψの関数として、前記数25を計算する。
M-ステップ:次の式にしたがって、更新パラメータを計算する。
U-ステップ:φをψ*に置換して、収束するまではE-ステップに戻る。
もし、閉じた部分で正確な最大化が得られなければ、数26に示す最大化が改善した計算に置き換わる。
また、別なalpha-GEMアルゴリズムでは、次のステップで処理が行なわれる。
初期設定:最初のサイクル用のパラメータφを設定する。
E-ステップ:φが与えられたときのψの関数として、前記数25を計算する。
M-ステップ:次の式にしたがって、更新パラメータを計算する。
U-ステップ:φをψ+に置換して、収束するまではE-ステップに戻る。alpha-HMMの近似されたバージョンが、GEMアルゴリズムと見なすことができる点に注目されたい。
隠れマルコフモデル(HMM)の公式化の問題は、一連のランダムな観測変数データを発生させる最良のマルコフモデルを推定することにある。観測データの変数列は、次の式であらわせる。この式におけるtはマルコフモデルを生成するためのデータの順序を表すものであり、繰り返しのインデックスとは意味が異なる。
ここでの「最良」とは、隠れマルコフモデルの尤度が最大化されることを意味する。隠れマルコフモデルを公式化するために、次に定義する各確率を推定する必要がある。
(a)状態遷移確率
状態がiからjに遷移する確率のことであり、次式のように定義される。
上記式において、si∈Σは有限な範囲での総和Σの中の1つの状態を示すアルファベットの記号である。状態skから状態slへのつながりが無い場合、状態遷移確率akl=0となる。これは、状態遷移における従前のトポロジーを反映している。
(b)出力確率
状態sjにおいて出力ykが発生する確率のことであり、次式のように定義される。
ここで、上記式における各要素は以下の通りである。
(c)初期状態確率
最初の状態がiである確率のことであり、次式のように定義される。
推定される上記(a)〜(c)の各要素は、以下のように示される。
つまり、ここで用いる記号θは、状態遷移確率,出力確率および初期状態確率を含ませたものである。
一方、隠れマルコフモデル推定の問題に関し、前記数28で示したランダムな観測データの変数列は、その値が次の式であらわされる。ここで、tはデータの順序を表すインデックスである。
与えられた変数列は、その変数列を発生できる最大尤度モデルθを推定する。この隠れマルコフモデル推定の問題は、以下の解釈をもたらす。
(a)数28で示す観測データが、前述した不完全データYとして与えられる。
(b)次の式に示す状態遷移は、消失データZである。
(c)完全データは、次式のように示される。
この確率は、次式で示される。
(d)前述の三重項θは、数36に示す完全データの確率構造を提供する。
この解釈は期待値最大化によって、実際にはalpha-EMアルゴリズムによって、前記HMM問題の公式化を可能にする。そして、以下の重要な定理が得られる。
定理:ψとφを、2つのマルコフモデルの各パラメータセットとする。そして、α≦1に対して、次の数38の式は、数39という意味を持つ。
ここで数39の不等式は、数38の不等式が成立する場合に成立するものである。
証明:基本の方程式である数24において、次式に示す不完全データのalpha-log尤度比を設定する。
数35に示す状態S(太字)が消失データであるとすると、X=(S,Y)(各記号は何れも太字)であるので、前記数25のQ-関数は次式で算出できる。
前記数18の特性と数24の基本的な方程式によって、この定理が保たれる。
数41の特性は、数40の形態で推定されるマルコフモデルψが、マルコフモデルφよりも高い確率で数列y(太字)を生成することを意味する。すなわち、 数33に示すマルコフモデルクラスでのalpha-EMアルゴリズムの繰り返しは、必然的な局所最適性のコストで最良の隠れマルコフモデルを生成する。したがって、数33に示すマルコフモデルのパラメータで数38を繰り返すことが、HMMアルゴリズムの理論版となる。ここでは、α=1の場合(後述する高速化パラメータの値が、β=3である場合)は、収束限界であるために実用的ではない。
alpha-EMアルゴリズムとしての主な問題は、ソフトウェアとして実現できる具体的なコンピュータアルゴリズムをいかに提供するのかにある。数38のalpha-HMMは理論上の形態であり、これはlog-EMアルゴリズムと同様に一般的なHMMである。ソフトウェアで実行可能なalpha-HMMとして、ここでは2つのバージョンを提示するが、第1のバージョンは正確であるものの、未来の情報の計算を必要とする。すなわち、得られたアルゴリズムはnon-causalである。一方、本発明で特に提案する第2のバージョンは、時間シフトと確率の近似を用いることによって、未来の情報を含めることなく計算が可能となる。なお、ここで提示するalpha-HMMアルゴリズムの生成において、全てはpmfであるので、積分が総和となる。以下、それぞれのバージョンについて説明する。
(A)第1の正確なバージョン
・non-causalな更新
正確なバージョンは、条件付き微分によって得られる。このバージョンで得られる更新がnon-causalであることは、既に述べている。この形態は次のバージョンのcausalな形態を得るための原点となるものである。
次式に示すように、全ての更新式は、x=(s,y)(各記号は何れも太字)に対する前記Q-関数の最大化から得られる。
ここで、tはalpha-HMMに対する繰り返しの指標(インデックス)である。そして、パラメータ集合は次式に示すようになる。なお、明細書中、小文字のπ、a、b等の記号は集合として表記する場合は大文字で示すようにしている。
・Πθt+1(Πは太字、t+1はθの添え字)の推定
初期状態確率πi|θt+1(πは太字、t+1はθの添え字)の更新は、次の修正されたQ-関数における微分の最大化によって得られる。
この計算は、次式の微分を用いることによって達成できる。
そして、次式を用いることで、数44のλを除去する。
すると、次の更新式を与える。
上記数47は観測データy(太字)に対して計算できるように見えるが、左辺に含まれるθt+1が右辺に含まれている。これは、non-causalityの存在による自己撞着性を示している。しかしながら、これは次のバージョンの近似によって解決されることとなる。したがって、確率のための2つの更なる更新式を得ることができた。
・Aθt+1(Aは太字、t+1はθの添え字)の推定
状態遷移確率aij|θt+1(t+1はθの添え字)の更新は、同様に微分によって得られる。
ここで次式を用いることで、λを除去する。
すると、次の更新式を与える。
ここでのNij(s)(sは太字)は、状態sにおいてiからjに移行する状態遷移の数である。数50も右辺にθt+1を含んでいる。これも次のバージョンで解決されることとなる。
・Bθt+1(Bは太字、t+1はθの添え字)の推定
出力確率bjk|θt+1(t+1はθの添え字)の更新は、再度の微分によって得られる。
ここで次式を用いることで、λを除去する。
すると、次の更新式を与える。
ここでのNbjk(s)(sは太字、jkはbの添え字)は、数列s(太字)によって生じる出力確率bjkに対する事象の数である。数53も右辺にθt+1を含んでいる。これは次のバージョンで解決されることとなる。
(B)第2のCausalな近似バージョン
前述したように、第1のバージョンに示す一連の更新式(数47,数50,数53)は、自己撞着性を解決するために近似を必要とする。
・時間のシフト
alpha-EMアルゴリズムの収束は、適切な収束判定基準によって、次の式のようになることを意味する。
これは、以下の近似式を与える。ここで、o(1)は高位の無限小を意味する。
これは、繰り返し指標tをシフトすることによって、数47,数50,数53に示す更新式の自己撞着性を解決することになる。
・αのシフト
前記時間のシフトは、尤度比の期待値のための確率環境が反転されることを意味する。したがって、元のパラメータαはαcausalに変換される。αとαcausalの関係は、以下の誘導によって得られる。ここで、数47,数50,数53における尤度比の期待値の核となる部分を考慮すると、数47の場合は以下の関係が成り立つ。
この誘導は尤度比の時系列を維持しながら、期待値が現在の環境P(s|y,θt)(s,yは太字)から、未来の環境P(s|y,θt+1)(s,yは太字)に変化することを意味する。上記数56を見るとθt+1の初期確率πiを左辺として計算するのに、第一式の右辺では同じ時系列のθt+1に依存する関数(項)をしなければならないという問題点があり、このままでは計算できないと言う自己撞着性を含んでいた。そこで、まず第二式では第一式の分母分子を{}内の項のα乗でそれぞれ除算する。さらに第三式ではテーラー展開によって微小近似すると共に時間シフトの概念を導入しており、これによって自己撞着矛盾を解決している。言い換えると、右辺には左辺と同時系列のθt+1に依存する関数(項)を用いない形になっている。この詳細な考え方を以下に示す。
a)未来に関する尤度比を現在の環境で計算すると、次式のようになる。
b)未来に関する尤度比を未来の環境で計算すると、次式のようになる。
したがって、時間シフト前後の関係は以下の通りである。
そこから得られる関係式は、次式のようになる。
したがって、α=−1のlog-EMとlog-HMMは、ここでのαcausal=1の場合に相当する。なお、符号の使用を簡単化するために、以後、高速化パラメータの値となるαcausalを、次式のようにβとしてあらわす。すると、数60の右側の条件式はβ≦3を示すことになり、従来技術であるlog-EMとlog-HMMはβ=1の場合に対応することがわかる。(ただし、β=3だと後述の実験結果より発散してしまうことがわかっている。)
・テーラー展開
前記全ての更新式(数47,数50,数53)は、尤度比の計算に能力を必要とする。これは計算に時間がかかるものであり、Tが増加するに従って扱いにくくなる。テーラー展開は、alpha-logアルゴリズムにおける尤度比の長所を失うことなく、計算を簡単化することができるもので、尤度比は次式のように近似化される。
数55と数62を数47,数50,数53に適用すると、過去の情報を利用するcausalで、かつ計算上効率的なalpha-HMMアルゴリズムが得られる。注目すべきは、これが数27のalpha-GEMアルゴリズムに対応するということである。
(C)ソフトウェアで実行可能なalpha-HMMアルゴリズム
・Πθt+1(Πは太字、t+1はθの添え字)の推定
高次の項o(1)を放棄することによって数55と数62の近似式を結合し、数47に適用すると以下の更新式が与えられる。
・Aθt+1(Aは太字、t+1はθの添え字)の推定
数55と数62の近似式を適用すると、数50は以下の更新式を与える。
ここで、Naij|θt(ijはaの添え字、tはθの添え字)は、条件付き確率P(s|y,θt)(s,yは太字)の下でのNij(s)(sは太字)の期待値であり、不完全データy(yは太字)によって生じる状態遷移のカウント数により推定される。
・Bθt+1(Bは太字、t+1はθの添え字)の推定
数55と数62の近似式を適用すると、数53は以下の更新式を与える。
ここで、NbjK|θt(jkはbの添え字、tはθの添え字)は、条件付き確率P(s|y,θt)(s,yは太字)の下でのNbjK(s)(sは太字)の期待値であり、不完全データy(yは太字)によって生じる状態遷移のカウント数によって推定される。
・従来のlog-HMMとの比較
上記の数63,数64,数65は、ソフトウェアとして実行可能である。これらの更新式は、以下に示す共通の特性を有する。
(1)β=1すなわちα=−1の場合は、log-EMアルゴリズムから得られた従来のHMMアルゴリズムと一致する。
(2)数60と数61との間の関係により、αが行なうのと同様に、パラメータβは更新用の現在と過去の情報の各要素を調整する。
(3)βに依存する、すなわちαに依存するすべての過去の更新値は、格納データを参照することによってのみ得ることができる。
(4)現在の項は、log-HMMアルゴリズムに対する前方−後方アルゴリズムを用いることで効率的に計算できる。
上記特性(1)〜(4)は、alpha-HMMアルゴリズムの各更新サイクルがlog-HMMアルゴリズムのオーバーヘッドより僅かに増えてだけであることを示唆している。したがって、設計パラメータα、あるいは高速化パラメータβが適切に選択されれば、本発明のalpha-HMMアルゴリズムは、従来のlog-HMMアルゴリズムよりも高速に収束することが予測される。この違いが生じる理由は、数63、数64、数65の式において、従来のlog-EMから得られたHMMアルゴリズムに相当するβ=1を代入すると分母分子の第二項が消去され、分母分子の第一項が約分されて尤度比が残らなかったのに対し、β>1で設定すると分母分子の第二項が消去されず、単純に分母分子を約分できないために尤度比の影響を残すことができるようになったことだと考えられる。
以下、上記理論に基づく好ましい実施例について説明する。
実施例1、2では、観測データが離散的な場合において、観測データが単一配列(一本の観測データ)を有する場合、複数配列(M本の観測データ)を有する場合について、それぞれ説明する。
実施例3、4では、観測データが連続的な場合において、観測データが単一配列(一本の観測データ)を有する場合、複数配列(M本の観測データ)を有する場合について、それぞれ説明する。
実施例5、6では、観測データが半連続的な場合において、観測データが単一配列(一本の観測データ)を有する場合、複数配列(M本の観測データ)を有する場合について、それぞれ説明する。
実施例7では、観測データが離散的かつ連続的な場合、すなわち離散系列と連続系列が混在している場合において、単一配列を有する場合と複数配列を有する場合をまとめて説明する。
図1は、本発明の第1実施例において、上述したalpha-HMM再推定アルゴリズムを実行可能にするプログラムの処理手順をあらわしたものである。
同図において、1は合計でT個の観測データを格納する記憶手段としてのレジスタで、各観測データは時刻τが1から順に時系列に並んで格納される。ここでyτは個々のデータ値を示し、yはT個のデータ値の集合を示しており、本実施例では配列が1本の観測データをレジスタ1に格納している。推定装置10は、レジスタ1に記憶される一列のデータ値が、どのような確率構造(モデル)を有しているのかを、HMMの未知パラメータを算出する以下のステップS1〜ステップS10の各手順に従って推定解析するものである。
ステップS1は、前述した高速化パラメータの値βを設定する部分である。従来のHMM推定アルゴリズムは、β=1の場合に相当する。つまり、β=1という特殊な条件下であれば、従来のalpha-HMM再推定アルゴリズムであっても、上記確率モデルの計算は可能である。しかし、本実施例で提案する新たなalpha-HMM再推定アルゴリズムでは、確率モデルの計算が可能な高速化パラメータの値βを1≦β<3の範囲に拡張することができる。これは、従来のβ=1を特例として含むものである。なお、βが大きな値であるほど、推定装置10としての処理の高速性は増すが、収束性を保持するにはβ<3でなければならない。
ステップS2は、初期確率と収束判定値を決める部分である。これは、後述するステップS3〜ステップS8の手順を繰り返す前に行なわれる。推定装置10が最終的に算出しようとするHMMの確率構造は、次式のようにあらわせる。これは、前述の数37と等しい。
上記式において、πS0(0はsの添え字)は最初の状態s0における確率であり、aSτ-1Sτ(τ−1およびτはsの添え字)は時刻τ−1の状態sτ-1が時刻τの状態sτに移る確率であり、bSτ(yτ)(τはsの添え字)は時刻τの状態sτに移行したときにyτが出力される確率である。推定装置10は、レジスタ1に格納された観測データを読み出して、その観測データが最も出現しやすくなる初期状態確率πと、状態遷移確率aと、出力確率bとを推定するが、このステップS2では、次式に示すように、それらの確率π,a,bの初期値を決定する。
上式において、θ0は数43にあるように、0回目の繰り返しにおける初期の確率π,a,bの組合せを示し、推定装置10は、その条件で状態iから状態jに移る状態遷移確率aij|θ0(0はθの添え字)の値と、状態jで状態kが出力される出力確率bjk|θ0(0はθの添え字)の値と、最初に状態iとなる初期状態確率πi|θ0(0はθの添え字)の値とをそれぞれ決定する。
またステップS2では、対数尤度に基づく収束範囲を決めることで、ステップS3〜ステップS9の繰り返しを終了させるための収束判定値を決定する。この収束判定値はステップS8で用いられ、具体的には後述する数83で示される。
これらの値をステップS1,S2で決定すると、ステップS3の手順に移行して、確率量とカウント値の設定が行なわれる。ステップS3は、最初にステップS2で求めた初期確率を利用し、それ以降はステップS9で更新された確率量とカウント値を利用して、数68に示す各確率値と、数69に示す各カウント値を、実際に推定装置10のメモリ(図示せず)に設定する部分である。
次のステップS4は、レジスタ1から観測データを読み出して、ステップS3で設定された上式の確率量から、数70と数71に示す前向き確率を計算する部分である。なお、ここに示すαは確率値であり、前述したパラメータの値αとは異なる。ここでは、数70に示す確率値αが計算され、そこから数71に示す確率値(尤度)P(y|θt)が計算される。
次のステップS5は、レジスタ1から観測データを読み出して、ステップS3で設定された上式の確率量から、数72に示す前向き確率を計算する部分である。なお、ここに示すβは確率値であり、前述したパラメータの値βとは異なる。ここでは、数72に示す確率値βが計算される。
上記前向き確率と後向き確率は、計算回数を減らすために、既存のHMM推定アルゴリズムにも組み込まれていたものである。
続くステップS6は、ステップS4で計算された確率値αと、ステップS5で計算された確率値βを用い、レジスタ1から観測データを読み出して、数73に示す状態遷移のカウント値と、数74に示す出力のカウント値をそれぞれ計算する部分である。
次のステップS7は、t+1回目に繰り返される確率π,a,bの組合せを条件として、状態iとなる初期状態確率πi|θt+1(t+1はθの添え字)の値と、状態iから状態jに移る状態遷移確率aij|θt+1(t+1はθの添え字)の値と、状態jで状態kが出力される出力確率bjk|θt+1(t+1はθの添え字)の値とを、ステップS4で求めたt回目の確率P(y|θt)と、それよりも1回前の確率P(y|θt-1)を利用し、且つステップS5で求めたt回目のカウント値Naij|θt,Nbjk|θt(ijはaの添え字、jkはbの添え字、tはθの添え字)と、それよりも1回前のカウント値Naij|θt-1,Nbjk|θt-1(ijはaの添え字、jkはbの添え字、t-1はθの添え字)を利用して計算する部分であり、これは前記数63,数64,65に対応している。
このステップS7では、前記ステップS1で設定したパラメータの値βが用いられ、本アルゴリズムの主要な特徴部分となっている。特に、添え字t+1の部分の確率を計算する上で、時間シフトした添え字tや添え字t-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み(高速化パラメータの値)が1≦β<3に拡張したことが、その効果として挙げられる。
本実施例では、配列が1本の観測データをレジスタ1に格納しているが、その場合のステップS7における確率用の更新式は、次式のようにあらわせる。
ステップS8は、計算の収束を判定する部分である。ここでは、直前のステップS7で計算された新たな確率量に基づく尤度Pを用い、その尤度Pと前記ステップS2における収束判定値との比較により、ステップS3〜ステップS9の計算が収束したか否かを判定する。近似的には、繰り返しの回数を指定して、指定回数に達したら、計算が収束したと判定してもよい。
ステップS8において計算が収束していないと判定された場合、ステップS9に移行して手順が繰り返される。ステップS9は、計算された確率量と事象のカウント値を、繰り返し計算のために更新する部分である。本実施例では、次式に示すように、2つの過去値をシフトさせる。
ここでは、繰り返しの回数が1つ増えることにより、t回目の確率量およびカウント値がt-1回目の確率量およびカウント値に更新され、t+1回目の確率量およびカウント値がt回目の確率量およびカウント値に更新される。なお、ステップS3とステップS9を一纏めにして、更新された確率量およびカウント値をそのまま次の繰り返しのために、推定装置10のメモリに設定記憶させてもよい。
一方、ステップS8において計算が収束していると判定された場合、ステップS10に移行して、次式に示す計算された三組の確率量を用いたHMMを採用する。
上記数79で計算した各値を用いて、推定装置10は数66に示すHMMの確率構造を出力することができる。
図2は、図1に示すアルゴリズムを実現する推定装置10の構成を示している。同図において、推定装置10は前記レジスタ1を内蔵する入力手段12と、レジスタ1に格納した観測データを適宜読み出し、上記ステップS1〜ステップS10の手順を実行して数66に示すHMMの確率構造を推定する推定手段14と、この推定手段14で得たHMMの確率構造を出力する出力手段16とにより概ね構成される。本実施例の入力手段12は、1本の配列を有する観測データをレジスタ1に格納できる構造となっているが、複数本の配列を有する観測データを取り扱う推定装置100については、次の実施例で詳しく説明する。
推定装置10のハードウェア構成は、例えば演算処理部としてCPUを備えたコンピュータで実現することができる。その場合、図1に示す推定装置10の処理手順を実行するプログラムが、メモリなどの記録媒体に記憶される。当該プログラムをどこに記憶するのかは限定せず、例えば通信手段を介してプログラムがコンピュータにダウンロードされる構成であってもよい。
推定手段14は、前記ステップS1における高速化パラメータの値βや、ステップS2における初期の確率値および収束判定の条件を設定する初期設定手段22と、ステップS3における確率量およびカウント値の設定や、ステップS9における確率量およびカウント値の更新を行なう更新設定手段24と、ステップS4における前向き確率の計算や、ステップS5における後ろ向き確率の計算を行なうと共に、そこからステップS6における状態遷移のカウント値と出力のカウント値をそれぞれ計算し、さらにステップS7における新たな確率量の計算を行なう演算手段26と、ステップS8における計算の収束判定を行ない、計算が収束していなければ、前記更新設定手段24による確率量およびカウント値の更新を行なわせる一方で、計算が収束していれば、ステップS10において、演算手段26が直前に計算した新たな確率量を、HMMの最終的なパラメータ値として確定させる判定手段28とを備えている。
また推定装置10には、前記高速化パラメータの値や、初期の確率値および収束判定の条件や、更新設定される確率量およびカウント値の値を読み書き可能に格納するメモリ30の他に、必要に応じて高速化パラメータの値や収束判定の条件を操作入力するキーボードやマウスなどの操作手段32が、推定手段14に接続して設けられる。これらの各装置構成により、上述したステップS1〜ステップS10の手順が実行される。
次に、本実施例における試験用の観測データでの実験結果を説明する。
(A)生成されたデータを用いた速度の評価
観測データの集合が強い凸性でなければ(そして強い凹面でもなければ)、全ての最適化アルゴリズムは局所的に最適となる。存在するデータの殆ど全ては、そのような性質を持っている。従来のあるいはlog-HMMアルゴリズムは例外的ではないし、alpha-HMMアルゴリズムも例外的ではない。したがって、そうしたHMMアルゴリズムを適用する前に、発生のメカニズムが分かっている人工的なデータ集合を生成することが必要になってくる。
ここでは、数33に示す数列をレジスタ1に格納する入力データとし、そのパラメータを次式のように指定する。
このマルコフ連鎖を使用し、[0, 1]上の一様な乱数によって、1000個のサンプルを有するデータ集合を生成した。完全な推定には、以下の困難があることを先に述べる。
(a)各データ集合の殆どが、局所的最適を与える対象である。
(b)たとえ局所的最適が避けられたとしても、正確な推定には無限の長さの数列が必要である。
(c)コンピュータが生成した乱数は単なる擬似的な乱数にすぎない。
したがって、以下が満たされれば満足するものとする。
(1)数33に示すパラメータ集合の数値が、数80からずれないこと
(2)alpha-HMMアルゴリズムの収束したパラメータ集合がlog-HMMのそれと近いこと
ステップS2における初期値は、次のように設定する。これにより、初期状態を除いて事前の情報は必要としないものとする。
以下の理由のため初期状態は固定した。
(a)分布された初期確率は、不安定な局所最適性を得やすい。これは、alpha-HMMアルゴリズムの各収束速度を比較するのに不要な曖昧さをもたらす。
(b)HMM推定の後の認識問題において、一つの初期状態はビタビ(Viterbi)アルゴリズムによって選択される。
図3は、対数尤度(log-likelihood)、すなわち繰り返し回数を考慮した尤度P(y|θt)に関する収束の傾向を示している。ここでは、数60および数61を満たす多様な値で実験した。
この図から以下のことがわかる。
(1)β=αcausalが増加するにつれて対数尤度の立ち上がりが早くなる。すなわち、収束するまでの繰り返し回数が少なくなる。
(2)alpha-EMアルゴリズムの限界(α=1)は、alpha-HMMアルゴリズムではβ=αcausal=3に相当するが、この場合、発散してしまう。
さらに2つの視点で、alpha-HMMアルゴリズムの性能を比較する必要がある。一つは局所最適の比較であり、もう一つはCPU時間である。局所最適の比較では、次式に示すように、θに対する全パラメータと共に、それに付随した対数尤度の値を一覧にする必要がある。
上記の実験結果の一覧を見るとβを1〜2.75の範囲で変化させても対数尤度(LL)の値がほとんど変化せずに正確さを保っていることがわかる。この実験ではβ≧3に設定した実験は発散する結果になるために省略しているが、3に近い2.75に設定しても実験結果にそれほど影響なく高速化できることがわかる。ただし、このような実験結果の収束限界となるβの値は入力データによって若干異なる。
次は、収束するまでのCPU時間の比較である。この比較はソフトウェアの実行手段とハードウェアの手法に依存する。ここでは、alpha-HMMアルゴリズムを標準的なPOSIX環境におけるC‐コードとして実行する。これらのコードは標準的なシングルコアプロセッサで走らせている。前記ステップS8で行われる収束は、以下の判定基準で測定した。
上式において、Pnewは今回の計算で得た確率Pの値であり、Poldは前回の計算で得た確率Pの値である。この収束判定基準は更新前後で確率の値がほとんど変化していなければ収束したと判断するものである。この収束判定式結果は次の表に示す通りであり、ここでは繰り返し回数とCPU時間による収束の比較を示している。
局所最適性とCPU時間とを合せて比較すると、以下のことがわかる。
(a)収束する各対数尤度は非常に接近しているものの、何れも局所最大値が存在する。
(b)対数のまたは従来のHMMアルゴリズムは、常に最良の対数尤度をもたらす訳ではない。log-HMMアルゴリズムよりも良好な局所最大値が、alpha-HMMアルゴリズムで得られた。
(c)改良したalpha-HMMアルゴリズムによって、繰り返し回数とCPU時間の高速化が成し遂げられたのは明らかである。
(c)繰り返し回数の改善が、直接CPU時間に影響している。これは前記数63,数64および数65に示すように、予め計算されてメモリに格納される過去の情報を利用しているからであり、従来のlog-HMMアルゴリズムに比べて、ソフトウェアの複雑さを増やしてはいない。
(d)一方、αcausal=β=3すなわちα=1の場合は、発散する。これは、alpha-EMアルゴリズムの許容範囲と一致する。したがって、数62による近似は許容される。αcausal=β=1すなわちα=−1の場合において、数63,数64および数65はlog-HMMアルゴリズムの正確な形態であることが重要である。この意味において、数63,数64および数65の更新式はlog-HMMアルゴリズムの拡張版ということができる。
前記図3を参照すると、従来のalpha-HMMアルゴリズムで唯一ステップS6からステップS7への計算が可能であったβ=1の結果に対して、今回提案するalpha-HMMアルゴリズムで計算が可能になったβ=2.75の結果では、少ない繰り返し回数で推定装置10の処理が収束していることが判る。具体的には、上記表1において、β=1の場合は繰り返し回数が263回で、CPUの使用時間が0.253秒となっているが、β=2.75の場合は繰り返し回数が70回、またCPUの使用時間が0.068秒に減少する。表1の最右列にあるのは、β=1のときの繰り返し回数とCPUの使用時間を1としたときの速度比を計算したもので、例えばβ=2.75に設定すれば繰り返し回数は3.76倍、またCPUの使用時間は3.72倍に高速化する。
以上の実験結果から、本実施例では従来よりも高速なHMMアルゴリズムが示された。本推定方法は、基になっているalpha -EMアルゴリズムを反映しており、改良したalpha-HMMアルゴリズムと呼ぶことができる。改良したalpha-HMMアルゴリズムは、従来のalpha-HMMアルゴリズムまたはlog-HMMアルゴリズムより優れている。計算上の複雑さの増加が非常に少ないので、繰り返し回数の抑制が直接的に期待したCPUの高速化を実現した。
また、改良したalpha-HMMアルゴリズムで採用する数63,数64および数65の更新式は、メモリに格納した過去の情報を追加的に必要とするだけである。したがって、既存のコードからのソフトウェアのバージョンアップは難しくない。この場合、alpha-対数の曲率を制御するαcausal=α+2が設計パラメータとなる。
さらに、より少ないサンプルで未来の情報を利用する方法について、確認を行なった。この方法は、テーラー展開と共にnon-causalな形態を使用したことにある。更新式の一例は、次式のように表せる。
ここで、θ^t+1とθ^t(^はθの上部に記される)は、より少ないサンプル(例えば1000個のサンプルから200個を取り出す)から統計量が推定されることを示している。すなわち、 サンプルの一部が未来を考えるのに使われる。この方法は数63,数64および数65の更新式よりも大きなオーバーヘッドを必要とするが、高速化の利益の方がこの余計な重荷よりも勝っている。もし、入力源がよく混合されている、あるいはエルゴードであると演繹的知識をユーザーが持っていれば、この方法を利用できる。実験では、高速化が数63,数64および数65の方法よりも悪くなかった。
以上のように本実施例では、入力される観測データを時系列に格納する記憶手段としてのレジスタ1と、観測データがどのような確率モデルであるのかを、HMMの未知パラメータを算出することで推定する推定手段14とを備えたHMMの推定装置10において、推定手段14は、HMMの高速化パラメータの値βを設定する初期設定手段22と、前記HMMの未知パラメータとして、状態遷移確率a,出力確率b,初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび出力の期待値Nbを更新設定する更新設定手段24と、更新設定手段24で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる(数68,数69を参照)と共に、前記レジスタ1から読み出した観測データと、初期設定手段22で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する(数63,数64,数65,数75,数76,数77を参照)演算手段26と、演算手段26による計算の収束を判定し、計算が収束していなければ、演算手段26で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段26で計算した新たな各確率量を最終的な値として出力させる判定手段28とを備えている。
このようにすれば、HMMの未知パラメータとして、状態遷移確率a,出力確率b,初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび出力の期待値Nbを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なHMM推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
また、ここでの演算手段26は、t+1回目の状態iとなる初期状態確率πi|θt+1(t+1はθの添え字)の値を、前記数63の式で計算し、t+1回目の状態iから状態jに移る状態遷移確率aij|θt+1(t+1はθの添え字)の値を、t回目およびt-1回目における状態遷移の期待値Naij|θtおよびNaij|θt-1(ijはaの添え字、tまたはt-1はθの添え字)を利用して、前記数64の式で計算し、状態jで状態kが出力される出力確率bjk|θt+1(t+1はθの添え字)の値を、t回目およびt-1回目における前記出力の期待値Nbjk|θtおよびNbjk|θt-1(jkはbの添え字、tまたはt-1はθの添え字)を利用して、前記数65の式で計算する構成となっている。
そのため、特に高速化パラメータの値βが1でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
さらに、ここでの更新設定手段24は、高速化パラメータの値を1<β<3に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが3未満の範囲で、βを1よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
なお上述した本実施例の作用効果は、初期設定手段22としての動作を実行する初期設定ステップと、更新設定手段24としての動作を実行する更新設定ステップと、演算手段26としての動作を実行する演算ステップと、判定手段28としての動作を実行する判定ステップとを備えたHMMの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段14として、コンピュータに機能させるHMMの推定プログラムであっても、同様に発揮される。
図4は、本発明の第2実施例におけるプログラムの処理手順をあらわしたものである。
同図において、101は合計でT個の観測データを格納する記憶装置としてのレジスタで、各観測データのそれぞれは、時刻τが1から順に時系列に並んで格納される。ここでyτは個々のデータ値を示し、yはT個のデータ値の集合を示し、nは何本目の配列かをあらわすインデックスを示しており、本実施例では複数の配列を有するM(=2以上の整数)本の観測データをレジスタ101に格納している。推定装置100は、レジスタ101に記憶されるM列のデータ値が、どのような確率構造(モデル)を有しているのかを、以下のステップS11〜ステップS20の各手順に従って推定解析するものである。
ステップS11は、前述した高速化パラメータの値βを設定する部分である。これは上記第1実施例のステップS1に相当するもので、高速化パラメータの値βは、従来のβ=1を特例として含む1≦β<3の範囲に拡張される。なお、βが大きな値であるほど、推定装置100としての処理の高速性は増すが、収束性を保持するにはβ<3でなければならない。
ステップS12は、初期確率と収束判定値を決める部分である。これは、前記ステップS2に相当するもので、後述するステップS13〜ステップS18の手順を繰り返す前に行なわれる。推定装置100は、レジスタ101に格納された観測データを読み出して、その観測データが最も出現しやすくなる初期状態確率πと、状態遷移確率aと、出力確率bとを推定するが、このステップS12では、次式に示すように、それらの確率π,a,bの初期値を決定する。
上式において、θ0は数43にあるように、0回目の繰り返しにおける初期の確率π,a,bの組合せを示し、推定装置100は、その条件で状態iから状態jに移る状態遷移確率aij|θ0(0はθの添え字)の値と、状態jで状態kが出力される出力確率bjk|θ0(0はθの添え字)の値と、最初に状態iとなる初期状態確率πi|θ0(0はθの添え字)の値とをそれぞれ決定する。
またステップS12では、対数尤度に基づく収束範囲を決めることで、ステップS13〜ステップS19の繰り返しを終了させるための収束判定値を決定する。この収束判定値はステップS18で用いられ、具体的には前記数83で示される。
これらの値をステップS11,S12で決定すると、ステップS13の手順に移行して、確率量とカウント値の設定が行なわれる。ステップS13は、最初にステップS12で求めた初期確率を利用し、それ以降はステップS19で更新された確率量とカウント値を利用して、数86に示す各確率値と、数87に示す各カウント値を、実際に推定装置100のメモリ(図示せず)に設定する部分である。
次のステップS14は、レジスタ101から観測データを読み出して、ステップS3で設定された上式の確率量から、数88と数89に示す前向き確率を計算する部分である。なお、ここに示すαは確率値であり、前述したパラメータの値αとは異なる。ここでは、数88に示す確率値αが計算され、そこから数89に示す確率値(尤度)P(y(n)|θt)が計算される。
次のステップS15は、レジスタ101から観測データを読み出して、ステップS13で設定された上式の確率量から、数90に示す前向き確率を計算する部分である。なお、ここに示すβは確率値であり、前述したパラメータの値βとは異なる。ここでは、数90に示す確率値βが計算される。
上記前向き確率と後向き確率は、計算回数を減らすために、既存のHMM推定アルゴリズムにも組み込まれていたものである。
続くステップS16は、ステップS14で計算された確率値αと、ステップS5で計算された確率値βを用い、レジスタ101から観測データを読み出して、数91に示す状態遷移のカウント値と、数92に示す出力のカウント値をそれぞれ計算する部分である。
次のステップS17は、t+1回目に繰り返される確率π,a,bの組合せを条件として、状態iとなる初期状態確率πi|θt+1(t+1はθの添え字)の値と、状態iから状態jに移る状態遷移確率aij|θt+1(t+1はθの添え字)の値と、状態jで状態kが出力される出力確率bjk|θt+1(t+1はθの添え字)の値とを、ステップS14で求めたt回目の確率P(y(n)|θt)と、それよりも1回前の確率P(y(n)|θt-1)を利用し、且つステップS5で求めたt回目のカウント値N(n) ajk|θt,N(n) bjk|θt(jkはaの添え字、tはθの添え字)と、それよりも1回前のカウント値N(n) ajk|θt-1,N(n) bjk|θt-1(jkはaの添え字、t-1はθの添え字)を利用して計算する部分であり、これは前記数63,数64,65に対応している。
このステップS17では、前記ステップS11で設定したパラメータの値βが用いられ、本アルゴリズムの主要な特徴部分となっている。特に、添え字t+1の部分の確率を計算する上で、時間シフトした添え字tや添え字t-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み(高速化パラメータの値)が1≦β<3に拡張したことが、その効果として挙げられる。
本実施例では、配列がM本の観測データをレジスタ101に格納しているが、その場合のステップS17における確率用の更新式は、次式のようにあらわせる。
ステップS18は、計算の収束を判定する部分である。ここでは、直前のステップS17で計算された新たな確率量に基づく尤度Pを用い、その尤度Pと前記ステップS12における収束判定値との比較により、ステップS13〜ステップS19の計算が収束したか否かを判定する。近似的には、繰り返しの回数を指定して、指定回数に達したら、計算が収束したと判定してもよい。
ステップS18において計算が収束していないと判定された場合、ステップS19に移行して手順が繰り返される。ステップS19は、計算された確率量と事象のカウント値を、繰り返し計算のために更新する部分である。本実施例では、次式に示すように、2つの過去値をシフトさせる。
ここでは、繰り返しの回数が1つ増えることにより、t回目の確率量およびカウント値がt-1回目の確率量およびカウント値に更新され、t+1回目の確率量およびカウント値がt回目の確率量およびカウント値に更新される。
一方、ステップS18において計算が収束していると判定された場合、ステップS20に移行して、次式に示す計算された三組の確率量を用いたHMMを採用する。
上記数97で計算した各値を用いて、推定装置10は数66に示すHMMの確率構造を出力することができる。
図5は、図4に示すアルゴリズムを実現する推定装置100の構成を示している。同図において、推定装置100は前記レジスタ101を内蔵する入力手段112と、レジスタ101に格納した観測データを適宜読み出し、上記ステップS11〜ステップS20の手順を実行して数66に示すHMMの確率構造を推定する推定手段114と、この推定手段114で得たHMMの確率構造を出力する出力手段116とにより概ね構成される。本実施例の入力手段112は、M本の配列を有する観測データをレジスタ101に格納できる構造となっている。推定装置100のハードウェア構成については、前記第1実施例の推定装置10と同様であるため、ここでは説明を省略する。
推定手段114は、初期設定手段122と、更新設定手段124と、演算手段126と、判定手段128とを備えている。これらは前記第1実施例の初期設定手段22,更新設定手段24,演算手段26および判定手段28にそれぞれ対応するもので、取り扱う観測データがM本の配列になった以外は、第1実施例と同様に機能する。
また推定装置110には、前記高速化パラメータの値や、初期の確率値および収束判定の条件や、更新設定される確率量およびカウント値の値を読み書き可能に格納するメモリ130の他に、必要に応じて高速化パラメータの値や収束判定の条件を操作入力するキーボードやマウスなどの操作手段132が、推定手段114に接続して設けられる。これらの各装置構成により、上述したステップS11〜ステップS20の手順が実行される。
以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ101と、観測データがどのような確率モデルであるのかを、HMMの未知パラメータを算出することで推定する推定手段114とを備えたHMMの推定装置100において、推定手段114は、HMMの高速化パラメータの値βを設定する初期設定手段122と、前記HMMの未知パラメータとして、状態遷移確率a,出力確率b,初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび出力の期待値Nbを更新設定する更新設定手段124と、更新設定手段124で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる(数68,数69を参照)と共に、前記レジスタ101から読み出した観測データと、初期設定手段122で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する(数63,数64,数65,数93,数94,数95を参照)演算手段126と、演算手段126による計算の収束を判定し、計算が収束していなければ、演算手段126で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段126で計算した新たな各確率量を最終的な値として出力させる判定手段128とを備えている。
このようにすれば、HMMの未知パラメータとして、状態遷移確率a,出力確率b,初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび出力の期待値Nbを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なHMM推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
また、ここでの演算手段126は、t+1回目の状態iとなる初期状態確率πi|θt+1(t+1はθの添え字)の値を、前記数63の式で計算し、t+1回目の状態iから状態jに移る状態遷移確率aij|θt+1(t+1はθの添え字)の値を、t回目およびt-1回目における状態遷移の期待値Naij|θtおよびNaij|θt-1(ijはaの添え字、tまたはt-1はθの添え字)を利用して、前記数64の式で計算し、状態jで状態kが出力される出力確率bjk|θt+1(t+1はθの添え字)の値を、t回目およびt-1回目における前記出力の期待値Nbjk|θtおよびNbjk|θt-1(jkはbの添え字、tまたはt-1はθの添え字)を利用して、前記数65の式で計算する構成となっている。
そのため、特に高速化パラメータの値βが1でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
さらに、ここでの更新設定手段124は、高速化パラメータの値を1<β<3に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが3未満の範囲で、βを1よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
なお上述した本実施例の作用効果は、初期設定手段122としての動作を実行する初期設定ステップと、更新設定手段124としての動作を実行する更新設定ステップと、演算手段126としての動作を実行する演算ステップと、判定手段128としての動作を実行する判定ステップとを備えたHMMの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段114として、コンピュータに機能させるHMMの推定プログラムであっても、同様に発揮される。
上記の実施例1、2では、観測データが離散的な場合について述べた。本実施例では、観測データが単一配列で連続的な場合の実施形態について説明する。
この観測データが連続的な場合とは、観測データのデータ列の各々の値が波のように変動しているような場合である。
具体的には、図6の観測データが離散的な場合のモデル図における出力確率bik,bjkが、図7の観測データが連続的な場合のモデル図のように確率密度関数である分岐確率cik,cjkに置き換わるというものである。
(A)前記第1実施例と前記第2実施例の式の変形
ここで、本実施例の具体的な内容に入る前に、離散的な場合について上記実施例と同様な説明を重複して行なう。
これは、実施例1、2と実施例3〜8との間において、数式の簡略化や符号の変更等が入っており、その点を踏まえて整合性の取れた説明をするという便宜上の理由のためである。なお、上記実施例1、2と実質的な差異はない。
まず、ソースデータとマルコフモデルの各種パラメータについて説明すると、HMMによるモデル化において、次式のようなソースデータ列が与えられる。ここで、tはデータの順序を表すインデックスである。
各ytは、スカラーあるいはベクトルである。上式は、単一配列の場合であるが、M本の複数配列で与えられるならば、ソースデータ列は次式の通りである。記号{}は集合を示している。
まず、単一配列では、数98の観測データy(yは太字)が与えられる。この場合のHMMの課題は、次式の最尤推定法(Maximum Likelihood Estimation:MLE)の認識において、最良なモデルを見つけることである。
ここで、S(Sは太字)は、次式のような状態遷移系列を意味している。
ランダムな変数とそれらの値のために、小文字のy(yは太字)およびs(sは太字)は同様に用いられる(通常、ランダム変数は大文字で示される)。
各確率は以下の数102〜数105の通りである。
・初期状態確率
・状態遷移確率
ここで、もし、状態siから状態sjへのつながりが無ければ、状態iから状態jに移る状態遷移確率はaij=0となる。これは、状態遷移のための前のトポロジーを反映している。
・出力確率
ここで、上式の各要素は次式の通りである。
上記数102〜数105に示した確率の集合をまとめてθとして次式に示す。
そして、確率的なデータ構造について以下のように解釈している。不完全なデータはy(yは太字)であり、推定される消失データはs(sは太字)であり、完全なデータは、数100の確率量と数106のパラメータを有しており、次式に示されるx(xは太字)として定義する。
・Alpha-EMアルゴリズム
不完全なデータ、消失データ、完全なデータによるHMMの解釈は、EMアルゴリズムに匹敵する。本紙の目的は、新たなHMM推定アルゴリズムを見つけることであるため、alpha-EMアルゴリズムから始まる道筋を示している。
観測データy(yは太字)は推定すべきパラメータの全体を表すψによってパラメータ化されるため、Py|ψ(y|ψ)(yは太字)を確率密度あるいは確率量とする。x∈X(xとXは太字)を、消失データを含む理想的な観測結果である完全なデータ、あるいは拡張されたデータとする。そして、不完全なデータ確率密度関数(pdf)、あるいは確率量関数(pmf)は次式で示される。
ここで、積分する範囲を次式とする。
pdfの積分は、pmfの総和となる。そして、条件付きのpdfあるいはpmfは次式である。
alpha-EMアルゴリズムでは、次式のようなalpha対数が用いられる。
ここで、α=−1の場合は対数、すなわち、次式である。
alpha-EMアルゴリズムは、alpha対数に関して、不完全なデータの尤度比を考慮する必要があるため、次式のようになる。
ここで、φとψは、繰返し最大化ステップにおける数106に対して、古いモデルのパラメータと新しいモデルのパラメータを示している。それから、alpha-EMアルゴリズムの基本的な方程式が次式として得られる。
ここで、D(α)は、常に負とならない2つの条件付きの確率Px|y,φ(x|y,φ)(xとyは太字)とPx|y,ψ(x|y,ψ)(xとyは太字)の間のalphaダイバージェンスである。上式において、次式のQ関数が重要である。Eは定量化(quantity)によって表される最大化のための目的関数を示している。
数114のために、もし、このQ関数が正ならば、a<1の範囲で、数114の左辺である不完全なデータのalpha対数尤度比も正となる。したがって、alpha-EMアルゴリズムとその変化形であるalpha-GEMアルゴリズムについて以下に示す。
[Alpha-EMアルゴリズム]
初期化: 数106の初期値を選択し、φとして用いる。
E-Step: 数115の計算を実行する。
M-step: 更新パラメータを次式で算出する。
なお、"arg max"とは最大値を与える変数を意味する記号である。すなわち、次式は、Q(α)(x|y、φ)(ψ|φ)(xとyは太字であり、(x|y、φ)はQの下付き)を変数ψの関数と考えたとき、その最大値を与える変数の値をψ*とする、という意味である。
U-step: φをψ*によって置き換え、収束することを確認する。収束しない場合は、E-stepに戻って繰り返しが繰返される。
[Alpha-GEMアルゴリズム]
これは、上記M-stepを次式のようにQ関数を非負とするψ+の算出に置き換えたアルゴリズムである。
まず、alpha-HMMアルゴリズムの近似バージョンがalpha-GEMアルゴリズムであることに着目すべきである。数114を有するalpha-EMアルゴリズムの長所によるalpha-HMM推定アルゴリズムの基本的な特性について説明する。
完全なデータをx=(s,y)(x、s、yは太字)とし、s(sは太字)を消失データとし、y(yは太字)を不完全なデータとし、数115は次式に相当するものとする。
もし、上式の値が負にならなければ、以下の不等式が成り立つ。
ここで、alpha-HMMの抽象バージョンに留意すべきである。
・単一配列のalpha-HMM
Non-Causalな更新方程式
離散型のアルファベット系列y(yは太字)の場合において、最大化されるQ関数は次式で示される。なお、lはアルゴリズムの繰返しの指標(インデックス)である。一方、上記文章中に出てくるtはデータの順番に対応する指標(インデックス)である。したがって、lとtは区別されるものである。
まず、状態遷移確率aijの更新方程式について説明する。状態遷移確率aij は、その更新後に確率量とする必要があるので、ラグランジュの未定係数法(Lagrange multiplier)を使用する必要がある。したがって、次式のように最大値を算出するための微分を行なう。
そうすると、上式から次式が得られる。
ここで、Nij(s)(sは太字)は、iからjへの状態遷移の数と位置である。次のステップでは、ソフトウェアで実行可能なアルゴリズムを目指すため、以下の問題を解決する必要がある。
(a)Non-causalityの存在: 右辺にθl+1(θは太字)が含まれている。このままでは、自己撞着性のため計算できない。
(b)右辺の算出には、状態数がNでデータ数がTのとき、O(NT)すなわちNのT乗のオーダーの演算が必要である。
上記の問題は、次のセクションで解決されることとなるが、その前に、出力確率bjk|θl+1(l+1はθの添え字)と初期状態確率πi|θl+1(l+1はθの添え字)の2つの更なる更新式について次式に示す。なお、Njk(s)(sは太字)は、状態st=jでの出力yt=kの発生回数である。
Causal近似と拡張系列: 離散出力の場合
更新方程式である数122〜数124のコアとなる部分は、繰返し指標(インデックス)のシフトによるCausal近似によって、次式のように変換することができる。
したがって、P(y|θl)=P(y|θl-1)+o(1)(yは太字)の領域で次の等式を得る。
ここで、αcausalという表記は、以後、βと表記する。
いま、数125が算出可能であるが、演算の複雑度のために、もう一つの近似が必要となる。このため、我々は次式の系列拡張を用いている。
そして、causal近似式である数125の適用と数127の系列拡張は、以下に示す遷移確率の更新方程式を与える。なお、これらの確率の更新式である上記数122〜数124は、それぞれ上記実施例1の数75〜77に対応するものである。
単一配列の離散型alpha-HMMの遷移確率aij|θl+1(l+1は、θの添え字)は、次式である。
ここで、以下の各特性について理解しておくことが重要である。
[特性1]
β=αcausal=1の場合は、従来のlog-HMMの方法に帰趨する。
[特性2]
数128の分子は、現在と過去の更新項の重み付けられた総和である。分母も同様である。
[特性3]
数128の第2行及び第3行は、これらの確率計算に伴う複雑度を抑えた従来の前方-後方アルゴリズムの方法に匹敵する。
[特性4]
alpha-HMMに唯一、追加する必要があるものは、θl-1での更新項を記憶することである。これは、CPU時間に直接現れているように繰返しを減らすことを意味する。実際の実験結果でも、この予想通りになっている。
次に、残りの2つの更新方程式を示す。
単一配列の離散型alpha-HMMの出力確率bjk|θl+1(l+1は、θの添え字)は、次式である。
単一配列の離散型alpha-HMMの初期状態確率πi|θl+1(l+1は、θの添え字)は、次式である。
(変形例)
以下、近似式の変形例について説明する。なお、以下の変形例は本発明の他の実施例でも同様に適用可能なものである。上述した隠れマルコフモデルの推定法では、有限個の過去情報を用いて各確率量をすべて次のような計算に基づいて計算している。
この高速化パラメータβは、利用者が選ぶデザインパラメータであり、現在値による繰り返し値に対し、1回前の繰り返し値を重みづけする効果を有している。この過去の値に重みづけを行って加算することは、有限の過去にまでさかのぼって行うこともできる。すなわち、ετ≧0をデザインパラメータとして設定し、次式のように置換してもよい。
そうすると、数131は次式となる。
次に展開による近似式を用いる場合について説明すると、上記数131や数134は、次式のように展開できる(例えばテーラー展開)。
ここで、o(1)は高位の無限小を意味する記号であり、上式の右辺第一項は,従来の隠れマルコフモデル推定アルゴリズムに相当している。そして、モーメンタム項に相当する右辺第二項は、展開により近似される高速化項の意味をもつ。
なお、T=1かつε1=P(y|θl-1)/P(y|θl)(yはそれぞれ太字),すなわち重みを用いた場合は、重みを用いていない場合に比べて、収束速度はわずかに遅くなるものの、β=3付近での収束性能の安定性が良くなる。換言すると、図3のような性能曲線において、計算が収束する立ち上がり部分では、高速化に伴い、波状の上下の変動が生じやすくなるが、その変動が少なくなる。
・複数配列のalpha-HMM
もし、HMMの予め設計されたトポロジーは、エルゴードなものであれば、単一の長いトレーニング系列y(yは太字)は、十分である。もし、選択されたトポロジーが、吸収状態を有しているならば、複数のトレーニング系列を使うことが望まれる。このような複数系列に対するalpha-HMM推定の更新式は、上記の単一系列に対する方法を利用する形で得ることができる。
離散的なシンボルの場合において、S(Sは太字)を次式のようにM個の状態遷移系列の集合とする。
すると、複数配列のQ関数は次式である。
ここで、Pは、次式に示されるMarkov過程の確率である。
上式において、初期状態確率πs0(n)(0はsの下付きの添え字であり、(n)は上付きの添え字である)、状態遷移確率ast-1(n)st(n) (t-1とtは、sの下付きの添え字であり、(n)は上付きの添え字である)、及び出力確率bst(n)(yt (n)) (tは、sの下付きの添え字であり、(n)は上付きの添え字である)の形態が、系列指標(インデックス)nから独立している点に着目すべきである。
初期状態確率πi|θl+1(l+1は、θの添え字)によって、初期状態の更新方程式の誘導は、Q関数の数137の微分から始まる。上述の単一配列との違いは、上式からも明らかなように微分がn回現れることである。そして、次式のようにnon-causalな方程式が得られる。
ここで、
であり、上式の分母の関数f(n)は次式である。
そして、causalの繰返し指標(インデックス)のシフトと系列の拡張は、次の更新方程式を与える。
同様に、状態遷移と出力の更新方程式が得られる。
ここで、
及び
また、次式の出力確率が得られる。
ここで、上式で省略された各項はそれぞれ次式となる。
上記数147のΣの下にある記号「t:yt=k」の意味は、yt=kが成り立つtのみを対象にする、という意味である。すなわち、上記数147の場合、y={yt}T(t=1)(左辺のyは太字であり、t=1は{yt}の下付きの添え字)において、yt=kが成り立つtのみについて加算する、という意味である。
(B)連続系列の単一配列
ここから、連続的な出力系列の場合の本実施例の具体的な説明を行なう。
もし、出力系列y={yt}T(t=1)(出力系列yは太字であり、t=1は{yt}の下付き)が、連続的な多変数の観測結果として現れれば、数128〜数130と似ているが、少しだけ異なる更新方程式が得られる。この場合、ytは太字の文字で示されていないが、適切な次元のユークリッド空間におけるベクトルである。
このような連続的なアルファベットの場合、最尤推定法(MLE)の問題は、以下の尤度を最大化することである。
ここで、状態stで、出力確率bstkt(yt)(stktのtは、それぞれs,kの下付きの添え字)において、kt番目の枝へ遷移する確率を特定する分岐確率Cstkt(stktのtは、それぞれs,kの下付きの添え字)は、ytの確率密度関数である。我々は、これをガウス密度(Gaussian density)とみなしている。そして、状態jで状態kが出力される出力確率bjk(yt)は次式となる。
ここで、μjkは平均値ベクトルであり、Σjkは共分散行列であり、総和と混同してはいけない。平均値ベクトルとは、分岐確率の確率密度関数の平均値を示している(以後、同様)。そして、状態jでの出力確率密度関数(pdf)は、bj(yt)がpdfとなるように次式となる。
ここで、上式の記号について説明する。
まず、離散系列の場合、kを出力記号の種類を表すインデックスとすると、出力ykのとる値は有限種類であるため、bj(yk)はbjkと略記できる。しかし、連続系列の場合は、離散値ではなく、連続値であるため、そのような表記はできない。そのため、時間を表すインデックスをtとして、ytはそのまま記載する必要がある。
また、N(yt;μjk,Σjk)は、平均値ベクトルがμjkであり、Σjkを共分散行列とする多次元正規確率密度関数(多次元ガウス確率密度関数)で、ytを分布の変数とするという意味である。
実際、このようなガウス混合モデルは、「L. A Liporace, “Maximum likelihood estimation for multivariate observations of Markov sources” IEEE Trans. IT, vol. 28, pp. 729-734, 1982.」、及び「B.-H. Juang, “Maximum-likelihood estimation for mixture multivariate stochastic observations of Markov chains,” AT & T Tech. J., vol. 64, pp. 1235-1245, 1985.」等にも記載されているモデルであり、log-HMMでもここまでが導出可能な場合にとどまっている。
図6は、離散的なアルファベットの場合を図示している。また、図7は、ガウス混合モデル(Gaussian mixture model)(bottom)の場合、すなわち連続的なアルファベットの場合を図示している。これらの図面を参照すると、図7中の分岐確率cjkの矢印は、図6中の出力確率bjkの矢印と対応していることがわかる。 数149と数100を参照すると、数149の分岐確率cstkt(stktのtは、それぞれs,kの下付きの添え字)と出力確率bstkt(yt)(stktのtは、それぞれs,kの下付きの添え字)の積が、数100の出力確率bst(yt)(stのtは、sの下付きの添え字)に対応しているように見えるかもしれない。しかし、数149において離散値シンボルとみなされるのは、出力確率bstkt(yt)(stktのtは、それぞれs,kの下付きの添え字)を除いた部分であるため、図6、図7の対応関係との矛盾は生じない。
混合確率の場合のため、消失データはs(sは太字)とc(cは太字)である。そのため、Q関数は次式である。
数149、数152から明らかなように、初期確率と状態遷移確率の更新方程式は、それぞれ数128、数130と同様である。分岐確率cjkの更新方程式は、状態遷移確率aijの場合と同様にラグランジュの未定係数法(Lagrange multiplier)によって得られる。
そして、右辺は、状態遷移確率の場合と同様にcausalで計算可能となり、次式のように変形される。
次式は、平均値ベクトルμjkの更新方程式である。μjkについて数152の直接微分として、次のnon-causalの方程式が得られる。
そして、繰返し指標(インデックス)のシフト、系列の拡張、総和の変更により、次式が得られる。
数156が、θl-1の過去情報が十分に利用できる点を示していることに着目すべきである。
共分散行列の更新には、行列微分が必要である。Q関数の数152を共分散行列の逆行列Σ-1 jk(jkはΣの添え字)について微分することにより、次のnon-causalな方程式が得られる。
そして、繰返し指標(インデックス)のシフト、系列の拡張、および総和の変更により、次の更新式が得られる。
ここで、上式の各項は以下の通りである。
共分散行列の更新式である数158には、効果的に十分に利用できる過去情報の形態を有していることに着目すべきである。なお、記号Σ-1 jk|θl;1とは、l+1回目の状態jで分岐kに移行したときの出力の共分散行列を表している(以後、同様)。
ここで、ガウス混合(Gaussian mixture)alpha-HMM、すなわち単一系列の連続型alpha-HMMの更新方法について簡単に述べる。
[単一系列の連続型alpha-HMMの初期状態確率]
更新式は、数130である。
[単一系列の連続型alpha-HMMの状態遷移確率]
更新式は、数128である。
[単一系列の連続型alpha-HMMの分岐確率]
更新式は、数154である。
[単一系列の連続型alpha-HMMの平均値ベクトル]
更新式は、数156である。
[単一系列の連続型alpha-HMMの共分散行列]
更新式は、数158であり、その要素は、数159、数160である。
ここで、θl-1のすべての情報をメモリに記憶する点を改めて強調することが重要である。θlによって指標(インデックス)を付された項の算出は、log-HMMと同等である。
図8は、本発明の第3実施例におけるプログラムの処理手順をあらわしたものである。このフローチャート全体の流れは、前記第1実施例の場合とほぼ同様である。ただし、本実施例では、観測データyが、連続的な多変数の観測結果として現れる単一配列の連続系列データである点で異なる。また、初期値として設定され、繰返し計算される未知パラメータには、初期状態確率、状態遷移確率、分岐確率、平均値ベクトル、共分散行列があるため、各ステップにおける計算式も異なる。
そのため、置き換わる数式の対応関係について説明する。まず、繰返し更新される未知パラメータの集合の組み合わせを次式とする。
推定装置10が最終的に算出しようとするHMMの確率構造は、数66の代わりに数149となる。ステップS32で決定する未知パラメータの初期値は、数67の代わりに次式となる。
ステップS32における収束判定値の決定方法は、数83と同じである。
ステップS33では、最初にステップS32で決定した未知パラメータの初期値を利用し、それ以降はステップS39で更新された未知パラメータとカウント値を利用して、次の数163に示す各確率値と、数164に示す各カウント値を、実際に推定装置10のメモリ(図示せず)に設定する。なお、前記実施例1、2と異なり、指標(インデックス)は、データの順番のtではなく、アルゴリズムの繰り返しのlを付して設定する。
その後のステップS34〜S36の演算動作は、前記実施例1の数70〜数74の一部の符号を置き換えて適用することによって行う。具体的には、出力確率bj.yτ+1(τ+1はyの下付き添え字)を分岐確率cj.yτ+1(τ+1はyの下付添え字)と置き換え、指標(インデックス)であるtをlと置き換える。
ステップS37では、l+1回目に繰り返される確率π,a,cと、平均値ベクトルμ,及び共分散行列Σの組み合わせを条件として、状態iとなる初期状態確率
初期状態確率πi|θl+1(l+1はθの添え字)の値と、状態iから状態jに移る状態遷移確率aij|θl+1(l+1はθの添え字)の値と、状態jで状態kが出力される出力確率bjk(yt)において、分岐kに移行する枝へ遷移する確率を特定する分岐確率cjk|l+1(l+1はθの添え字)の値と、平均値ベクトルμjk|θl+1(l+1はθの添え字)の値と,共分散行列Σjk|θl+1(l+1はθの添え字)の値と,ステップS34で求めたl回目の確率P(y|θl)と、それよりも1回前の確率P(y|θl-1)を利用し、且つステップS35で求めたl回目のカウント値Naij|θl,Ncjk|θl(ijはaの添え字、jkはcの添え字、lはθの添え字)と、それよりも1回前のカウント値Naij|θl-1,Ncjk|θl-1(ijはaの添え字、jkはcの添え字、l-1はθの添え字)を利用して計算する部分であり、これは前記数128、数130、数154、数156、数158〜数160の各更新式にそれぞれ対応している。
このステップS37においても高速化パラメータの値βが用いられ、添え字l+1の部分の確率を計算する上で、時間シフトした添え字lや添え字l-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み(高速化パラメータの値)が1≦β<3に拡張したこと等は、前記実施例1と同様にその効果として挙げられる。
本実施例では、配列が1本の観測データをレジスタ1に格納しているが、その場合のステップS37における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数130、状態遷移確率の更新式は数128、分岐確率の更新式は数154、平均値ベクトルの更新式は数156、共分散行列の更新式は数158〜数160である。
ステップS38では、前記実施例1と同様に、前記ステップ37で計算された新たな確率量に基づく尤度Pを用いて、その尤度PとステップS32における収束判定値との比較により、ステップS33〜S39の計算が収束したか否かを判定する。収束してなければ、ステップS39に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように2つの過去値をシフトさせる。
ここでは、繰り返しの回数が1つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。なお、ステップS33とステップS39を一纏めにして、更新された確率量およびカウント値をそのまま次の繰り返しのために、推定装置10のメモリに設定記憶させてもよい。
一方、ステップS38において計算が収束していると判定された場合、ステップS40に移行して、次式に示す計算された5組の確率量を用いたHMMを採用する。
上記数166で計算した各値を用いて、推定装置10は数149に示すHMMの確率構造を出力することができる。
図8に示すアルゴリズムを実現する推定装置の構成は、図2に示した前記第1実施例の推定装置10と同様であるため、ここでは説明を省略する。
以上のように本実施例では、入力される観測データを時系列に格納する記憶手段としてのレジスタ1と、観測データがどのような確率モデルであるのかを、HMMの未知パラメータを算出することで推定する推定手段14とを備えたHMMの推定装置10において、推定手段14は、HMMの高速化パラメータの値βを設定する初期設定手段22と、前記HMMの未知パラメータとして、状態遷移確率a,分岐確率c,平均値ベクトルμ,共分散行列Σ,初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび分岐の期待値Ncを更新設定する更新設定手段24と、更新設定手段24で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる(数163,数164を参照)と共に、前記レジスタ1から読み出した観測データと、初期設定手段22で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する(数128,数130,数154,数156,数158〜数160を参照)演算手段26と、演算手段26による計算の収束を判定し、計算が収束していなければ、演算手段26で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段26で計算した新たな各確率量を最終的な値として出力させる判定手段28とを備えている。
このようにすれば、HMMの未知パラメータとして、状態遷移確率a,分岐確率c,平均値ベクトルμ,共分散行列Σ,初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび分岐の期待値Ncを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なHMM推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
また、ここでの演算手段26は、l+1回目の状態iとなる初期状態確率πi|θl+1(l+1はθの添え字)の値を、前記数130の式で計算し、l+1回目の状態iから状態jに移る状態遷移確率aij|θl+1(l+1はθの添え字)の値を、l回目およびl-1回目における状態遷移の期待値Naij|θlおよびNaij|θl-1(ijはaの添え字、lまたはl-1はθの添え字)を利用して、前記数128の式で計算し、状態jで状態kが出力される出力確率bjk|θl+1(l+1はθの添え字)の際に分岐kに移行する枝へ遷移する確率を特定する分岐確率cjk|θl+1の値と、平均値ベクトルμjk|θl+1(l+1はθの添え字)の値と、共分散行列Σjk|θl+1の値とを、l回目およびl-1回目における前記出力の期待値Ncjk|θlおよびNcjk|θl-1(jkはcの添え字、lまたはl-1はθの添え字)を利用して、前記数154、数156、数158〜数160の式で計算する構成となっている。
そのため、特に高速化パラメータの値βが1でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
さらに、ここでの更新設定手段24は、高速化パラメータの値を1<β<3に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが3未満の範囲で、βを1よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
なお上述した本実施例の作用効果は、初期設定手段22としての動作を実行する初期設定ステップと、更新設定手段24としての動作を実行する更新設定ステップと、演算手段26としての動作を実行する演算ステップと、判定手段28としての動作を実行する判定ステップとを備えたHMMの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段14として、コンピュータに機能させるHMMの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。
本実施例では、連続的なシンボルで複数配列のデータ系列の場合の実施形態について説明する。
本実施例においても、上記実施例と同様に、繰返し更新される未知パラメータの集合の組み合わせを数161とし、連続的なシンボルのalpha-HMMの更新方程式が得られる。初期状態確率と状態遷移確率の更新方程式は、それぞれ数142〜数145と同様である。しかし、出力の更新方程式の集合は、離散的なシンボルの場合と異なる。我々は、分岐確率、平均値ベクトル、及び共分散行列の更新方程式を必要としている。
分岐確率については、初期状態確率と状態遷移確率の場合と同様にラグランジュの未定係数法を用いて算出可能である。そして、次式の更新方程式が得られる。
ここで、各要素は次式である。
Q関数の数137での直接ベクトル微分であるμjk|θl+1(l+1はθの下付きの添え字)が適用される。その更新方程式は次式である。
ここで、
数170、数171と同様に、共分散行列の更新方程式は、Σ-1 jk|θl(jk|θlはΣの下付き添え字であり、さらにlはθの下付き添え字)について行列微分を用いることによって得られる。
ここで、G(n) θl(θlはGの下付き添え字)は次式である。
図9は、本発明の第4実施例におけるプログラムの処理手順をあらわしたものである。このフローチャート全体の流れは、前記第3実施例の場合とほぼ同様である。ただし、本実施例では、観測データyが、連続的な多変数の観測結果として現れる複数配列の連続系列データである点で異なる。
置き換わる数式の対応関係について説明する。まず、推定装置100が最終的に算出しようとするHMMの確率構造は、数66の代わりに数149となる。ステップS42で決定する未知パラメータの初期値は、数162と同じである。
ステップS42における収束判定値の決定方法は、数83と同じである。
ステップS43では、最初にステップS42で決定した未知パラメータの初期値を利用し、それ以降はステップS49で更新された未知パラメータとカウント値を利用して、数174に示す各確率値と、数175に示す各カウント値を、実際に推定装置100のメモリ(図示せず)に設定する。なお、指標(インデックス)は、データの順番のtではなく、アルゴリズムの繰り返しのlを付して設定する。
その後のステップS44〜S46の演算動作は、前記実施例2の数88〜数92の一部の符号を置き換えて適用することによって行う。具体的には、出力確率bj.yτ+1(τ+1はyの下付添え字)を分岐確率cj.yτ+1(τ+1はyの下付添え字)と置き換え、指標(インデックス)であるtをlと置き換える。
ステップS47では、l+1回目に繰り返される確率π,a,cと、平均値ベクトルμ,及び共分散行列Σの組み合わせを条件として、状態iとなる初期状態確率πi|θl+1(l+1はθの添え字)の値と、状態iから状態jに移る状態遷移確率aij|θl+1(l+1はθの添え字)の値と、状態jで状態kが出力される出力確率bjk(yt)において、分岐kに移行する枝へ遷移する確率を特定する分岐確率cjk|l+1(l+1はθの添え字)の値と、平均値ベクトルμjk|θl+1(l+1はθの添え字)の値と,共分散行列Σjk|θl+1(l+1はθの添え字)の値とを、ステップS34で求めたl回目の確率P(y(n)|θl)と、それよりも1回前の確率P(y(n)|θl-1)を利用し、且つステップS35で求めたl回目のカウント値N(n) aij|θl,N(n) cjk|θl(ijはaの添え字、jkはcの添え字、lはθの添え字)と、それよりも1回前のカウント値N(n) aij|θl-1,N(n) cjk|θl-1(ijはaの添え字、jkはcの添え字、l-1はθの添え字)を利用して計算する部分であり、これは前記数142〜数145、数167〜数173に対応している。
このステップS47においても高速化パラメータの値βが用いられ、添え字l+1の部分の確率を計算する上で、時間シフトした添え字lや添え字l-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み(高速化パラメータの値)が1≦β<3に拡張したこと等は、前記実施例2と同様にその効果として挙げられる。
本実施例では、配列がM本の観測データをレジスタ101に格納しているが、その場合のステップS47における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数142、状態遷移確率の更新式は数143〜数145、分岐確率の更新式は数167〜数169、平均値ベクトルの更新式は数170、数171、共分散行列の更新式は数172、数173である。
ステップS48では、前記実施例2と同様に、前記ステップS47で計算された新たな確率量に基づく尤度Pを用いて、その尤度PとステップS42における収束判定値との比較により、ステップS43〜S49の計算が収束したか否かを判定する。収束してなければ、ステップS49に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように2つの過去値をシフトさせる。
ここでは、繰り返しの回数が1つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。
一方、ステップS48において計算が収束していると判定された場合、ステップS50に移行して、次式に示す計算された5組の確率量を用いたHMMを採用する。
上記数177で計算した各値を用いて、推定装置101は数149に示すHMMの確率構造を出力することができる。
図9に示すアルゴリズムを実現する推定装置の構成は、図5に示した前記第2実施例の推定装置101と同様であるため、ここでは説明を省略する。
以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ101と、観測データがどのような確率モデルであるのかを、HMMの未知パラメータを算出することで推定する推定手段114とを備えたHMMの推定装置100において、推定手段114は、HMMの高速化パラメータの値βを設定する初期設定手段122と、前記HMMの未知パラメータとして、状態遷移確率a,分岐確率c,平均値ベクトルμ,共分散行列Σ,初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび分岐の期待値Ncを更新設定する更新設定手段124と、更新設定手段124で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる(数174,数175を参照)と共に、前記レジスタ101から読み出した観測データと、初期設定手段122で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する(数142〜数145,数167〜数173を参照)演算手段126と、演算手段126による計算の収束を判定し、計算が収束していなければ、演算手段126で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段126で計算した新たな各確率量を最終的な値として出力させる判定手段128とを備えている。
このようにすれば、HMMの未知パラメータとして、状態遷移確率a,分岐確率c,平均値ベクトルμ,共分散行列Σ,初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび分岐の期待値Ncを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なHMM推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
また、ここでの演算手段126は、l+1回目の状態iとなる初期状態確率πi|θl+1(l+1はθの添え字)の値を、前記数142の式で計算し、l+1回目の状態iから状態jに移る状態遷移確率aij|θl+1(t+1はθの添え字)の値を、l回目およびl-1回目における状態遷移の期待値Naij|θlおよびNaij|θl-1(ijはaの添え字、lまたはl-1はθの添え字)を利用して、前記数143〜数145の式で計算し、状態jで状態kが出力される出力確率bjk|θl+1(l+1はθの添え字)の際に分岐kの枝へ遷移する確率を特定する分岐確率cjk|θl+1の値と、平均値ベクトルμjk|θl+1(l+1はθの添え字)の値と、共分散行列Σjk|θl+1の値とを、l回目およびl-1回目における前記分岐の期待値Ncjk|θlおよびNcjk|θl-1(jkはcの添え字、lまたはl-1はθの添え字)を利用して、前記数167〜数173の式で計算する構成となっている。
そのため、特に高速化パラメータの値βが1でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
さらに、ここでの更新設定手段124は、高速化パラメータの値を1<β<3に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが3未満の範囲で、βを1よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
なお上述した本実施例の作用効果は、初期設定手段122としての動作を実行する初期設定ステップと、更新設定手段124としての動作を実行する更新設定ステップと、演算手段126としての動作を実行する演算ステップと、判定手段128としての動作を実行する判定ステップとを備えたHMMの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段114として、コンピュータに機能させるHMMの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。
本実施例では、半連続的なシンボルで単一配列のデータ系列の場合の実施形態について説明する。
図7において、ガウス混合(Gaussian mixture)alpha-HMM(log-HMMも)のグラフ構造を見直すことにより、以下を実現する。
(a)ガウス混合HMMにおいて、各ガウス(Gaussian)pdfは、到着する状態jに依存する。すべてのN×Kのガウス密度の学習には、多様な長いトレーニング系列を必要とする。
(b)離散型の場合のbjkの役割を連続モードのcjkに割り当てる。また、平均値ベクトルと共分散行列は遷移状態jに依存しないという場合を考慮する。そして、図6の離散型の場合をまっすぐ延長した構造となる。これを半連続HMMと呼んでいる。この構造のもう一つの解釈は、MLE-VQ HMM (Maximum Likelihood Vector Quantization HMM)である。
半連続のalpha-HMMモデルは、数149を変形し、次式となる。
したがって、半連続のalpha-HMMの更新方程式は以下のようになる。
[単一配列の半連続のalpha-HMMの初期状態確率]
更新式は、数130である。
[単一配列の半連続のalpha-HMMの状態遷移確率]
更新式は、数128である。
[単一配列の半連続のalpha-HMMの分岐確率]
更新式は、数154である。
[単一配列の半連続のalpha-HMMの平均値ベクトル]
次式のように、数156において、μjk|θl+1(l+1はθの下付きの添え字)をμj|θl+1(l+1はθの下付きの添え字)とし、右辺のkt=kの項を除去したものである。これは、状態jに依存しないためである(以後、同様である)。
[単一配列の半連続のalpha-HMMの共分散行列]
次式のように、数158〜数160において、Σjk|θl+1(l+1はθの添え字),μjk|θl(lはθの添え字),μjk|θl-1(l-1はθの添え字)を、それぞれΣj|θl+1(l+1はθの添え字),μj|θl(lはθの添え字),μj|θl-1(l-1はθの添え字)とし、右辺のkt=kの項を除去したものである。
ここで、各要素は次式である。
本発明の第5実施例におけるプログラムの処理手順は、前記第3実施例の場合とほぼ同様であるため、図8を用いて説明する。ただし、本実施例では、観測データyが、連続的な多変数の観測結果として現れる単一配列の半連続系列データである点で異なる。
本実施例は、前記第3実施例の場合の特例と考えることもできるので、置き換わる数式の対応関係に絞って説明する。まず、繰返し更新される未知パラメータの集合の組み合わせを数161とし、推定装置10が最終的に算出しようとするHMMの確率構造は、数66の代わりに数178とする。ステップS32で決定する未知パラメータの初期値は、数67の代わりに次式となる。
ステップS32における収束判定値の決定方法は、数83と同じである。
ステップS33では、最初にステップS32で決定した未知パラメータの初期値を利用し、それ以降はステップS39で更新された未知パラメータとカウント値を利用して、次の数184に示す各確率値と、数185に示す各カウント値を、実際に推定装置10のメモリ(図示せず)に設定する。
その後のステップS34〜S36の演算動作は、前記実施例3と同様に前記実施例1の数70〜数74の一部の符号を置き換えて適用することによって行なう。
ステップS37では、l+1回目に繰り返される確率π,a,cと、平均値ベクトルμ,及び共分散行列Σの組み合わせを条件として、状態iとなる初期状態確率πi|θl+1(l+1はθの添え字)の値と、状態iから状態jに移る状態遷移確率aij|θl+1(l+1はθの添え字)の値と、状態jで状態kが出力される出力確率bjk(yt)において、分岐kに移行する枝へ遷移する確率を特定する分岐確率cjk|l+1(l+1はθの添え字)の値と、平均値ベクトルμj|θl+1(l+1はθの添え字)の値と,共分散行列Σj|θl+1(l+1はθの添え字)の値とを、ステップS34で求めたl回目の確率P(y|θl)と、それよりも1回前の確率P(y|θl-1)を利用し、且つステップS35で求めたl回目のカウント値Naij|θl,Ncjk|θl(ijはaの添え字、jkはcの添え字、lはθの添え字)と、それよりも1回前のカウント値Naij|θl-1,Ncjk|θl-1(ijはaの添え字、jkはcの添え字、l-1はθの添え字)を利用して計算する部分であり、これは前記数130、数128、数154、数179〜数182に対応している。
このステップS37においても高速化パラメータの値βが用いられ、添え字l+1の部分の確率を計算する上で、時間シフトした添え字lや添え字l-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み(高速化パラメータの値)が1≦β<3に拡張したこと等は、前記実施例3と同様にその効果として挙げられる。
本実施例では、配列が1本の観測データをレジスタ1に格納しているが、その場合のステップS37における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数130、状態遷移確率の更新式は数128、分岐確率の更新式は数154、平均値ベクトルの更新式は数179、共分散行列の更新式は数180〜数182である。
ステップS38では、前記実施例3と同様に、前記ステップ37で計算された新たな確率量に基づく尤度Pを用いて、その尤度PとステップS32における収束判定値との比較により、ステップS33〜S39の計算が収束したか否かを判定する。収束してなければ、ステップS39に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように2つの過去値をシフトさせる。
ここでは、繰り返しの回数が1つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。なお、ステップS33とステップS39を一纏めにして、更新された確率量およびカウント値をそのまま次の繰り返しのために、推定装置10のメモリに設定記憶させてもよい。
一方、ステップS38において計算が収束していると判定された場合、ステップS40に移行して、次式に示す計算された5組の確率量を用いたHMMを採用する。
上記数187で計算した各値を用いて、推定装置10は数178に示すHMMの確率構造を出力することができる。
図8に示すアルゴリズムを実現する推定装置の構成は、図2に示した前記第3実施例の推定装置10と同様であるため、ここでは説明を省略する。
以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ1と、観測データがどのような確率モデルであるのかを、HMMの未知パラメータを算出することで推定する推定手段14とを備えたHMMの推定装置10において、推定手段14は、HMMの高速化パラメータの値βを設定する初期設定手段22と、前記HMMの未知パラメータとして、状態遷移確率a,分岐確率c,平均値ベクトルμ,共分散行列Σ,初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび分岐の期待値Ncを更新設定する更新設定手段24と、更新設定手段24で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる(数184,数185を参照)と共に、前記レジスタ1から読み出した観測データと、初期設定手段22で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する(数130,数128,数154,数179〜数182を参照)演算手段26と、演算手段26による計算の収束を判定し、計算が収束していなければ、演算手段26で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段26で計算した新たな各確率量を最終的な値として出力させる判定手段28とを備えている。
このようにすれば、HMMの未知パラメータとして、状態遷移確率a,分岐確率c,平均値ベクトルμ,共分散行列Σ,初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび分岐の期待値Ncを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なHMM推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
また、ここでの演算手段26は、l+1回目に状態iとなる初期状態確率πi|θl+1(l+1はθの添え字)の値を、前記数130の式で計算し、l+1回目の状態iから状態jに移る状態遷移確率aij|θl+1(t+1はθの添え字)の値を、l回目およびl-1回目における状態遷移の期待値Naij|θlおよびNaij|θl-1(ijはaの添え字、lまたはl-1はθの添え字)を利用して、前記数128の式で計算し、状態jで状態kが出力される出力確率bjk|θl+1(l+1はθの添え字)の際に分岐kに移行する枝へ遷移する確率を特定する分岐確率cjk|θl+1の値と、平均値ベクトルμjk|θl+1(l+1はθの添え字)の値と、共分散行列Σjk|θl+1の値とを、l回目およびl-1回目における前記出力の期待値Ncjk|θlおよびNcjk|θl-1(jkはcの添え字、lまたはl-1はθの添え字)を利用して、前記数154、数179〜数182の式で計算する構成となっている。
そのため、特に高速化パラメータの値βが1でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
さらに、ここでの更新設定手段24は、高速化パラメータの値を1<β<3に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが3未満の範囲で、βを1よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
なお上述した本実施例の作用効果は、初期設定手段22としての動作を実行する初期設定ステップと、更新設定手段24としての動作を実行する更新設定ステップと、演算手段26としての動作を実行する演算ステップと、判定手段28としての動作を実行する判定ステップとを備えたHMMの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段14として、コンピュータに機能させるHMMの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。
本実施例では、半連続的なシンボルで複数配列のデータ系列の場合の実施形態について説明する。
この場合の複数配列の更新方程式は、平均ベクトルと共分散の状態依存性を制限することによって得られる。初期状態確率と状態遷移確率と分岐確率の更新方程式は、数142〜数145、および数167〜数169と同様である。したがって、次式のようになる。ここで、各記号の上にバー(-)を付し、バー(-)のない場合と記号を区別している。これは、例えば、数171と以下の数189では数式が異なっており、同じ記号が使えないためである。すなわち、バー(-)それ自体に特別な意味はない。
また、本実施例では、前記第5実施例において添え字がjだった記号は、添え字がkになっている。これは、kt=kの項の代わりにst=jの項が除去されたためである。(以後、同様である)。
上記数188と同様に、共分散行列の更新方程式は、状態依存性を除去することによって次式のように得られる。
本発明の第6実施例におけるプログラムの処理手順は、前記第4実施例の場合とほぼ同様であるため、図9を用いて説明する。ただし、本実施例では、観測データyが、連続的な多変数の観測結果として現れる複数配列の半連続系列データである点で異なる。
本実施例は、前記第4実施例の場合の特例と考えることもできるので、置き換わる数式の対応関係に絞って説明する。まず、繰返し更新される未知パラメータの集合の組み合わせを数161とし、推定装置100が最終的に算出しようとするHMMの確率構造は、数66の代わりに数178とする。ステップS42で決定する未知パラメータの初期値は、数182と同じである。
ステップS42における収束判定値の決定方法は、数83と同じである。
ステップS43では、最初にステップS42で決定した未知パラメータの初期値を利用し、それ以降はステップS49で更新された未知パラメータとカウント値を利用して、次の数192に示す各確率値と、数193に示す各カウント値を、実際に推定装置100のメモリ(図示せず)に設定する。
その後のステップS44〜S46の演算動作は、前記実施例4と同様に前記実施例2の数88〜数92の一部の符号を置き換えて適用することによって行なう。
ステップS47では、l+1回目に繰り返される確率π,a,cと、平均値ベクトルμ,及び共分散行列Σの組み合わせを条件として、状態iとなる初期状態確率πi|θl+1(l+1はθの添え字)の値と、状態iから状態jに移る状態遷移確率aij|θl+1(l+1はθの添え字)の値と、状態jで状態kが出力される出力確率bjk(yt)において、分岐kに移行する枝へ遷移する確率を特定する分岐確率cjk|l+1(l+1はθの添え字)の値と、平均値ベクトルμ-k|θl+1(l+1はθの添え字であり、-はμの上に付く)の値と,共分散行列Σk|θl+1(l+1はθの添え字であり、-はΣの上に付く)の値とを、ステップS34で求めたl回目の確率P(y(n)|θl)と、それよりも1回前の確率P(y(n)|θl-1)を利用し、且つステップS35で求めたl回目のカウント値N(n) aij|θl,N(n) cjk|θl(ijはaの添え字、jkはcの添え字、lはθの添え字)と、それよりも1回前のカウント値N(n) aij|θl-1,N(n) cjk|θl-1(ijはaの添え字、jkはcの添え字、l-1はθの添え字)を利用して計算する部分であり、これは前記数142〜数145、数167〜数169、数188〜数191に対応している。
このステップS47においても高速化パラメータの値βが用いられ、添え字l+1の部分の確率を計算する上で、時間シフトした添え字lや添え字l-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み(高速化パラメータの値)が1≦β<3に拡張したこと等は、前記実施例4と同様にその効果として挙げられる。
本実施例では、配列がM本の観測データをレジスタ101に格納しているが、その場合のステップS47における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数142、状態遷移確率の更新式は数143〜数145、分岐確率の更新式は数167〜数169、平均値ベクトルの更新式は数188、数189、共分散行列の更新式は数190、数191である。
ステップS48では、前記実施例4と同様に、前記ステップS47で計算された新たな確率量に基づく尤度Pを用いて、その尤度PとステップS42における収束判定値との比較により、ステップS43〜S49の計算が収束したか否かを判定する。収束してなければ、ステップS49に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように2つの過去値をシフトさせる。
ここでは、繰り返しの回数が1つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。
一方、ステップS48において計算が収束していると判定された場合、ステップS50に移行して、次式に示す計算された5組の確率量を用いたHMMを採用する。
上記数195で計算した各値を用いて、推定装置101は数178に示すHMMの確率構造を出力することができる。
図9に示すアルゴリズムを実現する推定装置の構成は、図5に示した前記第2実施例の推定装置101と同様であるため、ここでは説明を省略する。
以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ101と、観測データがどのような確率モデルであるのかを、HMMの未知パラメータを算出することで推定する推定手段114とを備えたHMMの推定装置100において、推定手段114は、HMMの高速化パラメータの値βを設定する初期設定手段122と、前記HMMの未知パラメータとして、状態遷移確率a,分岐確率c,平均値ベクトルμ-(-はμの上に付く),共分散行列Σ-(-はΣの上に付く),初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび分岐の期待値Ncを更新設定する更新設定手段124と、更新設定手段124で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる(数192,数193を参照)と共に、前記レジスタ101から読み出した観測データと、初期設定手段122で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する(数142〜数145,数167〜数169,数188〜数191を参照)演算手段126と、演算手段126による計算の収束を判定し、計算が収束していなければ、演算手段126で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段126で計算した新たな各確率量を最終的な値として出力させる判定手段128とを備えている。
このようにすれば、HMMの未知パラメータとして、状態遷移確率a,分岐確率c,平均値ベクトルμ-(-はμの上に付く),共分散行列Σ-(-はμの上に付く),初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび分岐の期待値Ncを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なHMM推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
前記数167〜数173の式で計算する構成となっている。
また、ここでの演算手段126は、l+1回目の状態iとなる初期状態確率πi|θl+1(l+1はθの添え字)の値を、前記数142の式で計算し、l+1回目の状態iから状態jに移る状態遷移確率aij|θl+1(t+1はθの添え字)の値を、l回目およびl-1回目における状態遷移の期待値Naij|θlおよびNaij|θl-1(ijはaの添え字、lまたはl-1はθの添え字)を利用して、前記数143〜数145の式で計算し、状態jで状態kが出力される出力確率bjk|θl+1(l+1はθの添え字)の際に分岐kに移行する枝へ遷移する確率を特定する分岐確率cjk|θl+1の値と、平均値ベクトルμ-k|θl+1(l+1はθの添え字であり、-はμの上に付く)の値と、共分散行列Σ-k|θl+1(l+1はθの添え字であり、-はΣの上に付く)の値とを、l回目およびl-1回目における前記分岐の期待値Ncjk|θlおよびNcjk|θl-1(jkはcの添え字、lまたはl-1はθの添え字)を利用して、前記数167〜数169、数188〜数191の式で計算する構成となっている。
そのため、特に高速化パラメータの値βが1でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
さらに、ここでの更新設定手段124は、高速化パラメータの値を1<β<3に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが3未満の範囲で、βを1よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
なお上述した本実施例の作用効果は、初期設定手段122としての動作を実行する初期設定ステップと、更新設定手段124としての動作を実行する更新設定ステップと、演算手段126としての動作を実行する演算ステップと、判定手段128としての動作を実行する判定ステップとを備えたHMMの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段114として、コンピュータに機能させるHMMの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。
上記実施例1〜6では、alpha-HMMの6つのタイプ、すなわち、{離散的,連続的,半連続的}×{単一配列,複数配列}について述べた。本実施例は、区切りのある連続アルファベットの場合である。これは、離散的な文字と連続的な文字が混在している場合の解釈である。この離散と連続が混在している更新方程式は、数154、数169にある総和を2段の総和に変形し、ゼロ出力の許容によって得られる。例えば、数154の分母の総和は、次式のような2段の総和に分割される。
ここで、Dgは、部分が重複しない集合である。Gは、その部分の濃度である。もし、Dgが連続的なアルファベットの特定のサブクラスに相当するならば、このことは、そのアルファベットが離散的なシンボルを伴うものと見なすことができる。換言すると、データが連続している各区間において、それぞれ総和を取り、さらに当該総和どうしについて総和を取っている。この離散的で連続的なアルファベットの場合にも、上記実施例のように単一配列の場合と複数配列の場合がある。したがって、この申請は、全部で8つのタイプのalpha-HMMを有している。
なお、上式の(・)という記号は、Σの中を省略し、総和の範囲を示すΣだけについて考察するために用いたものである。
この場合には、離散値は連続値を出す集団をグループ化した際のラベルと考えればよい。すなわち、データ系列は次式のように示される。
ここで、ラベルは、次式として考える。
上式のGは、G≦K(分岐枝の数)である。そして、新たにyt-(-はytの上に付く)をytと書き直せば、ΣT t=1(・)(TはΣの直上に付き、t=1は直下に付く)がグループ化に対応する部分を含んでいる場合には,これをΣt∈Dg(・)(t∈DgがΣの直下に付く)に変更するだけでよい.ただし、Dg={t|ct=g}である。
以下、各パラメータの更新式について述べる。
・離散連続混在系列の単数配列の場合
[単一配列の離散連続混在系列のalpha-HMMの初期状態確率]
初期状態確率は数130と同じである。
[単一配列の離散連続混在系列のalpha-HMMの状態遷移確率]
状態遷移確率は数128と同じである。
[単一配列の離散連続混在系列のalpha-HMMの分岐確率]
分岐確率はグループ化情報を併せもち、数154に代わって次式となる。以後、これを単にグループ化確率とよぶ。
[単一配列の離散連続混在系列のalpha-HMMの平均値ベクトル]
各グループに対する平均値ベクトルは数156に代わって次式となる。
[単一配列の離散連続混在系列のalpha-HMMの共分散行列]
各グループに対する共分散行列は、数158〜数160に代わって、次式となる。
各要素は次式である。
・離散連続混在系列の複数配列の場合
[複数配列の離散連続混在系列のalpha-HMMの初期状態確率]
初期状態確率は数142と同じである。
[複数配列の離散連続混在系列のalpha-HMMの状態遷移確率]
状態遷移確率は数143〜数145と同じである。
[複数配列の離散連続混在系列のalpha-HMMの分岐確率]
分岐確率すなわちグループ化確率は、数167〜数169の変形として次式のようになる。
各要素は次式である。
[複数配列の離散連続混在系列のalpha-HMMの平均値ベクトル]
平均値ベクトルは、数170、数171の変形として次式となる。
各要素は次式である。
[複数配列の離散連続混在系列のalpha-HMMの共分散行列]
共分散行列は、数172、数173の変形として次式となる。
各要素は次式である。
本実施例も上記第3〜6実施例と同様に、符号や未知パラメータの更新式が異なるだけであり、それらの数式を、単一配列の場合には、図2に示すハードウェア構成と図8に示すフローチャートに適用し、複数配列の場合には、図5に示すハードウェア構成と図9に示すフローチャートに適用するものであるため、ここでは、詳細な説明は省略する。なお、繰返し更新される未知パラメータの集合の組み合わせは数161とし、各更新式は、上述したように単一配列の場合は、数130、数128、数199〜203であり、複数配列の場合は、数142〜数145、数204〜数210である。
なお、本発明は上記実施例に限定されるものではなく、本発明の趣旨を逸脱しない範囲で変更可能である。例えば、上記各実施例の近似計算において、第3実施例で示した変形例のように遡る過去値を増やしたり、高次の近似にしたりすることも可能である。
以上のように本実施例(単一配列の場合もほとんど同様であるため、ここでは省略する)においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ101と、観測データがどのような確率モデルであるのかを、HMMの未知パラメータを算出することで推定する推定手段114とを備えたHMMの推定装置100において、推定手段114は、HMMの高速化パラメータの値βを設定する初期設定手段122と、前記HMMの未知パラメータとして、状態遷移確率a,分岐確率c,平均値ベクトルμ,共分散行列Σ,初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび分岐の期待値Ncを更新設定する更新設定手段124と、更新設定手段124で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記レジスタ101から読み出した観測データと、初期設定手段122で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する(数142〜数145,数204〜数210を参照)演算手段126と、演算手段126による計算の収束を判定し、計算が収束していなければ、演算手段126で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段126で計算した新たな各確率量を最終的な値として出力させる判定手段128とを備えている。
このようにすれば、HMMの未知パラメータとして、状態遷移確率a,分岐確率c,平均値ベクトルμ,共分散行列Σ,初期状態確率πおよび尤度P(y|θ)の各確率量と、状態遷移の期待値Naおよび分岐の期待値Ncを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なHMM推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
また、ここでの演算手段126は、l+1回目の状態iとなる初期状態確率πi|θl+1(l+1はθの添え字)の値を、前記数142の式で計算し、l+1回目の状態iから状態jに移る状態遷移確率aij|θl+1(t+1はθの添え字)の値を、l回目およびl-1回目における状態遷移の期待値Naij|θlおよびNaij|θl-1(ijはaの添え字、lまたはl-1はθの添え字)を利用して、前記数143〜数145の式で計算し、状態jでのグループ化確率cjg|θl+1(l+1はθの添え字)と、そのグループ化を反映した出力確率bjg|θl+1の値と、平均値ベクトルμjg|θl+1(l+1はθの添え字)の値と、共分散行列Σjg|θl+1(l+1はθの添え字)の値とを、l回目およびl-1回目における前記分岐の期待値Ncjg|θlおよびNcjg|θl-1(jgはcの添え字、lまたはl-1はθの添え字)を利用して、前記数203〜数209の式で計算する構成となっている。
そのため、特に高速化パラメータの値βが1でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
さらに、ここでの更新設定手段124は、高速化パラメータの値を1<β<3に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが3未満の範囲で、βを1よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
なお上述した本実施例の作用効果は、初期設定手段122としての動作を実行する初期設定ステップと、更新設定手段124としての動作を実行する更新設定ステップと、演算手段126としての動作を実行する演算ステップと、判定手段128としての動作を実行する判定ステップとを備えたHMMの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段114として、コンピュータに機能させるHMMの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。