JP5709179B2 - 隠れマルコフモデルの推定方法，推定装置および推定プログラム - Google Patents

Description

本発明は、あるデータ系列が与えられたときに、それを生成する隠れマルコフモデル（ＨＭＭ：Hidden Markov Model）を高速で推定するための隠れマルコフモデルの推定方法，推定装置および推定プログラムに関する。

与えられたデータ系列から隠れマルコフモデルを推定するアルゴリズム（ＨＭＭ再推定アルゴリズム）は、1970年前後にほぼ同時に提案したエル．バウム（L.Baum）氏とエル．ウェルチ（L.Welch）氏の名前をとって、Baum-Welchアルゴリズムと呼ばれている。後に、これはマルコフ性の計算方式に特定したＥＭ（期待値最大化：Expectation-Maximization）アルゴリズムになっていることが分かっている。ＥＭアルゴリズムとは、観測データが現れる尤度（不完全データ尤度：incomplete-data likelihood）が最大となるような確率モデルを得るための手法であり、仮の完全データモデルに基づく対数尤度（log-likelihood）に対して、観測データに基づく条件付期待値を求めるステップ（E-step）と、その期待値の最大化を達成するようにモデルの更新を行なうステップ（M-step）とからなる二つの手順を、交互に繰り返すことで達成される。

一方、ＥＭアルゴリズムは反復的であり、対数よりも有能な代理関数の使用が、汎用的で高速なアルゴリズム構造を提供すると期待される。そこで本願発明者の一人は、ＥＭアルゴリズムを部分集合として含むalpha-ＥＭアルゴリズムを、非特許文献１において提案している。これは、log（対数）関数の拡張版であるalpha-logを用いるものであり、ＨＭＭ再推定アルゴリズムについても、alpha-ＥＭアルゴリズム版（alpha-ＨＭＭ再推定アルゴリズム）があるものと予想されていた。

松山 泰男（Ｙ．Matsuyama），alpha-ＥＭアルゴリズム：alpha-log情報測定を利用した代理の尤度最大化法（The alpha-EM algorithm: Surrogate likelihood maximization using alpha-logarithmic information measures），IEEE情報理論会議（IEEE Trans. on Inform. Theory），第49巻（vol.49），pp.692-706，2003年

しかし、上述したalpha-ＥＭアルゴリズムを用いた通常の導出法では、繰り返しの更新計算時に自分自身の値が必要となって計算できない、言い換えると更新計算する際に、その更新計算の結果と同時系列の項を使用して計算しなければならない、という自己撞着形となってしまい、特殊な条件下以外では、最終的に求めようとする確率モデルの計算ができない。そのため従来は、alpha-logを代理関数としたalpha-ＨＭＭ再推定アルゴリズムは不可能であるとみなされていた。

そこで本発明は、上記問題点を解決して、特殊な条件下以外であっても未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることができる隠れマルコフモデルの推定方法，推定装置および推定プログラムを提供することを目的とする。

本発明は、観測データを時系列に格納する記憶手段と、前記観測データがどのような確率モデルであるのかを、隠れマルコフモデルの未知パラメータを算出することで推定する推定手段とを備えた隠れマルコフモデルの推定装置において、前記推定手段は、前記隠れマルコフモデルの高速化パラメータを設定する初期設定手段と、前記隠れマルコフモデルの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を更新設定する更新設定手段と、前記更新設定手段で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定手段で設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算手段と、前記演算手段による計算の収束を判定し、計算が収束していなければ、前記演算手段で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、前記演算手段で計算した新たな各確率量を最終的な値として出力させる判定手段とを備えて構成される。

この場合の演算手段は、前記新たな各確率量として、ｔ+1回目の状態ｉとなる初期状態確率π_ｉ｜θt+1（t+1はθの添え字）の値を、後述する数６３の式で計算し、ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎ_ａij｜_θtおよびＮ_ａij｜_θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、数６４の式で計算し、状態ｊで状態ｋが出力される出力確率ｂ_jk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎ_ｂjk｜_θtおよびＮ_ｂjk｜_θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、数６５の式で計算する構成とする。

さらに更新設定手段は、前記高速化パラメータの値を１＜β＜３に設定する構成とする。

本発明は、記憶手段に観測データを時系列に格納し、前記観測データがどのような確率モデルであるのかを、推定手段が隠れマルコフモデルの未知パラメータを算出することで推定する隠れマルコフモデルの推定方法において、前記未知パラメータの算出は、前記隠れマルコフモデルの高速化パラメータを設定する初期設定ステップと、前記隠れマルコフモデルの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を更新設定する更新設定ステップと、前記更新設定ステップで更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定ステップで設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算ステップと、前記演算ステップによる計算の収束を判定し、計算が収束していなければ、前記演算ステップで計算した新たな各確率量および各期待値を前記更新設定ステップで設定更新させ、計算が収束していれば、前記演算ステップで計算した新たな各確率量を最終的な値として出力させる判定ステップとにより行なわれる。

この場合の演算ステップは、前記新たな各確率量として、ｔ+1回目の状態ｉとなる初期状態確率π_ｉ｜θt+1（t+1はθの添え字）の値を計算し、ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎ_ａij｜_θtおよびＮ_ａij｜_θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して計算し、状態ｊで状態ｋが出力される出力確率ｂ_jk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎ_ｂjk｜_θtおよびＮ_ｂjk｜_θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して計算する。

また、観測データが離散系列でない場合には、演算ステップは、アルゴリズムの繰り返し指標（インデックス）を示すｌ(エル）を付与した演算式により、状態遷移や分岐の期待値を利用して、さらに出力確率の代わりに分岐確率と、平均値ベクトルと、共分散行列とを計算する。

また、観測データが複数本ある場合には、演算ステップは、その観測データが何本目の配列かを示すインデックスを付して、同様に計算する。

さらに更新設定ステップは、前記高速化パラメータの値を１＜β＜３に設定する。

本発明は、記憶手段に時系列に格納された観測データがどのような確率モデルであるのかを、隠れマルコフモデルの未知パラメータを算出することで推定する推定手段として、コンピュータを機能させる隠れマルコフモデルの推定プログラムにおいて、前記推定手段を、前記隠れマルコフモデルの高速化パラメータを設定する初期設定手段と、前記隠れマルコフモデルの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を更新設定する更新設定手段と、前記更新設定手段で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定手段で設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算手段と、前記演算手段による計算の収束を判定し、計算が収束していなければ、前記演算手段による新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、前記演算手段による新たな各確率量を最終的な値として確定させる判定手段として機能させている。

この場合の演算手段は、前記新たな各確率量として、ｔ+1回目の状態ｉとなる初期状態確率π_ｉ｜θt+1（t+1はθの添え字）の値を計算し、ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎ_ａij｜_θtおよびＮ_ａij｜_θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して計算し、状態ｊで状態ｋが出力される出力確率ｂ_jk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎ_ｂjk｜_θtおよびＮ_ｂjk｜_θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して計算する。

また、観測データが離散系列でない場合には、演算手段は、アルゴリズムの繰り返し指標（インデックス）を示すｌを付与した演算式により、状態遷移や分岐の期待値を利用して、さらに出力確率の代わりに分岐確率と、平均値ベクトルと、共分散行列とを計算する。

また、観測データが複数本ある場合には、演算手段は、その観測データが何本目の配列かを示すインデックスを付して、同様に計算する。

さらに更新設定手段は、前記高速化パラメータの値を１＜β＜３に設定している。

請求項１，４，１４の発明によれば、ＨＭＭの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることができる。

請求項２，５〜１２，１５〜２２の発明によれば、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。

請求項３，１３，２３の発明によれば、繰り返しの計算が発散しない範囲で、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。

本発明の第１実施例を示す推定装置の処理手順を示すフローチャートである。 同上、推定装置の構成を示すブロック図である。 同上、alpha-ＨＭＭアルゴリズムの収束速度について、繰り返し回数と尤度との相関関係を示すグラフである。 本発明の第２実施例を示す推定装置の処理手順を示すフローチャートである。 同上、推定装置の構成を示すブロック図である。 離散型アルファベットの場合のモデル図である。 連続型アルファベットの場合のモデル図である。 本発明の第３実施例、第５実施例、及び第７実施例の単一配列のデータの場合に対応する推定装置の処理手順を示すフローチャートである。 本発明の第４実施例、第６実施例、及び第７実施例の複数配列のデータの場合に対応する推定装置の処理手順を示すフローチャートである。

１，１０１ レジスタ（記憶手段）
１４，１１４ 推定手段
２２，１２２ 初期設定手段
２４，１２４ 更新設定手段
２６，１２６ 演算手段
２８，１２８ 判定手段

以下、添付図面を参照しながら、本発明における隠れマルコフモデルの推定方法，推定装置および推定プログラムの好ましい各実施例について説明する。

先ず、本発明の具体的な実施例を提示する前に、当該実施例における独自のalpha-ＨＭＭ再推定アルゴリズムを導入するまでの理論的な経緯について、以下の数式を参照しながら説明する。

一般に、隠れマルコフモデルによる学習アルゴリズムは、凸ダイバージェンスの最適化から始まって生成されている。まず、次の式に示すように、２つの確率密度p，qの間の凸ダイバージェンスについて考える。

ここで、ＹはＫ次元ユークリッド空間である。関数ｆ(r）は、(０，∞）において凸である。関数f（r）の双対関数g（r）は、次式を満たす必要がある。

上記数１における「≧」という記号は、ｐ＝ｑを殆どどこでも保持する場合に限り保持する。ｆ（1）の正規化が任意であるため、ここでは、Ｄ_ｆ（ｐ・ｑ）とＤ_ｇ（ｑ・ｐ）がｐとｑとの間の擬似距離として作用するように、ｆ（1）を次のように選定する。

留意すべきなのは、次の２つの式が、カルバック・ライブラー・ダイバージェンス（Kullback-Leibler divergence）を発生させる、ということである。もしｐとｑが確率量の関数であれば、数１における積分は総和に置き換えられる。これは、以後説明する積分の全てにおいて同様である。

上記数１〜数５の関係は、対数より有能であるものの、対数と類似した関数であることを意味している。連続して２回微分可能な関数のクラスについて考えると、数２は以下の等式を与える。

また、ｒ＝１付近では以下の式がそれぞれ成り立つ。

ここで、o（1）は高次無限小の項を意味する。もし、乗算的に分離可能な凸関数を使用するのであれば、o（1）は不要であることに注目されたい。

以下の関数は、その一例である。

もし、上記の関数を、ｆ（1）＝０およびｆ”（1）＝１であるようにシフトして正規化すれば、ｋ＝−ｃ（１−ｃ）となる。前記数７のコアとなる項は以下の式を有している。

上式において、次の式は興味深い関数であり、これはパラメータｃによって凸性が調整できる単調関数である。

これは、以下の式の関係が保たれることを指し示すのに重要である。

したがって、Ｌ^(c)（r）はパラメータｃを有する法則化された対数と見なされる。つまり、ｃ-対数と呼ぶことができる。しかしながら、パラメータｃの増加はその数値を減少させる。そこで、ここでは以下の関係によって、法則化された対数をパラメータ化する。

このように、α-対数（alpha-log）は以下の式に示すような凸ダイバージェンスＤ_ｆ（ｐ・ｑ）から始まることによって得られる。

ここでの対数は次式のようになる。

注目すべきは、次の式に示す凸関数が、ｒ=ｐ/ｑとして上記数１に適用されるとαダイバージェンスが発生し、これはalpha-ＥＭアルゴリズムにおいて重要な役割を果たす。

αダイバージェンスは、次の式であらわせる。

上式において、ｐ＝ｑが殆どどこでも適用される場合に限り、等式のゼロが達成される。

次に、alpha-ＨＭＭ再推定アルゴリズムの原点となるalpha-ＥＭアルゴリズムの期待値と最大化の処理ステップについて説明する。ここではＰ_Ｙ｜ψ（ｙ｜ψ）を、ψによってパラメータ化されたＹにおける観測データｙの確率密度または確率量とする。集合Ｙは、不完全データの集合と見なされる。ｘ∈Ｘ（太字）は、未知または消失したデータを含む理想的観測の完全データあるいは増補データとする。そして、不完全データのpdf（probability density function：確率密度関数）あるいはpmf（probability mass function：確率量関数）は次式に示される。

ここでの積分範囲は、次式のようになり、pdfの積分はpmfの総和となる。

そして、条件付きのpdfまたはpmfは次式の通りとなる。

pdfまたはpmfの不完全なデータのalpha-log尤度比は、次式の通りである。

ここで、次の式におけるαダイバージェンスの計算は、alpha-ＥＭアルゴリズムの基本的な関係を示している。

その式は、以下の通りである。

したがって、最大化のための目的関数は以下の式のような定量化（quantity）によって表される。

上記説明によって、alpha-ＥＭアルゴリズムがalpha-ＨＭＭ再推定アルゴリズムの開始点と位置付けられることがわかる。

alpha-ＥＭアルゴリズムでは、次のステップで処理が行なわれる。

初期設定：最初のサイクル用のパラメータφを設定する。

Ｅ-ステップ：φが与えられたときのψの関数として、前記数２５を計算する。

Ｍ-ステップ：次の式にしたがって、更新パラメータを計算する。

Ｕ-ステップ：φをψ^*に置換して、収束するまではＥ-ステップに戻る。
もし、閉じた部分で正確な最大化が得られなければ、数２６に示す最大化が改善した計算に置き換わる。

また、別なalpha-ＧＥＭアルゴリズムでは、次のステップで処理が行なわれる。

初期設定：最初のサイクル用のパラメータφを設定する。

Ｅ-ステップ：φが与えられたときのψの関数として、前記数２５を計算する。

Ｍ-ステップ：次の式にしたがって、更新パラメータを計算する。

Ｕ-ステップ：φをψ^＋に置換して、収束するまではＥ-ステップに戻る。alpha-ＨＭＭの近似されたバージョンが、ＧＥＭアルゴリズムと見なすことができる点に注目されたい。

隠れマルコフモデル（ＨＭＭ）の公式化の問題は、一連のランダムな観測変数データを発生させる最良のマルコフモデルを推定することにある。観測データの変数列は、次の式であらわせる。この式におけるtはマルコフモデルを生成するためのデータの順序を表すものであり、繰り返しのインデックスとは意味が異なる。

ここでの「最良」とは、隠れマルコフモデルの尤度が最大化されることを意味する。隠れマルコフモデルを公式化するために、次に定義する各確率を推定する必要がある。

（ａ）状態遷移確率
状態がｉからｊに遷移する確率のことであり、次式のように定義される。

上記式において、ｓ_i∈Σは有限な範囲での総和Σの中の１つの状態を示すアルファベットの記号である。状態ｓ_kから状態ｓ_lへのつながりが無い場合、状態遷移確率ａ_kl＝０となる。これは、状態遷移における従前のトポロジーを反映している。

（ｂ）出力確率
状態ｓ_jにおいて出力ｙ_kが発生する確率のことであり、次式のように定義される。

ここで、上記式における各要素は以下の通りである。

（ｃ）初期状態確率
最初の状態がｉである確率のことであり、次式のように定義される。

推定される上記（ａ）〜（ｃ）の各要素は、以下のように示される。

つまり、ここで用いる記号θは、状態遷移確率，出力確率および初期状態確率を含ませたものである。

一方、隠れマルコフモデル推定の問題に関し、前記数２８で示したランダムな観測データの変数列は、その値が次の式であらわされる。ここで、tはデータの順序を表すインデックスである。

与えられた変数列は、その変数列を発生できる最大尤度モデルθを推定する。この隠れマルコフモデル推定の問題は、以下の解釈をもたらす。
（ａ）数２８で示す観測データが、前述した不完全データＹとして与えられる。
（ｂ）次の式に示す状態遷移は、消失データＺである。

（ｃ）完全データは、次式のように示される。

この確率は、次式で示される。

（ｄ）前述の三重項θは、数３６に示す完全データの確率構造を提供する。

この解釈は期待値最大化によって、実際にはalpha-ＥＭアルゴリズムによって、前記ＨＭＭ問題の公式化を可能にする。そして、以下の重要な定理が得られる。

定理：ψとφを、２つのマルコフモデルの各パラメータセットとする。そして、α≦１に対して、次の数３８の式は、数３９という意味を持つ。

ここで数３９の不等式は、数３８の不等式が成立する場合に成立するものである。

証明：基本の方程式である数２４において、次式に示す不完全データのalpha-log尤度比を設定する。

数３５に示す状態Ｓ（太字）が消失データであるとすると、Ｘ＝（Ｓ，Ｙ）（各記号は何れも太字）であるので、前記数２５のＱ-関数は次式で算出できる。

前記数１８の特性と数２４の基本的な方程式によって、この定理が保たれる。

数４１の特性は、数４０の形態で推定されるマルコフモデルψが、マルコフモデルφよりも高い確率で数列ｙ（太字）を生成することを意味する。すなわち、 数３３に示すマルコフモデルクラスでのalpha-ＥＭアルゴリズムの繰り返しは、必然的な局所最適性のコストで最良の隠れマルコフモデルを生成する。したがって、数３３に示すマルコフモデルのパラメータで数３８を繰り返すことが、ＨＭＭアルゴリズムの理論版となる。ここでは、α＝１の場合（後述する高速化パラメータの値が、β＝３である場合)は、収束限界であるために実用的ではない。

alpha-ＥＭアルゴリズムとしての主な問題は、ソフトウェアとして実現できる具体的なコンピュータアルゴリズムをいかに提供するのかにある。数３８のalpha-ＨＭＭは理論上の形態であり、これはlog-ＥＭアルゴリズムと同様に一般的なＨＭＭである。ソフトウェアで実行可能なalpha-ＨＭＭとして、ここでは２つのバージョンを提示するが、第１のバージョンは正確であるものの、未来の情報の計算を必要とする。すなわち、得られたアルゴリズムはnon-causalである。一方、本発明で特に提案する第２のバージョンは、時間シフトと確率の近似を用いることによって、未来の情報を含めることなく計算が可能となる。なお、ここで提示するalpha-ＨＭＭアルゴリズムの生成において、全てはpmfであるので、積分が総和となる。以下、それぞれのバージョンについて説明する。

（Ａ）第１の正確なバージョン
・non-causalな更新
正確なバージョンは、条件付き微分によって得られる。このバージョンで得られる更新がnon-causalであることは、既に述べている。この形態は次のバージョンのcausalな形態を得るための原点となるものである。

次式に示すように、全ての更新式は、ｘ＝（ｓ，ｙ）（各記号は何れも太字）に対する前記Ｑ-関数の最大化から得られる。

ここで、ｔはalpha-ＨＭＭに対する繰り返しの指標（インデックス）である。そして、パラメータ集合は次式に示すようになる。なお、明細書中、小文字のπ、a、b等の記号は集合として表記する場合は大文字で示すようにしている。

・Π_θt+1（Πは太字、t+1はθの添え字）の推定
初期状態確率π_ｉ｜θt+1（πは太字、t+1はθの添え字）の更新は、次の修正されたＱ-関数における微分の最大化によって得られる。

この計算は、次式の微分を用いることによって達成できる。

そして、次式を用いることで、数４４のλを除去する。

すると、次の更新式を与える。

上記数４７は観測データｙ（太字）に対して計算できるように見えるが、左辺に含まれるθ_t+1が右辺に含まれている。これは、non-causalityの存在による自己撞着性を示している。しかしながら、これは次のバージョンの近似によって解決されることとなる。したがって、確率のための２つの更なる更新式を得ることができた。

・Ａ_θt+1（Ａは太字、t+1はθの添え字）の推定
状態遷移確率ａ_ij｜θt+1（t+1はθの添え字）の更新は、同様に微分によって得られる。

ここで次式を用いることで、λを除去する。

すると、次の更新式を与える。

ここでのＮ_ij（ｓ）（ｓは太字）は、状態ｓにおいてｉからｊに移行する状態遷移の数である。数５０も右辺にθ_t+1を含んでいる。これも次のバージョンで解決されることとなる。

・Ｂ_θt+1（Ｂは太字、t+1はθの添え字）の推定
出力確率ｂ_jk｜θt+1（t+1はθの添え字）の更新は、再度の微分によって得られる。

ここで次式を用いることで、λを除去する。

すると、次の更新式を与える。

ここでのＮ_ｂjk（ｓ）（ｓは太字、jkはｂの添え字）は、数列ｓ（太字）によって生じる出力確率ｂ_jkに対する事象の数である。数５３も右辺にθ_t+1を含んでいる。これは次のバージョンで解決されることとなる。

（Ｂ）第２のCausalな近似バージョン
前述したように、第１のバージョンに示す一連の更新式（数４７，数５０，数５３）は、自己撞着性を解決するために近似を必要とする。

・時間のシフト
alpha-ＥＭアルゴリズムの収束は、適切な収束判定基準によって、次の式のようになることを意味する。

これは、以下の近似式を与える。ここで、o(1)は高位の無限小を意味する。

これは、繰り返し指標tをシフトすることによって、数４７，数５０，数５３に示す更新式の自己撞着性を解決することになる。

・αのシフト
前記時間のシフトは、尤度比の期待値のための確率環境が反転されることを意味する。したがって、元のパラメータαはα_causalに変換される。αとα_causalの関係は、以下の誘導によって得られる。ここで、数４７，数５０，数５３における尤度比の期待値の核となる部分を考慮すると、数４７の場合は以下の関係が成り立つ。

この誘導は尤度比の時系列を維持しながら、期待値が現在の環境Ｐ（ｓ｜ｙ，θ_t）（ｓ，ｙは太字）から、未来の環境Ｐ（ｓ｜ｙ，θ_t+1）（ｓ，ｙは太字）に変化することを意味する。上記数５６を見るとθ_ｔ＋１の初期確率π_iを左辺として計算するのに、第一式の右辺では同じ時系列のθ_ｔ＋１に依存する関数（項）をしなければならないという問題点があり、このままでは計算できないと言う自己撞着性を含んでいた。そこで、まず第二式では第一式の分母分子を｛｝内の項のα乗でそれぞれ除算する。さらに第三式ではテーラー展開によって微小近似すると共に時間シフトの概念を導入しており、これによって自己撞着矛盾を解決している。言い換えると、右辺には左辺と同時系列のθ_ｔ＋１に依存する関数（項）を用いない形になっている。この詳細な考え方を以下に示す。

ａ）未来に関する尤度比を現在の環境で計算すると、次式のようになる。

ｂ）未来に関する尤度比を未来の環境で計算すると、次式のようになる。

したがって、時間シフト前後の関係は以下の通りである。

そこから得られる関係式は、次式のようになる。

したがって、α＝−１のlog-ＥＭとlog-ＨＭＭは、ここでのα_causal＝１の場合に相当する。なお、符号の使用を簡単化するために、以後、高速化パラメータの値となるα_causalを、次式のようにβとしてあらわす。すると、数６０の右側の条件式はβ≦３を示すことになり、従来技術であるlog-ＥＭとlog-ＨＭＭはβ＝１の場合に対応することがわかる。（ただし、β＝３だと後述の実験結果より発散してしまうことがわかっている。）

・テーラー展開
前記全ての更新式（数４７，数５０，数５３）は、尤度比の計算に能力を必要とする。これは計算に時間がかかるものであり、Ｔが増加するに従って扱いにくくなる。テーラー展開は、alpha-logアルゴリズムにおける尤度比の長所を失うことなく、計算を簡単化することができるもので、尤度比は次式のように近似化される。

数５５と数６２を数４７，数５０，数５３に適用すると、過去の情報を利用するcausalで、かつ計算上効率的なalpha-ＨＭＭアルゴリズムが得られる。注目すべきは、これが数２７のalpha-ＧＥＭアルゴリズムに対応するということである。

（Ｃ）ソフトウェアで実行可能なalpha-ＨＭＭアルゴリズム
・Π_θt+1（Πは太字、t+1はθの添え字）の推定
高次の項ｏ（1）を放棄することによって数５５と数６２の近似式を結合し、数４７に適用すると以下の更新式が与えられる。

・Ａ_θt+1（Ａは太字、t+1はθの添え字）の推定
数５５と数６２の近似式を適用すると、数５０は以下の更新式を与える。

ここで、Ｎ_ａij｜θt（ijはａの添え字、tはθの添え字）は、条件付き確率Ｐ（ｓ｜ｙ，θ_t）（ｓ，ｙは太字）の下でのＮ_ij（ｓ）（ｓは太字）の期待値であり、不完全データｙ（ｙは太字）によって生じる状態遷移のカウント数により推定される。

・Ｂ_θt+1（Ｂは太字、t+1はθの添え字）の推定
数５５と数６２の近似式を適用すると、数５３は以下の更新式を与える。

ここで、Ｎ_ｂjK｜θt（jkはｂの添え字、tはθの添え字）は、条件付き確率Ｐ（ｓ｜ｙ，θ_t）（ｓ，ｙは太字）の下でのＮ_ｂjK（ｓ）（ｓは太字）の期待値であり、不完全データｙ（ｙは太字）によって生じる状態遷移のカウント数によって推定される。

・従来のlog-ＨＭＭとの比較
上記の数６３，数６４，数６５は、ソフトウェアとして実行可能である。これらの更新式は、以下に示す共通の特性を有する。
（１）β＝１すなわちα＝−１の場合は、log-ＥＭアルゴリズムから得られた従来のＨＭＭアルゴリズムと一致する。
（２）数６０と数６１との間の関係により、αが行なうのと同様に、パラメータβは更新用の現在と過去の情報の各要素を調整する。
（３）βに依存する、すなわちαに依存するすべての過去の更新値は、格納データを参照することによってのみ得ることができる。
（４）現在の項は、log-ＨＭＭアルゴリズムに対する前方−後方アルゴリズムを用いることで効率的に計算できる。

上記特性（１）〜（４）は、alpha-ＨＭＭアルゴリズムの各更新サイクルがlog-ＨＭＭアルゴリズムのオーバーヘッドより僅かに増えてだけであることを示唆している。したがって、設計パラメータα、あるいは高速化パラメータβが適切に選択されれば、本発明のalpha-ＨＭＭアルゴリズムは、従来のlog-ＨＭＭアルゴリズムよりも高速に収束することが予測される。この違いが生じる理由は、数６３、数６４、数６５の式において、従来のlog-ＥＭから得られたＨＭＭアルゴリズムに相当するβ＝１を代入すると分母分子の第二項が消去され、分母分子の第一項が約分されて尤度比が残らなかったのに対し、β＞1で設定すると分母分子の第二項が消去されず、単純に分母分子を約分できないために尤度比の影響を残すことができるようになったことだと考えられる。

以下、上記理論に基づく好ましい実施例について説明する。

実施例１、２では、観測データが離散的な場合において、観測データが単一配列（一本の観測データ）を有する場合、複数配列（Ｍ本の観測データ）を有する場合について、それぞれ説明する。

実施例３、４では、観測データが連続的な場合において、観測データが単一配列（一本の観測データ）を有する場合、複数配列（Ｍ本の観測データ）を有する場合について、それぞれ説明する。

実施例５、６では、観測データが半連続的な場合において、観測データが単一配列（一本の観測データ）を有する場合、複数配列（Ｍ本の観測データ）を有する場合について、それぞれ説明する。

実施例７では、観測データが離散的かつ連続的な場合、すなわち離散系列と連続系列が混在している場合において、単一配列を有する場合と複数配列を有する場合をまとめて説明する。

図１は、本発明の第１実施例において、上述したalpha-ＨＭＭ再推定アルゴリズムを実行可能にするプログラムの処理手順をあらわしたものである。

同図において、１は合計でＴ個の観測データを格納する記憶手段としてのレジスタで、各観測データは時刻τが１から順に時系列に並んで格納される。ここでｙ_τは個々のデータ値を示し、ｙはＴ個のデータ値の集合を示しており、本実施例では配列が１本の観測データをレジスタ１に格納している。推定装置１０は、レジスタ１に記憶される一列のデータ値が、どのような確率構造（モデル）を有しているのかを、ＨＭＭの未知パラメータを算出する以下のステップＳ１〜ステップＳ１０の各手順に従って推定解析するものである。

ステップＳ１は、前述した高速化パラメータの値βを設定する部分である。従来のＨＭＭ推定アルゴリズムは、β＝１の場合に相当する。つまり、β＝１という特殊な条件下であれば、従来のalpha-ＨＭＭ再推定アルゴリズムであっても、上記確率モデルの計算は可能である。しかし、本実施例で提案する新たなalpha-ＨＭＭ再推定アルゴリズムでは、確率モデルの計算が可能な高速化パラメータの値βを１≦β＜３の範囲に拡張することができる。これは、従来のβ＝１を特例として含むものである。なお、βが大きな値であるほど、推定装置１０としての処理の高速性は増すが、収束性を保持するにはβ＜３でなければならない。

ステップＳ２は、初期確率と収束判定値を決める部分である。これは、後述するステップＳ３〜ステップＳ８の手順を繰り返す前に行なわれる。推定装置１０が最終的に算出しようとするＨＭＭの確率構造は、次式のようにあらわせる。これは、前述の数３７と等しい。

上記式において、π_S0（０はｓの添え字）は最初の状態ｓ_０における確率であり、ａ_Sτ-1Sτ（τ−１およびτはｓの添え字）は時刻τ−１の状態ｓ_τ-1が時刻τの状態ｓ_τに移る確率であり、ｂ_Sτ（ｙ_τ）（τはｓの添え字）は時刻τの状態ｓ_τに移行したときにｙ_τが出力される確率である。推定装置１０は、レジスタ１に格納された観測データを読み出して、その観測データが最も出現しやすくなる初期状態確率πと、状態遷移確率ａと、出力確率ｂとを推定するが、このステップＳ２では、次式に示すように、それらの確率π，ａ，ｂの初期値を決定する。

上式において、θ_０は数４３にあるように、０回目の繰り返しにおける初期の確率π，ａ，ｂの組合せを示し、推定装置１０は、その条件で状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θ０（０はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂ_jk｜θ0（０はθの添え字）の値と、最初に状態ｉとなる初期状態確率π_ｉ｜θ0（０はθの添え字）の値とをそれぞれ決定する。

またステップＳ２では、対数尤度に基づく収束範囲を決めることで、ステップＳ３〜ステップＳ９の繰り返しを終了させるための収束判定値を決定する。この収束判定値はステップＳ８で用いられ、具体的には後述する数８３で示される。

これらの値をステップＳ１，Ｓ２で決定すると、ステップＳ３の手順に移行して、確率量とカウント値の設定が行なわれる。ステップＳ３は、最初にステップＳ２で求めた初期確率を利用し、それ以降はステップＳ９で更新された確率量とカウント値を利用して、数６８に示す各確率値と、数６９に示す各カウント値を、実際に推定装置１０のメモリ（図示せず）に設定する部分である。

次のステップＳ４は、レジスタ１から観測データを読み出して、ステップＳ３で設定された上式の確率量から、数７０と数７１に示す前向き確率を計算する部分である。なお、ここに示すαは確率値であり、前述したパラメータの値αとは異なる。ここでは、数７０に示す確率値αが計算され、そこから数７１に示す確率値（尤度）Ｐ（ｙ｜θ_t）が計算される。

次のステップＳ５は、レジスタ１から観測データを読み出して、ステップＳ３で設定された上式の確率量から、数７２に示す前向き確率を計算する部分である。なお、ここに示すβは確率値であり、前述したパラメータの値βとは異なる。ここでは、数７２に示す確率値βが計算される。

上記前向き確率と後向き確率は、計算回数を減らすために、既存のＨＭＭ推定アルゴリズムにも組み込まれていたものである。

続くステップＳ６は、ステップＳ４で計算された確率値αと、ステップＳ５で計算された確率値βを用い、レジスタ１から観測データを読み出して、数７３に示す状態遷移のカウント値と、数７４に示す出力のカウント値をそれぞれ計算する部分である。

次のステップＳ７は、t+1回目に繰り返される確率π，ａ，ｂの組合せを条件として、状態ｉとなる初期状態確率π_ｉ｜θt+1（t+1はθの添え字）の値と、状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θt+1（t+1はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂ_jk｜θt+1（t+1はθの添え字）の値とを、ステップＳ４で求めたt回目の確率Ｐ（ｙ｜θ_t）と、それよりも１回前の確率Ｐ（ｙ｜θ_t-1）を利用し、且つステップＳ５で求めたt回目のカウント値Ｎ_ａij｜_θt，Ｎ_ｂjk｜_θt（ijはａの添え字、jkはｂの添え字、tはθの添え字）と、それよりも１回前のカウント値Ｎ_ａij｜_θt-1，Ｎ_ｂjk｜_θt-1（ijはａの添え字、jkはｂの添え字、t-1はθの添え字）を利用して計算する部分であり、これは前記数６３，数６４，６５に対応している。

このステップＳ７では、前記ステップＳ１で設定したパラメータの値βが用いられ、本アルゴリズムの主要な特徴部分となっている。特に、添え字ｔ+1の部分の確率を計算する上で、時間シフトした添え字tや添え字t-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したことが、その効果として挙げられる。

本実施例では、配列が１本の観測データをレジスタ１に格納しているが、その場合のステップＳ７における確率用の更新式は、次式のようにあらわせる。

ステップＳ８は、計算の収束を判定する部分である。ここでは、直前のステップＳ７で計算された新たな確率量に基づく尤度Ｐを用い、その尤度Ｐと前記ステップＳ２における収束判定値との比較により、ステップＳ３〜ステップＳ９の計算が収束したか否かを判定する。近似的には、繰り返しの回数を指定して、指定回数に達したら、計算が収束したと判定してもよい。

ステップＳ８において計算が収束していないと判定された場合、ステップＳ９に移行して手順が繰り返される。ステップＳ９は、計算された確率量と事象のカウント値を、繰り返し計算のために更新する部分である。本実施例では、次式に示すように、２つの過去値をシフトさせる。

ここでは、繰り返しの回数が１つ増えることにより、t回目の確率量およびカウント値がt-1回目の確率量およびカウント値に更新され、t+1回目の確率量およびカウント値がt回目の確率量およびカウント値に更新される。なお、ステップＳ３とステップＳ９を一纏めにして、更新された確率量およびカウント値をそのまま次の繰り返しのために、推定装置１０のメモリに設定記憶させてもよい。

一方、ステップＳ８において計算が収束していると判定された場合、ステップＳ１０に移行して、次式に示す計算された三組の確率量を用いたＨＭＭを採用する。

上記数７９で計算した各値を用いて、推定装置１０は数６６に示すＨＭＭの確率構造を出力することができる。

図２は、図１に示すアルゴリズムを実現する推定装置１０の構成を示している。同図において、推定装置１０は前記レジスタ１を内蔵する入力手段１２と、レジスタ１に格納した観測データを適宜読み出し、上記ステップＳ１〜ステップＳ１０の手順を実行して数６６に示すＨＭＭの確率構造を推定する推定手段１４と、この推定手段１４で得たＨＭＭの確率構造を出力する出力手段１６とにより概ね構成される。本実施例の入力手段１２は、１本の配列を有する観測データをレジスタ１に格納できる構造となっているが、複数本の配列を有する観測データを取り扱う推定装置１００については、次の実施例で詳しく説明する。

推定装置１０のハードウェア構成は、例えば演算処理部としてＣＰＵを備えたコンピュータで実現することができる。その場合、図１に示す推定装置１０の処理手順を実行するプログラムが、メモリなどの記録媒体に記憶される。当該プログラムをどこに記憶するのかは限定せず、例えば通信手段を介してプログラムがコンピュータにダウンロードされる構成であってもよい。

推定手段１４は、前記ステップＳ１における高速化パラメータの値βや、ステップＳ２における初期の確率値および収束判定の条件を設定する初期設定手段２２と、ステップＳ３における確率量およびカウント値の設定や、ステップＳ９における確率量およびカウント値の更新を行なう更新設定手段２４と、ステップＳ４における前向き確率の計算や、ステップＳ５における後ろ向き確率の計算を行なうと共に、そこからステップＳ６における状態遷移のカウント値と出力のカウント値をそれぞれ計算し、さらにステップＳ７における新たな確率量の計算を行なう演算手段２６と、ステップＳ８における計算の収束判定を行ない、計算が収束していなければ、前記更新設定手段２４による確率量およびカウント値の更新を行なわせる一方で、計算が収束していれば、ステップＳ１０において、演算手段２６が直前に計算した新たな確率量を、ＨＭＭの最終的なパラメータ値として確定させる判定手段２８とを備えている。

また推定装置１０には、前記高速化パラメータの値や、初期の確率値および収束判定の条件や、更新設定される確率量およびカウント値の値を読み書き可能に格納するメモリ３０の他に、必要に応じて高速化パラメータの値や収束判定の条件を操作入力するキーボードやマウスなどの操作手段３２が、推定手段１４に接続して設けられる。これらの各装置構成により、上述したステップＳ１〜ステップＳ１０の手順が実行される。

次に、本実施例における試験用の観測データでの実験結果を説明する。

（Ａ）生成されたデータを用いた速度の評価
観測データの集合が強い凸性でなければ（そして強い凹面でもなければ）、全ての最適化アルゴリズムは局所的に最適となる。存在するデータの殆ど全ては、そのような性質を持っている。従来のあるいはlog-ＨＭＭアルゴリズムは例外的ではないし、alpha-ＨＭＭアルゴリズムも例外的ではない。したがって、そうしたＨＭＭアルゴリズムを適用する前に、発生のメカニズムが分かっている人工的なデータ集合を生成することが必要になってくる。

ここでは、数３３に示す数列をレジスタ１に格納する入力データとし、そのパラメータを次式のように指定する。

このマルコフ連鎖を使用し、[0, 1]上の一様な乱数によって、1000個のサンプルを有するデータ集合を生成した。完全な推定には、以下の困難があることを先に述べる。
（ａ）各データ集合の殆どが、局所的最適を与える対象である。
（ｂ）たとえ局所的最適が避けられたとしても、正確な推定には無限の長さの数列が必要である。
（ｃ）コンピュータが生成した乱数は単なる擬似的な乱数にすぎない。

したがって、以下が満たされれば満足するものとする。
（１）数３３に示すパラメータ集合の数値が、数８０からずれないこと
（２）alpha-ＨＭＭアルゴリズムの収束したパラメータ集合がlog-ＨＭＭのそれと近いこと
ステップＳ２における初期値は、次のように設定する。これにより、初期状態を除いて事前の情報は必要としないものとする。

以下の理由のため初期状態は固定した。
（ａ）分布された初期確率は、不安定な局所最適性を得やすい。これは、alpha-ＨＭＭアルゴリズムの各収束速度を比較するのに不要な曖昧さをもたらす。
（ｂ）ＨＭＭ推定の後の認識問題において、一つの初期状態はビタビ（Viterbi）アルゴリズムによって選択される。

図３は、対数尤度（log-likelihood）、すなわち繰り返し回数を考慮した尤度Ｐ（ｙ｜θ_t）に関する収束の傾向を示している。ここでは、数６０および数６１を満たす多様な値で実験した。

この図から以下のことがわかる。
（１）β=α_causalが増加するにつれて対数尤度の立ち上がりが早くなる。すなわち、収束するまでの繰り返し回数が少なくなる。
（２）alpha-ＥＭアルゴリズムの限界(α＝１)は、alpha-ＨＭＭアルゴリズムではβ=α_causal＝３に相当するが、この場合、発散してしまう。

さらに２つの視点で、alpha-ＨＭＭアルゴリズムの性能を比較する必要がある。一つは局所最適の比較であり、もう一つはＣＰＵ時間である。局所最適の比較では、次式に示すように、θに対する全パラメータと共に、それに付随した対数尤度の値を一覧にする必要がある。

上記の実験結果の一覧を見るとβを1〜2.75の範囲で変化させても対数尤度（ＬＬ）の値がほとんど変化せずに正確さを保っていることがわかる。この実験ではβ≧３に設定した実験は発散する結果になるために省略しているが、３に近い2.75に設定しても実験結果にそれほど影響なく高速化できることがわかる。ただし、このような実験結果の収束限界となるβの値は入力データによって若干異なる。

次は、収束するまでのＣＰＵ時間の比較である。この比較はソフトウェアの実行手段とハードウェアの手法に依存する。ここでは、alpha-ＨＭＭアルゴリズムを標準的なPOSIX環境におけるＣ‐コードとして実行する。これらのコードは標準的なシングルコアプロセッサで走らせている。前記ステップＳ８で行われる収束は、以下の判定基準で測定した。

上式において、Ｐ_newは今回の計算で得た確率Ｐの値であり、Ｐ_oldは前回の計算で得た確率Ｐの値である。この収束判定基準は更新前後で確率の値がほとんど変化していなければ収束したと判断するものである。この収束判定式結果は次の表に示す通りであり、ここでは繰り返し回数とＣＰＵ時間による収束の比較を示している。

局所最適性とＣＰＵ時間とを合せて比較すると、以下のことがわかる。
（ａ）収束する各対数尤度は非常に接近しているものの、何れも局所最大値が存在する。
（ｂ）対数のまたは従来のＨＭＭアルゴリズムは、常に最良の対数尤度をもたらす訳ではない。log-ＨＭＭアルゴリズムよりも良好な局所最大値が、alpha-ＨＭＭアルゴリズムで得られた。
（ｃ）改良したalpha-ＨＭＭアルゴリズムによって、繰り返し回数とＣＰＵ時間の高速化が成し遂げられたのは明らかである。
（ｃ）繰り返し回数の改善が、直接ＣＰＵ時間に影響している。これは前記数６３，数６４および数６５に示すように、予め計算されてメモリに格納される過去の情報を利用しているからであり、従来のlog-ＨＭＭアルゴリズムに比べて、ソフトウェアの複雑さを増やしてはいない。
（ｄ）一方、α_causal＝β＝３すなわちα＝１の場合は、発散する。これは、alpha-ＥＭアルゴリズムの許容範囲と一致する。したがって、数６２による近似は許容される。α_causal＝β＝１すなわちα＝−１の場合において、数６３，数６４および数６５はlog-ＨＭＭアルゴリズムの正確な形態であることが重要である。この意味において、数６３，数６４および数６５の更新式はlog-ＨＭＭアルゴリズムの拡張版ということができる。

前記図３を参照すると、従来のalpha-ＨＭＭアルゴリズムで唯一ステップＳ６からステップＳ７への計算が可能であったβ＝１の結果に対して、今回提案するalpha-ＨＭＭアルゴリズムで計算が可能になったβ＝2.75の結果では、少ない繰り返し回数で推定装置１０の処理が収束していることが判る。具体的には、上記表１において、β＝１の場合は繰り返し回数が263回で、ＣＰＵの使用時間が0.253秒となっているが、β＝2.75の場合は繰り返し回数が70回、またＣＰＵの使用時間が0.068秒に減少する。表１の最右列にあるのは、β＝１のときの繰り返し回数とＣＰＵの使用時間を１としたときの速度比を計算したもので、例えばβ＝2.75に設定すれば繰り返し回数は3.76倍、またＣＰＵの使用時間は3.72倍に高速化する。

以上の実験結果から、本実施例では従来よりも高速なＨＭＭアルゴリズムが示された。本推定方法は、基になっているalpha -ＥＭアルゴリズムを反映しており、改良したalpha-ＨＭＭアルゴリズムと呼ぶことができる。改良したalpha-ＨＭＭアルゴリズムは、従来のalpha-ＨＭＭアルゴリズムまたはlog-ＨＭＭアルゴリズムより優れている。計算上の複雑さの増加が非常に少ないので、繰り返し回数の抑制が直接的に期待したＣＰＵの高速化を実現した。

また、改良したalpha-ＨＭＭアルゴリズムで採用する数６３，数６４および数６５の更新式は、メモリに格納した過去の情報を追加的に必要とするだけである。したがって、既存のコードからのソフトウェアのバージョンアップは難しくない。この場合、alpha-対数の曲率を制御するα_causal＝α＋２が設計パラメータとなる。

さらに、より少ないサンプルで未来の情報を利用する方法について、確認を行なった。この方法は、テーラー展開と共にnon-causalな形態を使用したことにある。更新式の一例は、次式のように表せる。

ここで、θ＾_ｔ＋１とθ＾_ｔ（＾はθの上部に記される）は、より少ないサンプル（例えば1000個のサンプルから200個を取り出す）から統計量が推定されることを示している。すなわち、 サンプルの一部が未来を考えるのに使われる。この方法は数６３，数６４および数６５の更新式よりも大きなオーバーヘッドを必要とするが、高速化の利益の方がこの余計な重荷よりも勝っている。もし、入力源がよく混合されている、あるいはエルゴードであると演繹的知識をユーザーが持っていれば、この方法を利用できる。実験では、高速化が数６３，数６４および数６５の方法よりも悪くなかった。

以上のように本実施例では、入力される観測データを時系列に格納する記憶手段としてのレジスタ１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１４とを備えたＨＭＭの推定装置１０において、推定手段１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，出力確率ｂ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび出力の期待値Ｎ_ｂを更新設定する更新設定手段２４と、更新設定手段２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数６８，数６９を参照）と共に、前記レジスタ１から読み出した観測データと、初期設定手段２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数６３，数６４，数６５，数７５，数７６，数７７を参照）演算手段２６と、演算手段２６による計算の収束を判定し、計算が収束していなければ、演算手段２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段２６で計算した新たな各確率量を最終的な値として出力させる判定手段２８とを備えている。

このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，出力確率ｂ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび出力の期待値Ｎ_ｂを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。

また、ここでの演算手段２６は、ｔ+1回目の状態ｉとなる初期状態確率π_ｉ｜θt+1（t+1はθの添え字）の値を、前記数６３の式で計算し、ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における状態遷移の期待値Ｎ_ａij｜_θtおよびＮ_ａij｜_θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、前記数６４の式で計算し、状態ｊで状態ｋが出力される出力確率ｂ_jk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎ_ｂjk｜_θtおよびＮ_ｂjk｜_θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、前記数６５の式で計算する構成となっている。

そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。

さらに、ここでの更新設定手段２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。

なお上述した本実施例の作用効果は、初期設定手段２２としての動作を実行する初期設定ステップと、更新設定手段２４としての動作を実行する更新設定ステップと、演算手段２６としての動作を実行する演算ステップと、判定手段２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。

図４は、本発明の第２実施例におけるプログラムの処理手順をあらわしたものである。

同図において、１０１は合計でＴ個の観測データを格納する記憶装置としてのレジスタで、各観測データのそれぞれは、時刻τが１から順に時系列に並んで格納される。ここでｙ_τは個々のデータ値を示し、ｙはＴ個のデータ値の集合を示し、ｎは何本目の配列かをあらわすインデックスを示しており、本実施例では複数の配列を有するＭ（＝２以上の整数）本の観測データをレジスタ１０１に格納している。推定装置１００は、レジスタ１０１に記憶されるＭ列のデータ値が、どのような確率構造（モデル）を有しているのかを、以下のステップＳ１１〜ステップＳ２０の各手順に従って推定解析するものである。

ステップＳ１１は、前述した高速化パラメータの値βを設定する部分である。これは上記第１実施例のステップＳ１に相当するもので、高速化パラメータの値βは、従来のβ＝１を特例として含む１≦β＜３の範囲に拡張される。なお、βが大きな値であるほど、推定装置１００としての処理の高速性は増すが、収束性を保持するにはβ＜３でなければならない。

ステップＳ１２は、初期確率と収束判定値を決める部分である。これは、前記ステップＳ２に相当するもので、後述するステップＳ１３〜ステップＳ１８の手順を繰り返す前に行なわれる。推定装置１００は、レジスタ１０１に格納された観測データを読み出して、その観測データが最も出現しやすくなる初期状態確率πと、状態遷移確率ａと、出力確率ｂとを推定するが、このステップＳ１２では、次式に示すように、それらの確率π，ａ，ｂの初期値を決定する。

上式において、θ_０は数４３にあるように、０回目の繰り返しにおける初期の確率π，ａ，ｂの組合せを示し、推定装置１００は、その条件で状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θ０（０はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂ_jk｜θ0（０はθの添え字）の値と、最初に状態ｉとなる初期状態確率π_ｉ｜θ0（０はθの添え字）の値とをそれぞれ決定する。

またステップＳ１２では、対数尤度に基づく収束範囲を決めることで、ステップＳ１３〜ステップＳ１９の繰り返しを終了させるための収束判定値を決定する。この収束判定値はステップＳ１８で用いられ、具体的には前記数８３で示される。

これらの値をステップＳ１１，Ｓ１２で決定すると、ステップＳ１３の手順に移行して、確率量とカウント値の設定が行なわれる。ステップＳ１３は、最初にステップＳ１２で求めた初期確率を利用し、それ以降はステップＳ１９で更新された確率量とカウント値を利用して、数８６に示す各確率値と、数８７に示す各カウント値を、実際に推定装置１００のメモリ（図示せず）に設定する部分である。

次のステップＳ１４は、レジスタ１０１から観測データを読み出して、ステップＳ３で設定された上式の確率量から、数８８と数８９に示す前向き確率を計算する部分である。なお、ここに示すαは確率値であり、前述したパラメータの値αとは異なる。ここでは、数８８に示す確率値αが計算され、そこから数８９に示す確率値（尤度）Ｐ（ｙ⁽ⁿ⁾｜θ_t）が計算される。

次のステップＳ１５は、レジスタ１０１から観測データを読み出して、ステップＳ１３で設定された上式の確率量から、数９０に示す前向き確率を計算する部分である。なお、ここに示すβは確率値であり、前述したパラメータの値βとは異なる。ここでは、数９０に示す確率値βが計算される。

上記前向き確率と後向き確率は、計算回数を減らすために、既存のＨＭＭ推定アルゴリズムにも組み込まれていたものである。

続くステップＳ１６は、ステップＳ１４で計算された確率値αと、ステップＳ５で計算された確率値βを用い、レジスタ１０１から観測データを読み出して、数９１に示す状態遷移のカウント値と、数９２に示す出力のカウント値をそれぞれ計算する部分である。

次のステップＳ１７は、t+1回目に繰り返される確率π，ａ，ｂの組合せを条件として、状態ｉとなる初期状態確率π_ｉ｜θt+1（t+1はθの添え字）の値と、状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θt+1（t+1はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂ_jk｜θt+1（t+1はθの添え字）の値とを、ステップＳ１４で求めたt回目の確率Ｐ（ｙ⁽ⁿ⁾｜θ_t）と、それよりも１回前の確率Ｐ（ｙ⁽ⁿ⁾｜θ_t-1）を利用し、且つステップＳ５で求めたt回目のカウント値Ｎ⁽ⁿ⁾ _ａjk｜_θt，Ｎ⁽ⁿ⁾ _ｂjk｜_θt（jkはａの添え字、tはθの添え字）と、それよりも１回前のカウント値Ｎ⁽ⁿ⁾ _ａjk｜_θt-1，Ｎ⁽ⁿ⁾ _ｂjk｜_θt-1（jkはａの添え字、t-1はθの添え字）を利用して計算する部分であり、これは前記数６３，数６４，６５に対応している。

このステップＳ１７では、前記ステップＳ１１で設定したパラメータの値βが用いられ、本アルゴリズムの主要な特徴部分となっている。特に、添え字ｔ+1の部分の確率を計算する上で、時間シフトした添え字tや添え字t-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したことが、その効果として挙げられる。

本実施例では、配列がＭ本の観測データをレジスタ１０１に格納しているが、その場合のステップＳ１７における確率用の更新式は、次式のようにあらわせる。

ステップＳ１８は、計算の収束を判定する部分である。ここでは、直前のステップＳ１７で計算された新たな確率量に基づく尤度Ｐを用い、その尤度Ｐと前記ステップＳ１２における収束判定値との比較により、ステップＳ１３〜ステップＳ１９の計算が収束したか否かを判定する。近似的には、繰り返しの回数を指定して、指定回数に達したら、計算が収束したと判定してもよい。

ステップＳ１８において計算が収束していないと判定された場合、ステップＳ１９に移行して手順が繰り返される。ステップＳ１９は、計算された確率量と事象のカウント値を、繰り返し計算のために更新する部分である。本実施例では、次式に示すように、２つの過去値をシフトさせる。

ここでは、繰り返しの回数が１つ増えることにより、t回目の確率量およびカウント値がt-1回目の確率量およびカウント値に更新され、t+1回目の確率量およびカウント値がt回目の確率量およびカウント値に更新される。

一方、ステップＳ１８において計算が収束していると判定された場合、ステップＳ２０に移行して、次式に示す計算された三組の確率量を用いたＨＭＭを採用する。

上記数９７で計算した各値を用いて、推定装置１０は数６６に示すＨＭＭの確率構造を出力することができる。

図５は、図４に示すアルゴリズムを実現する推定装置１００の構成を示している。同図において、推定装置１００は前記レジスタ１０１を内蔵する入力手段１１２と、レジスタ１０１に格納した観測データを適宜読み出し、上記ステップＳ１１〜ステップＳ２０の手順を実行して数６６に示すＨＭＭの確率構造を推定する推定手段１１４と、この推定手段１１４で得たＨＭＭの確率構造を出力する出力手段１１６とにより概ね構成される。本実施例の入力手段１１２は、Ｍ本の配列を有する観測データをレジスタ１０１に格納できる構造となっている。推定装置１００のハードウェア構成については、前記第１実施例の推定装置１０と同様であるため、ここでは説明を省略する。

推定手段１１４は、初期設定手段１２２と、更新設定手段１２４と、演算手段１２６と、判定手段１２８とを備えている。これらは前記第１実施例の初期設定手段２２，更新設定手段２４，演算手段２６および判定手段２８にそれぞれ対応するもので、取り扱う観測データがＭ本の配列になった以外は、第１実施例と同様に機能する。

また推定装置１１０には、前記高速化パラメータの値や、初期の確率値および収束判定の条件や、更新設定される確率量およびカウント値の値を読み書き可能に格納するメモリ１３０の他に、必要に応じて高速化パラメータの値や収束判定の条件を操作入力するキーボードやマウスなどの操作手段１３２が、推定手段１１４に接続して設けられる。これらの各装置構成により、上述したステップＳ１１〜ステップＳ２０の手順が実行される。

以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１０１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１１４とを備えたＨＭＭの推定装置１００において、推定手段１１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段１２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，出力確率ｂ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび出力の期待値Ｎ_ｂを更新設定する更新設定手段１２４と、更新設定手段１２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数６８，数６９を参照）と共に、前記レジスタ１０１から読み出した観測データと、初期設定手段１２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数６３，数６４，数６５，数９３，数９４，数９５を参照）演算手段１２６と、演算手段１２６による計算の収束を判定し、計算が収束していなければ、演算手段１２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段１２６で計算した新たな各確率量を最終的な値として出力させる判定手段１２８とを備えている。

このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，出力確率ｂ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび出力の期待値Ｎ_ｂを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。

また、ここでの演算手段１２６は、ｔ+1回目の状態ｉとなる初期状態確率π_ｉ｜θt+1（t+1はθの添え字）の値を、前記数６３の式で計算し、ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における状態遷移の期待値Ｎ_ａij｜_θtおよびＮ_ａij｜_θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、前記数６４の式で計算し、状態ｊで状態ｋが出力される出力確率ｂ_jk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎ_ｂjk｜_θtおよびＮ_ｂjk｜_θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、前記数６５の式で計算する構成となっている。

そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。

さらに、ここでの更新設定手段１２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。

なお上述した本実施例の作用効果は、初期設定手段１２２としての動作を実行する初期設定ステップと、更新設定手段１２４としての動作を実行する更新設定ステップと、演算手段１２６としての動作を実行する演算ステップと、判定手段１２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。

上記の実施例１、２では、観測データが離散的な場合について述べた。本実施例では、観測データが単一配列で連続的な場合の実施形態について説明する。

この観測データが連続的な場合とは、観測データのデータ列の各々の値が波のように変動しているような場合である。

具体的には、図６の観測データが離散的な場合のモデル図における出力確率ｂ_ik，ｂ_jkが、図７の観測データが連続的な場合のモデル図のように確率密度関数である分岐確率ｃ_ik，ｃ_jkに置き換わるというものである。

（Ａ）前記第１実施例と前記第２実施例の式の変形
ここで、本実施例の具体的な内容に入る前に、離散的な場合について上記実施例と同様な説明を重複して行なう。

これは、実施例１、２と実施例３〜８との間において、数式の簡略化や符号の変更等が入っており、その点を踏まえて整合性の取れた説明をするという便宜上の理由のためである。なお、上記実施例１、２と実質的な差異はない。

まず、ソースデータとマルコフモデルの各種パラメータについて説明すると、HMMによるモデル化において、次式のようなソースデータ列が与えられる。ここで、tはデータの順序を表すインデックスである。

各ｙ_tは、スカラーあるいはベクトルである。上式は、単一配列の場合であるが、M本の複数配列で与えられるならば、ソースデータ列は次式の通りである。記号｛｝は集合を示している。

まず、単一配列では、数９８の観測データｙ（ｙは太字）が与えられる。この場合のHMMの課題は、次式の最尤推定法（Maximum Likelihood Estimation：MLE）の認識において、最良なモデルを見つけることである。

ここで、S（Sは太字）は、次式のような状態遷移系列を意味している。

ランダムな変数とそれらの値のために、小文字のｙ（ｙは太字）およびｓ（ｓは太字）は同様に用いられる（通常、ランダム変数は大文字で示される）。

各確率は以下の数１０２〜数１０５の通りである。

・初期状態確率

・状態遷移確率

ここで、もし、状態ｓ_ｉから状態ｓ_ｊへのつながりが無ければ、状態ｉから状態ｊに移る状態遷移確率はａ_ij＝０となる。これは、状態遷移のための前のトポロジーを反映している。

・出力確率

ここで、上式の各要素は次式の通りである。

上記数１０２〜数１０５に示した確率の集合をまとめてθとして次式に示す。

そして、確率的なデータ構造について以下のように解釈している。不完全なデータはｙ（ｙは太字）であり、推定される消失データはｓ（ｓは太字）であり、完全なデータは、数１００の確率量と数１０６のパラメータを有しており、次式に示されるｘ（ｘは太字）として定義する。

・Alpha-EMアルゴリズム
不完全なデータ、消失データ、完全なデータによるHMMの解釈は、EMアルゴリズムに匹敵する。本紙の目的は、新たなHMM推定アルゴリズムを見つけることであるため、alpha-EMアルゴリズムから始まる道筋を示している。

観測データｙ（ｙは太字）は推定すべきパラメータの全体を表すψによってパラメータ化されるため、Ｐ_y|ψ（ｙ｜ψ）（ｙは太字）を確率密度あるいは確率量とする。ｘ∈Ｘ（ｘとＸは太字）を、消失データを含む理想的な観測結果である完全なデータ、あるいは拡張されたデータとする。そして、不完全なデータ確率密度関数(pdf)、あるいは確率量関数(pmf)は次式で示される。

ここで、積分する範囲を次式とする。

pdfの積分は、pmfの総和となる。そして、条件付きのpdfあるいはpmfは次式である。

alpha-EMアルゴリズムでは、次式のようなalpha対数が用いられる。

ここで、α＝−１の場合は対数、すなわち、次式である。

alpha-EMアルゴリズムは、alpha対数に関して、不完全なデータの尤度比を考慮する必要があるため、次式のようになる。

ここで、φとψは、繰返し最大化ステップにおける数１０６に対して、古いモデルのパラメータと新しいモデルのパラメータを示している。それから、alpha-EMアルゴリズムの基本的な方程式が次式として得られる。

ここで、D(α)は、常に負とならない２つの条件付きの確率P_ｘ|ｙ,φ（x|y,φ）（ｘとｙは太字）とP_ｘ|ｙ,ψ（x|y,ψ）（ｘとｙは太字）の間のalphaダイバージェンスである。上式において、次式のＱ関数が重要である。Ｅは定量化（quantity）によって表される最大化のための目的関数を示している。

数１１４のために、もし、このＱ関数が正ならば、a＜1の範囲で、数１１４の左辺である不完全なデータのalpha対数尤度比も正となる。したがって、alpha-EMアルゴリズムとその変化形であるalpha-GEMアルゴリズムについて以下に示す。

[Alpha-EMアルゴリズム]
初期化: 数１０６の初期値を選択し、φとして用いる。
E-Step: 数１１５の計算を実行する。
M-step: 更新パラメータを次式で算出する。

なお、"arg max"とは最大値を与える変数を意味する記号である。すなわち、次式は、Ｑ^（α）（ｘ|ｙ、φ）（ψ｜φ）（ｘとｙは太字であり、（ｘ|ｙ、φ）はＱの下付き）を変数ψの関数と考えたとき、その最大値を与える変数の値をψ^*とする、という意味である。

U-step: φをψ^*によって置き換え、収束することを確認する。収束しない場合は、E-stepに戻って繰り返しが繰返される。

[Alpha-GEMアルゴリズム]
これは、上記M-stepを次式のようにＱ関数を非負とするψ⁺の算出に置き換えたアルゴリズムである。

まず、alpha-HMMアルゴリズムの近似バージョンがalpha-GEMアルゴリズムであることに着目すべきである。数１１４を有するalpha-EMアルゴリズムの長所によるalpha-HMM推定アルゴリズムの基本的な特性について説明する。

完全なデータをｘ＝（s，ｙ）（ｘ、ｓ、ｙは太字）とし、ｓ（ｓは太字）を消失データとし、ｙ（ｙは太字）を不完全なデータとし、数１１５は次式に相当するものとする。

もし、上式の値が負にならなければ、以下の不等式が成り立つ。

ここで、alpha-HMMの抽象バージョンに留意すべきである。

・単一配列のalpha-HMM
Non-Causalな更新方程式
離散型のアルファベット系列ｙ（ｙは太字）の場合において、最大化されるQ関数は次式で示される。なお、ｌはアルゴリズムの繰返しの指標（インデックス）である。一方、上記文章中に出てくるｔはデータの順番に対応する指標（インデックス）である。したがって、ｌとｔは区別されるものである。

まず、状態遷移確率a_ijの更新方程式について説明する。状態遷移確率a_ij は、その更新後に確率量とする必要があるので、ラグランジュの未定係数法（Lagrange multiplier）を使用する必要がある。したがって、次式のように最大値を算出するための微分を行なう。

そうすると、上式から次式が得られる。

ここで、Ｎ_ij(ｓ)（ｓは太字）は、iからjへの状態遷移の数と位置である。次のステップでは、ソフトウェアで実行可能なアルゴリズムを目指すため、以下の問題を解決する必要がある。
(a)Non-causalityの存在: 右辺にθ_l+1(θは太字）が含まれている。このままでは、自己撞着性のため計算できない。
(b)右辺の算出には、状態数がＮでデータ数がＴのとき、Ｏ(Ｎ^T)すなわちＮのＴ乗のオーダーの演算が必要である。

上記の問題は、次のセクションで解決されることとなるが、その前に、出力確率b_jk|θl+1（l+1はθの添え字）と初期状態確率π_i|θl+1（l+1はθの添え字）の２つの更なる更新式について次式に示す。なお、Ｎ_jk(ｓ)（ｓは太字）は、状態ｓ_t＝ｊでの出力y_t＝kの発生回数である。

Causal近似と拡張系列: 離散出力の場合
更新方程式である数１２２〜数１２４のコアとなる部分は、繰返し指標（インデックス）のシフトによるCausal近似によって、次式のように変換することができる。

したがって、P（ｙ|θ_l）＝P（ｙ|θ_l-1）＋o(1)(yは太字）の領域で次の等式を得る。

ここで、α_causalという表記は、以後、βと表記する。

いま、数１２５が算出可能であるが、演算の複雑度のために、もう一つの近似が必要となる。このため、我々は次式の系列拡張を用いている。

そして、causal近似式である数１２５の適用と数１２７の系列拡張は、以下に示す遷移確率の更新方程式を与える。なお、これらの確率の更新式である上記数１２２〜数１２４は、それぞれ上記実施例１の数７５〜７７に対応するものである。

単一配列の離散型alpha-HMMの遷移確率a_ij|θl+1(l+1は、θの添え字)は、次式である。

ここで、以下の各特性について理解しておくことが重要である。
[特性１]
β＝α_causal＝１の場合は、従来のlog-HMMの方法に帰趨する。
[特性２]
数１２８の分子は、現在と過去の更新項の重み付けられた総和である。分母も同様である。
[特性３]
数１２８の第２行及び第３行は、これらの確率計算に伴う複雑度を抑えた従来の前方-後方アルゴリズムの方法に匹敵する。
[特性４]
alpha-HMMに唯一、追加する必要があるものは、θ_l-1での更新項を記憶することである。これは、CPU時間に直接現れているように繰返しを減らすことを意味する。実際の実験結果でも、この予想通りになっている。

次に、残りの２つの更新方程式を示す。

単一配列の離散型alpha-HMMの出力確率ｂ_jk|θl+1(l+1は、θの添え字)は、次式である。

単一配列の離散型alpha-HMMの初期状態確率π_i|θl+1(l+1は、θの添え字)は、次式である。

（変形例）
以下、近似式の変形例について説明する。なお、以下の変形例は本発明の他の実施例でも同様に適用可能なものである。上述した隠れマルコフモデルの推定法では、有限個の過去情報を用いて各確率量をすべて次のような計算に基づいて計算している。

この高速化パラメータβは、利用者が選ぶデザインパラメータであり、現在値による繰り返し値に対し、1回前の繰り返し値を重みづけする効果を有している。この過去の値に重みづけを行って加算することは、有限の過去にまでさかのぼって行うこともできる。すなわち、ε_τ≧０をデザインパラメータとして設定し、次式のように置換してもよい。

そうすると、数１３１は次式となる。

次に展開による近似式を用いる場合について説明すると、上記数１３１や数１３４は、次式のように展開できる（例えばテーラー展開）。

ここで、o(1)は高位の無限小を意味する記号であり、上式の右辺第一項は，従来の隠れマルコフモデル推定アルゴリズムに相当している。そして、モーメンタム項に相当する右辺第二項は、展開により近似される高速化項の意味をもつ。

なお、Ｔ＝１かつε_１＝Ｐ（ｙ｜θ_l-1）／Ｐ（ｙ｜θ_l）（ｙはそれぞれ太字），すなわち重みを用いた場合は、重みを用いていない場合に比べて、収束速度はわずかに遅くなるものの、β＝３付近での収束性能の安定性が良くなる。換言すると、図３のような性能曲線において、計算が収束する立ち上がり部分では、高速化に伴い、波状の上下の変動が生じやすくなるが、その変動が少なくなる。

・複数配列のalpha-HMM
もし、HMMの予め設計されたトポロジーは、エルゴードなものであれば、単一の長いトレーニング系列ｙ（ｙは太字）は、十分である。もし、選択されたトポロジーが、吸収状態を有しているならば、複数のトレーニング系列を使うことが望まれる。このような複数系列に対するalpha-HMM推定の更新式は、上記の単一系列に対する方法を利用する形で得ることができる。

離散的なシンボルの場合において、Ｓ（Ｓは太字）を次式のようにＭ個の状態遷移系列の集合とする。

すると、複数配列のＱ関数は次式である。

ここで、Ｐは、次式に示されるMarkov過程の確率である。

上式において、初期状態確率π_s0(n)(0はｓの下付きの添え字であり、（ｎ）は上付きの添え字である)、状態遷移確率a_st-1(n)st(n) (t-1とtは、ｓの下付きの添え字であり、（ｎ）は上付きの添え字である)、及び出力確率ｂ_st(n)(y_t ⁽ⁿ⁾) (tは、ｓの下付きの添え字であり、（ｎ）は上付きの添え字である)の形態が、系列指標（インデックス）ｎから独立している点に着目すべきである。

初期状態確率π_i|θl+1(l+1は、θの添え字)によって、初期状態の更新方程式の誘導は、Ｑ関数の数１３７の微分から始まる。上述の単一配列との違いは、上式からも明らかなように微分がｎ回現れることである。そして、次式のようにnon-causalな方程式が得られる。

ここで、

であり、上式の分母の関数ｆ(ｎ)は次式である。

そして、causalの繰返し指標（インデックス）のシフトと系列の拡張は、次の更新方程式を与える。

同様に、状態遷移と出力の更新方程式が得られる。

ここで、

及び

また、次式の出力確率が得られる。

ここで、上式で省略された各項はそれぞれ次式となる。

上記数１４７のΣの下にある記号「ｔ：ｙ_t＝ｋ」の意味は、ｙ_t＝ｋが成り立つｔのみを対象にする、という意味である。すなわち、上記数１４７の場合、ｙ＝｛ｙ_ｔ｝^T（ｔ＝１）（左辺のyは太字であり、ｔ＝１は｛ｙ_ｔ｝の下付きの添え字）において、ｙ_t＝ｋが成り立つｔのみについて加算する、という意味である。

（Ｂ）連続系列の単一配列
ここから、連続的な出力系列の場合の本実施例の具体的な説明を行なう。

もし、出力系列ｙ＝{ｙ_t}^Ｔ（ｔ＝１）（出力系列ｙは太字であり、ｔ＝１は｛ｙ_t｝の下付き）が、連続的な多変数の観測結果として現れれば、数１２８〜数１３０と似ているが、少しだけ異なる更新方程式が得られる。この場合、ｙ_tは太字の文字で示されていないが、適切な次元のユークリッド空間におけるベクトルである。

このような連続的なアルファベットの場合、最尤推定法（MLE）の問題は、以下の尤度を最大化することである。

ここで、状態s_tで、出力確率ｂ_ｓtｋt（ｙ_t）（_ｓtｋtのtは、それぞれｓ，ｋの下付きの添え字）において、k_t番目の枝へ遷移する確率を特定する分岐確率Ｃ_ｓtｋt（_ｓtｋtのtは、それぞれｓ，ｋの下付きの添え字）は、ｙ_tの確率密度関数である。我々は、これをガウス密度（Gaussian density）とみなしている。そして、状態ｊで状態ｋが出力される出力確率ｂ_ｊｋ（ｙ_t）は次式となる。

ここで、μ_ｊｋは平均値ベクトルであり、Σ_ｊｋは共分散行列であり、総和と混同してはいけない。平均値ベクトルとは、分岐確率の確率密度関数の平均値を示している（以後、同様）。そして、状態jでの出力確率密度関数(pdf)は、ｂ_ｊ（ｙ_t）がpdfとなるように次式となる。

ここで、上式の記号について説明する。

まず、離散系列の場合、ｋを出力記号の種類を表すインデックスとすると、出力y_kのとる値は有限種類であるため、ｂ_j(y_k)はｂ_ｊｋと略記できる。しかし、連続系列の場合は、離散値ではなく、連続値であるため、そのような表記はできない。そのため、時間を表すインデックスをｔとして、y_ｔはそのまま記載する必要がある。

また、N（ｙ_ｔ；μ_ｊｋ，Σ_ｊｋ）は、平均値ベクトルがμ_ｊｋであり、Σ_ｊｋを共分散行列とする多次元正規確率密度関数（多次元ガウス確率密度関数）で、ｙ_ｔを分布の変数とするという意味である。

実際、このようなガウス混合モデルは、「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でもここまでが導出可能な場合にとどまっている。

図６は、離散的なアルファベットの場合を図示している。また、図７は、ガウス混合モデル（Gaussian mixture model）(bottom)の場合、すなわち連続的なアルファベットの場合を図示している。これらの図面を参照すると、図７中の分岐確率ｃ_ｊｋの矢印は、図６中の出力確率ｂ_ｊｋの矢印と対応していることがわかる。 数１４９と数１００を参照すると、数１４９の分岐確率ｃ_ｓtｋt（_ｓtｋtのtは、それぞれｓ，ｋの下付きの添え字）と出力確率ｂ_ｓtｋt（ｙ_t）（_ｓtｋtのtは、それぞれｓ，ｋの下付きの添え字）の積が、数１００の出力確率ｂ_ｓt（ｙ_t）（_ｓtのtは、ｓの下付きの添え字）に対応しているように見えるかもしれない。しかし、数１４９において離散値シンボルとみなされるのは、出力確率ｂ_ｓtｋt（ｙ_t）（_ｓtｋtのtは、それぞれｓ，ｋの下付きの添え字）を除いた部分であるため、図６、図７の対応関係との矛盾は生じない。

混合確率の場合のため、消失データはｓ（ｓは太字）とｃ（ｃは太字）である。そのため、Ｑ関数は次式である。

数１４９、数１５２から明らかなように、初期確率と状態遷移確率の更新方程式は、それぞれ数１２８、数１３０と同様である。分岐確率ｃ_ｊｋの更新方程式は、状態遷移確率a_ijの場合と同様にラグランジュの未定係数法（Lagrange multiplier）によって得られる。

そして、右辺は、状態遷移確率の場合と同様にcausalで計算可能となり、次式のように変形される。

次式は、平均値ベクトルμ_ｊｋの更新方程式である。μ_ｊｋについて数１５２の直接微分として、次のnon-causalの方程式が得られる。

そして、繰返し指標（インデックス）のシフト、系列の拡張、総和の変更により、次式が得られる。

数１５６が、θ_l-1の過去情報が十分に利用できる点を示していることに着目すべきである。

共分散行列の更新には、行列微分が必要である。Ｑ関数の数１５２を共分散行列の逆行列Σ^-1 _ｊｋ（ｊｋはΣの添え字）について微分することにより、次のnon-causalな方程式が得られる。

そして、繰返し指標（インデックス）のシフト、系列の拡張、および総和の変更により、次の更新式が得られる。

ここで、上式の各項は以下の通りである。

共分散行列の更新式である数１５８には、効果的に十分に利用できる過去情報の形態を有していることに着目すべきである。なお、記号Σ^-1 _{ｊｋ｜θl;1}とは、l+1回目の状態ｊで分岐ｋに移行したときの出力の共分散行列を表している（以後、同様）。

ここで、ガウス混合（Gaussian mixture）alpha-HMM、すなわち単一系列の連続型alpha-HMMの更新方法について簡単に述べる。
[単一系列の連続型alpha-HMMの初期状態確率]
更新式は、数１３０である。
[単一系列の連続型alpha-HMMの状態遷移確率]
更新式は、数１２８である。
[単一系列の連続型alpha-HMMの分岐確率]
更新式は、数１５４である。
[単一系列の連続型alpha-HMMの平均値ベクトル]
更新式は、数１５６である。
[単一系列の連続型alpha-HMMの共分散行列]
更新式は、数１５８であり、その要素は、数１５９、数１６０である。

ここで、θ_l-1のすべての情報をメモリに記憶する点を改めて強調することが重要である。θ_lによって指標（インデックス）を付された項の算出は、log-HMMと同等である。

図８は、本発明の第３実施例におけるプログラムの処理手順をあらわしたものである。このフローチャート全体の流れは、前記第１実施例の場合とほぼ同様である。ただし、本実施例では、観測データｙが、連続的な多変数の観測結果として現れる単一配列の連続系列データである点で異なる。また、初期値として設定され、繰返し計算される未知パラメータには、初期状態確率、状態遷移確率、分岐確率、平均値ベクトル、共分散行列があるため、各ステップにおける計算式も異なる。

そのため、置き換わる数式の対応関係について説明する。まず、繰返し更新される未知パラメータの集合の組み合わせを次式とする。

推定装置１０が最終的に算出しようとするHMMの確率構造は、数６６の代わりに数１４９となる。ステップＳ３２で決定する未知パラメータの初期値は、数６７の代わりに次式となる。

ステップＳ３２における収束判定値の決定方法は、数８３と同じである。

ステップＳ３３では、最初にステップＳ３２で決定した未知パラメータの初期値を利用し、それ以降はステップＳ３９で更新された未知パラメータとカウント値を利用して、次の数１６３に示す各確率値と、数１６４に示す各カウント値を、実際に推定装置１０のメモリ（図示せず）に設定する。なお、前記実施例１、２と異なり、指標（インデックス）は、データの順番のｔではなく、アルゴリズムの繰り返しのｌを付して設定する。

その後のステップＳ３４〜Ｓ３６の演算動作は、前記実施例１の数７０〜数７４の一部の符号を置き換えて適用することによって行う。具体的には、出力確率ｂ_{ｊ．ｙτ+1}（τ＋１はｙの下付き添え字）を分岐確率ｃ_{ｊ．ｙτ+1}（τ＋１はｙの下付添え字）と置き換え、指標（インデックス）であるｔをｌと置き換える。

ステップＳ３７では、l+1回目に繰り返される確率π，a，ｃと、平均値ベクトルμ，及び共分散行列Σの組み合わせを条件として、状態iとなる初期状態確率
初期状態確率π_ｉ｜θl+1（l+1はθの添え字）の値と、状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θl+1（l+1はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂ_jk（ｙ_t）において、分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃ_ｊｋ|l+1（l+1はθの添え字）の値と、平均値ベクトルμ_jｋ|θl+1（l+1はθの添え字）の値と，共分散行列Σ_jｋ|θl+1（l+1はθの添え字）の値と，ステップＳ３４で求めたl回目の確率Ｐ（ｙ｜θ_l）と、それよりも１回前の確率Ｐ（ｙ｜θ_l-1）を利用し、且つステップＳ３５で求めたl回目のカウント値Ｎ_ａij｜_θl，Ｎ_cjk｜_θl（ijはａの添え字、jkはcの添え字、lはθの添え字）と、それよりも１回前のカウント値Ｎ_ａij｜_θl-1，Ｎ_cjk｜_θl-1（ijはａの添え字、jkはｃの添え字、l-1はθの添え字）を利用して計算する部分であり、これは前記数１２８、数１３０、数１５４、数１５６、数１５８〜数１６０の各更新式にそれぞれ対応している。

このステップＳ３７においても高速化パラメータの値βが用いられ、添え字l+1の部分の確率を計算する上で、時間シフトした添え字lや添え字l-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したこと等は、前記実施例１と同様にその効果として挙げられる。

本実施例では、配列が１本の観測データをレジスタ１に格納しているが、その場合のステップＳ３７における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数１３０、状態遷移確率の更新式は数１２８、分岐確率の更新式は数１５４、平均値ベクトルの更新式は数１５６、共分散行列の更新式は数１５８〜数１６０である。

ステップＳ３８では、前記実施例１と同様に、前記ステップ３７で計算された新たな確率量に基づく尤度Ｐを用いて、その尤度ＰとステップＳ３２における収束判定値との比較により、ステップＳ３３〜Ｓ３９の計算が収束したか否かを判定する。収束してなければ、ステップＳ３９に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように２つの過去値をシフトさせる。

ここでは、繰り返しの回数が１つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。なお、ステップＳ３３とステップＳ３９を一纏めにして、更新された確率量およびカウント値をそのまま次の繰り返しのために、推定装置１０のメモリに設定記憶させてもよい。

一方、ステップＳ３８において計算が収束していると判定された場合、ステップＳ４０に移行して、次式に示す計算された５組の確率量を用いたＨＭＭを採用する。

上記数１６６で計算した各値を用いて、推定装置１０は数１４９に示すＨＭＭの確率構造を出力することができる。

図８に示すアルゴリズムを実現する推定装置の構成は、図２に示した前記第１実施例の推定装置１０と同様であるため、ここでは説明を省略する。

以上のように本実施例では、入力される観測データを時系列に格納する記憶手段としてのレジスタ１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１４とを備えたＨＭＭの推定装置１０において、推定手段１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび分岐の期待値Ｎ_ｃを更新設定する更新設定手段２４と、更新設定手段２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数１６３，数１６４を参照）と共に、前記レジスタ１から読み出した観測データと、初期設定手段２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１２８，数１３０，数１５４，数１５６，数１５８〜数１６０を参照）演算手段２６と、演算手段２６による計算の収束を判定し、計算が収束していなければ、演算手段２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段２６で計算した新たな各確率量を最終的な値として出力させる判定手段２８とを備えている。

このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび分岐の期待値Ｎ_ｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。

また、ここでの演算手段２６は、l+1回目の状態ｉとなる初期状態確率π_ｉ｜θl+1（l+1はθの添え字）の値を、前記数１３０の式で計算し、l+1回目の状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における状態遷移の期待値Ｎ_ａij｜_θlおよびＮ_ａij｜_θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、前記数１２８の式で計算し、状態ｊで状態ｋが出力される出力確率ｂ_jk｜θl+1（l+1はθの添え字）の際に分岐ｋに移行する枝へ遷移する確率を特定する分岐確率ｃ_jk｜θl+1の値と、平均値ベクトルμ_jk｜θl+1（l+1はθの添え字）の値と、共分散行列Σ_jk｜θl+1の値とを、l回目およびl-1回目における前記出力の期待値Ｎ_cjk｜_θlおよびＮ_cjk｜_θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、前記数１５４、数１５６、数１５８〜数１６０の式で計算する構成となっている。

そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。

さらに、ここでの更新設定手段２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。

なお上述した本実施例の作用効果は、初期設定手段２２としての動作を実行する初期設定ステップと、更新設定手段２４としての動作を実行する更新設定ステップと、演算手段２６としての動作を実行する演算ステップと、判定手段２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。

また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。

本実施例では、連続的なシンボルで複数配列のデータ系列の場合の実施形態について説明する。

本実施例においても、上記実施例と同様に、繰返し更新される未知パラメータの集合の組み合わせを数１６１とし、連続的なシンボルのalpha-HMMの更新方程式が得られる。初期状態確率と状態遷移確率の更新方程式は、それぞれ数１４２〜数１４５と同様である。しかし、出力の更新方程式の集合は、離散的なシンボルの場合と異なる。我々は、分岐確率、平均値ベクトル、及び共分散行列の更新方程式を必要としている。

分岐確率については、初期状態確率と状態遷移確率の場合と同様にラグランジュの未定係数法を用いて算出可能である。そして、次式の更新方程式が得られる。

ここで、各要素は次式である。

Ｑ関数の数１３７での直接ベクトル微分であるμ_{ｊｋ｜θl+1}（l+1はθの下付きの添え字）が適用される。その更新方程式は次式である。

ここで、

数１７０、数１７１と同様に、共分散行列の更新方程式は、Σ^-1 _ｊｋ｜θl（ｊｋ｜θlはΣの下付き添え字であり、さらにlはθの下付き添え字）について行列微分を用いることによって得られる。

ここで、Ｇ⁽ⁿ⁾ _θl（θlはＧの下付き添え字）は次式である。

図９は、本発明の第４実施例におけるプログラムの処理手順をあらわしたものである。このフローチャート全体の流れは、前記第３実施例の場合とほぼ同様である。ただし、本実施例では、観測データｙが、連続的な多変数の観測結果として現れる複数配列の連続系列データである点で異なる。

置き換わる数式の対応関係について説明する。まず、推定装置１００が最終的に算出しようとするHMMの確率構造は、数６６の代わりに数１４９となる。ステップＳ４２で決定する未知パラメータの初期値は、数１６２と同じである。

ステップＳ４２における収束判定値の決定方法は、数８３と同じである。

ステップＳ４３では、最初にステップＳ４２で決定した未知パラメータの初期値を利用し、それ以降はステップＳ４９で更新された未知パラメータとカウント値を利用して、数１７４に示す各確率値と、数１７５に示す各カウント値を、実際に推定装置１００のメモリ（図示せず）に設定する。なお、指標（インデックス）は、データの順番のｔではなく、アルゴリズムの繰り返しのｌを付して設定する。

その後のステップＳ４４〜Ｓ４６の演算動作は、前記実施例２の数８８〜数９２の一部の符号を置き換えて適用することによって行う。具体的には、出力確率ｂ_{ｊ．ｙτ+1}（τ＋１はｙの下付添え字）を分岐確率ｃ_{ｊ．ｙτ+1}（τ＋１はｙの下付添え字）と置き換え、指標（インデックス）であるｔをｌと置き換える。

ステップＳ４７では、l+1回目に繰り返される確率π，a，ｃと、平均値ベクトルμ，及び共分散行列Σの組み合わせを条件として、状態iとなる初期状態確率π_ｉ｜θl+1（l+1はθの添え字）の値と、状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θl+1（l+1はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂ_jk（ｙ_t）において、分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃ_ｊｋ|l+1（l+1はθの添え字）の値と、平均値ベクトルμ_jｋ|θl+1（l+1はθの添え字）の値と，共分散行列Σ_jｋ|θl+1（l+1はθの添え字）の値とを、ステップＳ３４で求めたl回目の確率Ｐ（ｙ^（n）｜θ_l）と、それよりも１回前の確率Ｐ（ｙ^（n）｜θ_l-1）を利用し、且つステップＳ３５で求めたl回目のカウント値Ｎ^（n） _ａij｜_θl，Ｎ^（n） _cjk｜_θl（ijはａの添え字、jkはcの添え字、lはθの添え字）と、それよりも１回前のカウント値Ｎ^（n） _ａij｜_θl-1，Ｎ^（n） _cjk｜_θl-1（ijはａの添え字、jkはｃの添え字、l-1はθの添え字）を利用して計算する部分であり、これは前記数１４２〜数１４５、数１６７〜数１７３に対応している。

このステップＳ４７においても高速化パラメータの値βが用いられ、添え字l+1の部分の確率を計算する上で、時間シフトした添え字lや添え字l-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したこと等は、前記実施例２と同様にその効果として挙げられる。

本実施例では、配列がＭ本の観測データをレジスタ１０１に格納しているが、その場合のステップＳ４７における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数１４２、状態遷移確率の更新式は数１４３〜数１４５、分岐確率の更新式は数１６７〜数１６９、平均値ベクトルの更新式は数１７０、数１７１、共分散行列の更新式は数１７２、数１７３である。

ステップＳ４８では、前記実施例２と同様に、前記ステップＳ４７で計算された新たな確率量に基づく尤度Ｐを用いて、その尤度ＰとステップＳ４２における収束判定値との比較により、ステップＳ４３〜Ｓ４９の計算が収束したか否かを判定する。収束してなければ、ステップＳ４９に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように２つの過去値をシフトさせる。

ここでは、繰り返しの回数が１つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。

一方、ステップＳ４８において計算が収束していると判定された場合、ステップＳ５０に移行して、次式に示す計算された５組の確率量を用いたＨＭＭを採用する。

上記数１７７で計算した各値を用いて、推定装置１０１は数１４９に示すＨＭＭの確率構造を出力することができる。

図９に示すアルゴリズムを実現する推定装置の構成は、図５に示した前記第２実施例の推定装置１０１と同様であるため、ここでは説明を省略する。

以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１０１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１１４とを備えたＨＭＭの推定装置１００において、推定手段１１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段１２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび分岐の期待値Ｎ_ｃを更新設定する更新設定手段１２４と、更新設定手段１２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数１７４，数１７５を参照）と共に、前記レジスタ１０１から読み出した観測データと、初期設定手段１２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１４２〜数１４５，数１６７〜数１７３を参照）演算手段１２６と、演算手段１２６による計算の収束を判定し、計算が収束していなければ、演算手段１２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段１２６で計算した新たな各確率量を最終的な値として出力させる判定手段１２８とを備えている。

このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび分岐の期待値Ｎ_ｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。

また、ここでの演算手段１２６は、l+1回目の状態ｉとなる初期状態確率π_ｉ｜θl+1（l+1はθの添え字）の値を、前記数１４２の式で計算し、l+1回目の状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θl+1（t+1はθの添え字）の値を、l回目およびl-1回目における状態遷移の期待値Ｎ_ａij｜_θlおよびＮ_ａij｜_θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、前記数１４３〜数１４５の式で計算し、状態ｊで状態ｋが出力される出力確率ｂ_jk｜θl+1（l+1はθの添え字）の際に分岐ｋの枝へ遷移する確率を特定する分岐確率ｃ_jk｜θl+1の値と、平均値ベクトルμ_jk｜θl+1（l+1はθの添え字）の値と、共分散行列Σ_jk｜θl+1の値とを、l回目およびl-1回目における前記分岐の期待値Ｎ_cjk｜_θlおよびＮ_cjk｜_θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、前記数１６７〜数１７３の式で計算する構成となっている。

そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。

さらに、ここでの更新設定手段１２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。

なお上述した本実施例の作用効果は、初期設定手段１２２としての動作を実行する初期設定ステップと、更新設定手段１２４としての動作を実行する更新設定ステップと、演算手段１２６としての動作を実行する演算ステップと、判定手段１２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。

本実施例では、半連続的なシンボルで単一配列のデータ系列の場合の実施形態について説明する。

図７において、ガウス混合（Gaussian mixture）alpha-HMM（log-HMMも）のグラフ構造を見直すことにより、以下を実現する。
（a）ガウス混合HMMにおいて、各ガウス（Gaussian）pdfは、到着する状態jに依存する。すべてのN×Kのガウス密度の学習には、多様な長いトレーニング系列を必要とする。
(b)離散型の場合のb_jkの役割を連続モードのc_jkに割り当てる。また、平均値ベクトルと共分散行列は遷移状態jに依存しないという場合を考慮する。そして、図６の離散型の場合をまっすぐ延長した構造となる。これを半連続HMMと呼んでいる。この構造のもう一つの解釈は、MLE-VQ HMM (Maximum Likelihood Vector Quantization HMM)である。

半連続のalpha-HMMモデルは、数１４９を変形し、次式となる。

したがって、半連続のalpha-HMMの更新方程式は以下のようになる。
[単一配列の半連続のalpha-HMMの初期状態確率]
更新式は、数１３０である。
[単一配列の半連続のalpha-HMMの状態遷移確率]
更新式は、数１２８である。
[単一配列の半連続のalpha-HMMの分岐確率]
更新式は、数１５４である。
[単一配列の半連続のalpha-HMMの平均値ベクトル]
次式のように、数１５６において、μ_{ｊｋ｜θl+1}（l+1はθの下付きの添え字）をμ_ｊ｜θl+1（l+1はθの下付きの添え字）とし、右辺のk_t=kの項を除去したものである。これは、状態ｊに依存しないためである（以後、同様である）。

[単一配列の半連続のalpha-HMMの共分散行列]
次式のように、数１５８〜数１６０において、Σ_{ｊｋ｜θl+1}（l+1はθの添え字），μ_ｊｋ｜θl（lはθの添え字），μ_{ｊｋ｜θl-1}（l-1はθの添え字）を、それぞれΣ_ｊ｜θl+1（l+1はθの添え字），μ_ｊ｜θl（lはθの添え字），μ_ｊ｜θl-1（l-1はθの添え字）とし、右辺のk_t=kの項を除去したものである。

ここで、各要素は次式である。

本発明の第５実施例におけるプログラムの処理手順は、前記第３実施例の場合とほぼ同様であるため、図８を用いて説明する。ただし、本実施例では、観測データｙが、連続的な多変数の観測結果として現れる単一配列の半連続系列データである点で異なる。

本実施例は、前記第３実施例の場合の特例と考えることもできるので、置き換わる数式の対応関係に絞って説明する。まず、繰返し更新される未知パラメータの集合の組み合わせを数１６１とし、推定装置１０が最終的に算出しようとするHMMの確率構造は、数６６の代わりに数１７８とする。ステップＳ３２で決定する未知パラメータの初期値は、数６７の代わりに次式となる。

ステップＳ３２における収束判定値の決定方法は、数８３と同じである。

ステップＳ３３では、最初にステップＳ３２で決定した未知パラメータの初期値を利用し、それ以降はステップＳ３９で更新された未知パラメータとカウント値を利用して、次の数１８４に示す各確率値と、数１８５に示す各カウント値を、実際に推定装置１０のメモリ（図示せず）に設定する。

その後のステップＳ３４〜Ｓ３６の演算動作は、前記実施例３と同様に前記実施例１の数７０〜数７４の一部の符号を置き換えて適用することによって行なう。

ステップＳ３７では、l+1回目に繰り返される確率π，a，ｃと、平均値ベクトルμ，及び共分散行列Σの組み合わせを条件として、状態iとなる初期状態確率π_ｉ｜θl+1（l+1はθの添え字）の値と、状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θl+1（l+1はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂ_jk（ｙ_t）において、分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃ_ｊｋ|l+1（l+1はθの添え字）の値と、平均値ベクトルμ_j|θl+1（l+1はθの添え字）の値と，共分散行列Σ_j|θl+1（l+1はθの添え字）の値とを、ステップＳ３４で求めたl回目の確率Ｐ（ｙ｜θ_l）と、それよりも１回前の確率Ｐ（ｙ｜θ_l-1）を利用し、且つステップＳ３５で求めたl回目のカウント値Ｎ_ａij｜_θl，Ｎ_cjk｜_θl（ijはａの添え字、jkはcの添え字、lはθの添え字）と、それよりも１回前のカウント値Ｎ_ａij｜_θl-1，Ｎ_cjk｜_θl-1（ijはａの添え字、jkはｃの添え字、l-1はθの添え字）を利用して計算する部分であり、これは前記数１３０、数１２８、数１５４、数１７９〜数１８２に対応している。

このステップＳ３７においても高速化パラメータの値βが用いられ、添え字l+1の部分の確率を計算する上で、時間シフトした添え字lや添え字l-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したこと等は、前記実施例３と同様にその効果として挙げられる。

本実施例では、配列が１本の観測データをレジスタ１に格納しているが、その場合のステップＳ３７における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数１３０、状態遷移確率の更新式は数１２８、分岐確率の更新式は数１５４、平均値ベクトルの更新式は数１７９、共分散行列の更新式は数１８０〜数１８２である。

ステップＳ３８では、前記実施例３と同様に、前記ステップ３７で計算された新たな確率量に基づく尤度Ｐを用いて、その尤度ＰとステップＳ３２における収束判定値との比較により、ステップＳ３３〜Ｓ３９の計算が収束したか否かを判定する。収束してなければ、ステップＳ３９に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように２つの過去値をシフトさせる。

ここでは、繰り返しの回数が１つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。なお、ステップＳ３３とステップＳ３９を一纏めにして、更新された確率量およびカウント値をそのまま次の繰り返しのために、推定装置１０のメモリに設定記憶させてもよい。

一方、ステップＳ３８において計算が収束していると判定された場合、ステップＳ４０に移行して、次式に示す計算された５組の確率量を用いたＨＭＭを採用する。

上記数１８７で計算した各値を用いて、推定装置１０は数１７８に示すＨＭＭの確率構造を出力することができる。

図８に示すアルゴリズムを実現する推定装置の構成は、図２に示した前記第３実施例の推定装置１０と同様であるため、ここでは説明を省略する。

以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１４とを備えたＨＭＭの推定装置１０において、推定手段１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび分岐の期待値Ｎ_ｃを更新設定する更新設定手段２４と、更新設定手段２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数１８４，数１８５を参照）と共に、前記レジスタ１から読み出した観測データと、初期設定手段２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１３０，数１２８，数１５４，数１７９〜数１８２を参照）演算手段２６と、演算手段２６による計算の収束を判定し、計算が収束していなければ、演算手段２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段２６で計算した新たな各確率量を最終的な値として出力させる判定手段２８とを備えている。

このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび分岐の期待値Ｎ_ｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。

また、ここでの演算手段２６は、l+1回目に状態ｉとなる初期状態確率π_ｉ｜θl+1（l+1はθの添え字）の値を、前記数１３０の式で計算し、l+1回目の状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θl+1（t+1はθの添え字）の値を、l回目およびl-1回目における状態遷移の期待値Ｎ_ａij｜_θlおよびＮ_ａij｜_θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、前記数１２８の式で計算し、状態ｊで状態ｋが出力される出力確率ｂ_jk｜θl+1（l+1はθの添え字）の際に分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃ_jk｜θl+1の値と、平均値ベクトルμ_jk｜θl+1（l+1はθの添え字）の値と、共分散行列Σ_jk｜θl+1の値とを、l回目およびl-1回目における前記出力の期待値Ｎ_cjk｜_θlおよびＮ_cjk｜_θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、前記数１５４、数１７９〜数１８２の式で計算する構成となっている。

そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。

さらに、ここでの更新設定手段２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。

なお上述した本実施例の作用効果は、初期設定手段２２としての動作を実行する初期設定ステップと、更新設定手段２４としての動作を実行する更新設定ステップと、演算手段２６としての動作を実行する演算ステップと、判定手段２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。

本実施例では、半連続的なシンボルで複数配列のデータ系列の場合の実施形態について説明する。

この場合の複数配列の更新方程式は、平均ベクトルと共分散の状態依存性を制限することによって得られる。初期状態確率と状態遷移確率と分岐確率の更新方程式は、数１４２〜数１４５、および数１６７〜数１６９と同様である。したがって、次式のようになる。ここで、各記号の上にバー(-)を付し、バー(-）のない場合と記号を区別している。これは、例えば、数１７１と以下の数１８９では数式が異なっており、同じ記号が使えないためである。すなわち、バー(-）それ自体に特別な意味はない。

また、本実施例では、前記第５実施例において添え字がｊだった記号は、添え字がｋになっている。これは、k_t=kの項の代わりにｓ_ｔ＝ｊの項が除去されたためである。（以後、同様である）。

上記数１８８と同様に、共分散行列の更新方程式は、状態依存性を除去することによって次式のように得られる。

本発明の第６実施例におけるプログラムの処理手順は、前記第４実施例の場合とほぼ同様であるため、図９を用いて説明する。ただし、本実施例では、観測データｙが、連続的な多変数の観測結果として現れる複数配列の半連続系列データである点で異なる。

本実施例は、前記第４実施例の場合の特例と考えることもできるので、置き換わる数式の対応関係に絞って説明する。まず、繰返し更新される未知パラメータの集合の組み合わせを数１６１とし、推定装置１００が最終的に算出しようとするHMMの確率構造は、数６６の代わりに数１７８とする。ステップＳ４２で決定する未知パラメータの初期値は、数１８２と同じである。

ステップＳ４２における収束判定値の決定方法は、数８３と同じである。

ステップＳ４３では、最初にステップＳ４２で決定した未知パラメータの初期値を利用し、それ以降はステップＳ４９で更新された未知パラメータとカウント値を利用して、次の数１９２に示す各確率値と、数１９３に示す各カウント値を、実際に推定装置１００のメモリ（図示せず）に設定する。

その後のステップＳ４４〜Ｓ４６の演算動作は、前記実施例４と同様に前記実施例２の数８８〜数９２の一部の符号を置き換えて適用することによって行なう。

ステップＳ４７では、l+1回目に繰り返される確率π，a，ｃと、平均値ベクトルμ，及び共分散行列Σの組み合わせを条件として、状態iとなる初期状態確率π_ｉ｜θl+1（l+1はθの添え字）の値と、状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θl+1（l+1はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂ_jk（ｙ_t）において、分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃ_ｊｋ|l+1（l+1はθの添え字）の値と、平均値ベクトルμ-_k|θl+1（l+1はθの添え字であり、-はμの上に付く）の値と，共分散行列Σ_k|θl+1（l+1はθの添え字であり、-はΣの上に付く）の値とを、ステップＳ３４で求めたl回目の確率Ｐ（ｙ^（n）｜θ_l）と、それよりも１回前の確率Ｐ（ｙ^（n）｜θ_l-1）を利用し、且つステップＳ３５で求めたl回目のカウント値Ｎ^（n） _ａij｜_θl，Ｎ^（n） _cjk｜_θl（ijはａの添え字、jkはcの添え字、lはθの添え字）と、それよりも１回前のカウント値Ｎ^（n） _ａij｜_θl-1，Ｎ^（n） _cjk｜_θl-1（ijはａの添え字、jkはｃの添え字、l-1はθの添え字）を利用して計算する部分であり、これは前記数１４２〜数１４５、数１６７〜数１６９、数１８８〜数１９１に対応している。

このステップＳ４７においても高速化パラメータの値βが用いられ、添え字l+1の部分の確率を計算する上で、時間シフトした添え字lや添え字l-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したこと等は、前記実施例４と同様にその効果として挙げられる。

本実施例では、配列がＭ本の観測データをレジスタ１０１に格納しているが、その場合のステップＳ４７における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数１４２、状態遷移確率の更新式は数１４３〜数１４５、分岐確率の更新式は数１６７〜数１６９、平均値ベクトルの更新式は数１８８、数１８９、共分散行列の更新式は数１９０、数１９１である。

ステップＳ４８では、前記実施例４と同様に、前記ステップＳ４７で計算された新たな確率量に基づく尤度Ｐを用いて、その尤度ＰとステップＳ４２における収束判定値との比較により、ステップＳ４３〜Ｓ４９の計算が収束したか否かを判定する。収束してなければ、ステップＳ４９に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように２つの過去値をシフトさせる。

ここでは、繰り返しの回数が１つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。

一方、ステップＳ４８において計算が収束していると判定された場合、ステップＳ５０に移行して、次式に示す計算された５組の確率量を用いたＨＭＭを採用する。

上記数１９５で計算した各値を用いて、推定装置１０１は数１７８に示すＨＭＭの確率構造を出力することができる。

図９に示すアルゴリズムを実現する推定装置の構成は、図５に示した前記第２実施例の推定装置１０１と同様であるため、ここでは説明を省略する。

以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１０１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１１４とを備えたＨＭＭの推定装置１００において、推定手段１１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段１２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ-(-はμの上に付く），共分散行列Σ-(-はΣの上に付く），初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび分岐の期待値Ｎ_ｃを更新設定する更新設定手段１２４と、更新設定手段１２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数１９２，数１９３を参照）と共に、前記レジスタ１０１から読み出した観測データと、初期設定手段１２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１４２〜数１４５，数１６７〜数１６９，数１８８〜数１９１を参照）演算手段１２６と、演算手段１２６による計算の収束を判定し、計算が収束していなければ、演算手段１２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段１２６で計算した新たな各確率量を最終的な値として出力させる判定手段１２８とを備えている。

このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ-(-はμの上に付く），共分散行列Σ-(-はμの上に付く），初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび分岐の期待値Ｎ_ｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
前記数１６７〜数１７３の式で計算する構成となっている。

また、ここでの演算手段１２６は、l+1回目の状態ｉとなる初期状態確率π_ｉ｜θl+1（l+1はθの添え字）の値を、前記数１４２の式で計算し、l+1回目の状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θl+1（t+1はθの添え字）の値を、l回目およびl-1回目における状態遷移の期待値Ｎ_ａij｜_θlおよびＮ_ａij｜_θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、前記数１４３〜数１４５の式で計算し、状態ｊで状態ｋが出力される出力確率ｂ_jk｜θl+1（l+1はθの添え字）の際に分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃ_jk｜θl+1の値と、平均値ベクトルμ-_k｜θl+1（l+1はθの添え字であり、-はμの上に付く）の値と、共分散行列Σ-_k｜θl+1(l+1はθの添え字であり、-はΣの上に付く）の値とを、l回目およびl-1回目における前記分岐の期待値Ｎ_cjk｜_θlおよびＮ_cjk｜_θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、前記数１６７〜数１６９、数１８８〜数１９１の式で計算する構成となっている。

そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。

さらに、ここでの更新設定手段１２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。

なお上述した本実施例の作用効果は、初期設定手段１２２としての動作を実行する初期設定ステップと、更新設定手段１２４としての動作を実行する更新設定ステップと、演算手段１２６としての動作を実行する演算ステップと、判定手段１２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。

上記実施例１〜６では、alpha-HMMの６つのタイプ、すなわち、{離散的,連続的,半連続的}×{単一配列,複数配列}について述べた。本実施例は、区切りのある連続アルファベットの場合である。これは、離散的な文字と連続的な文字が混在している場合の解釈である。この離散と連続が混在している更新方程式は、数１５４、数１６９にある総和を２段の総和に変形し、ゼロ出力の許容によって得られる。例えば、数１５４の分母の総和は、次式のような２段の総和に分割される。

ここで、D_gは、部分が重複しない集合である。Ｇは、その部分の濃度である。もし、D_gが連続的なアルファベットの特定のサブクラスに相当するならば、このことは、そのアルファベットが離散的なシンボルを伴うものと見なすことができる。換言すると、データが連続している各区間において、それぞれ総和を取り、さらに当該総和どうしについて総和を取っている。この離散的で連続的なアルファベットの場合にも、上記実施例のように単一配列の場合と複数配列の場合がある。したがって、この申請は、全部で８つのタイプのalpha-HMMを有している。

なお、上式の（・）という記号は、Σの中を省略し、総和の範囲を示すΣだけについて考察するために用いたものである。

この場合には、離散値は連続値を出す集団をグループ化した際のラベルと考えればよい。すなわち、データ系列は次式のように示される。

ここで、ラベルは、次式として考える。

上式のＧは、Ｇ≦Ｋ（分岐枝の数）である。そして、新たにｙ_t-（-はｙ_tの上に付く）をｙ_tと書き直せば、Σ^T _t=1（・）（TはΣの直上に付き、ｔ＝１は直下に付く）がグループ化に対応する部分を含んでいる場合には，これをΣ_t∈Dg（・）（t∈D_gがΣの直下に付く）に変更するだけでよい．ただし、D_g＝｛t|c_t＝ｇ｝である。

以下、各パラメータの更新式について述べる。
・離散連続混在系列の単数配列の場合
[単一配列の離散連続混在系列のalpha-HMMの初期状態確率]
初期状態確率は数１３０と同じである。
[単一配列の離散連続混在系列のalpha-HMMの状態遷移確率]
状態遷移確率は数１２８と同じである。
[単一配列の離散連続混在系列のalpha-HMMの分岐確率]
分岐確率はグループ化情報を併せもち、数１５４に代わって次式となる。以後、これを単にグループ化確率とよぶ。

[単一配列の離散連続混在系列のalpha-HMMの平均値ベクトル]
各グループに対する平均値ベクトルは数１５６に代わって次式となる。

[単一配列の離散連続混在系列のalpha-HMMの共分散行列]
各グループに対する共分散行列は、数１５８〜数１６０に代わって、次式となる。

各要素は次式である。

・離散連続混在系列の複数配列の場合
[複数配列の離散連続混在系列のalpha-HMMの初期状態確率]
初期状態確率は数１４２と同じである。
[複数配列の離散連続混在系列のalpha-HMMの状態遷移確率]
状態遷移確率は数１４３〜数１４５と同じである。
[複数配列の離散連続混在系列のalpha-HMMの分岐確率]
分岐確率すなわちグループ化確率は、数１６７〜数１６９の変形として次式のようになる。

各要素は次式である。

[複数配列の離散連続混在系列のalpha-HMMの平均値ベクトル]
平均値ベクトルは、数１７０、数１７１の変形として次式となる。

各要素は次式である。

[複数配列の離散連続混在系列のalpha-HMMの共分散行列]
共分散行列は、数１７２、数１７３の変形として次式となる。

各要素は次式である。

本実施例も上記第３〜６実施例と同様に、符号や未知パラメータの更新式が異なるだけであり、それらの数式を、単一配列の場合には、図２に示すハードウェア構成と図８に示すフローチャートに適用し、複数配列の場合には、図５に示すハードウェア構成と図９に示すフローチャートに適用するものであるため、ここでは、詳細な説明は省略する。なお、繰返し更新される未知パラメータの集合の組み合わせは数１６１とし、各更新式は、上述したように単一配列の場合は、数１３０、数１２８、数１９９〜２０３であり、複数配列の場合は、数１４２〜数１４５、数２０４〜数２１０である。

なお、本発明は上記実施例に限定されるものではなく、本発明の趣旨を逸脱しない範囲で変更可能である。例えば、上記各実施例の近似計算において、第３実施例で示した変形例のように遡る過去値を増やしたり、高次の近似にしたりすることも可能である。

以上のように本実施例（単一配列の場合もほとんど同様であるため、ここでは省略する）においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１０１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１１４とを備えたＨＭＭの推定装置１００において、推定手段１１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段１２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび分岐の期待値Ｎ_ｃを更新設定する更新設定手段１２４と、更新設定手段１２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記レジスタ１０１から読み出した観測データと、初期設定手段１２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１４２〜数１４５，数２０４〜数２１０を参照）演算手段１２６と、演算手段１２６による計算の収束を判定し、計算が収束していなければ、演算手段１２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段１２６で計算した新たな各確率量を最終的な値として出力させる判定手段１２８とを備えている。

このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎ_ａおよび分岐の期待値Ｎ_ｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。

また、ここでの演算手段１２６は、l+1回目の状態ｉとなる初期状態確率π_ｉ｜θl+1（l+1はθの添え字）の値を、前記数１４２の式で計算し、l+1回目の状態ｉから状態ｊに移る状態遷移確率ａ_ij｜θl+1（t+1はθの添え字）の値を、l回目およびl-1回目における状態遷移の期待値Ｎ_ａij｜_θlおよびＮ_ａij｜_θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、前記数１４３〜数１４５の式で計算し、状態jでのグループ化確率ｃ_jg｜θl+1（l+1はθの添え字）と、そのグループ化を反映した出力確率ｂ_jg｜θl+1の値と、平均値ベクトルμ_jg｜θl+1（l+1はθの添え字）の値と、共分散行列Σ_jg｜θl+1(l+1はθの添え字）の値とを、l回目およびl-1回目における前記分岐の期待値Ｎ_cjg｜_θlおよびＮ_cjg｜_θl-1（jgはｃの添え字、lまたはl-1はθの添え字）を利用して、前記数２０３〜数２０９の式で計算する構成となっている。

そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。

さらに、ここでの更新設定手段１２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。

なお上述した本実施例の作用効果は、初期設定手段１２２としての動作を実行する初期設定ステップと、更新設定手段１２４としての動作を実行する更新設定ステップと、演算手段１２６としての動作を実行する演算ステップと、判定手段１２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。

本発明で適用するalpha-ＨＭＭアルゴリズムは、音声認識と合成，ロボットの動作認識と生成，神経情報認識，および生命情報配列の認識などの非常に広範な応用性を有している。具体的には、コンピュータによる音声認識において、観測データとなる音声の特徴パラメータの時間的な変化と確率的な変動とを統計的に扱うために、ここで提案したalpha-ＨＭＭアルゴリズムを使用することができる。

その中で、alpha-ＨＭＭアルゴリズムを用いた音声認識における学習では、複数の音声を入力としてＨＭＭのパラメータを決定する。そして、実際の音声認識動作のときは、認識対象の音声よりパラメータが決定されたＨＭＭを用いてその確率を計算することが可能になる。また、音声認識等の場合は、離散的なデータに基づいて各確率量を算出するが、連続的なデータに基づいて各確率量を算出すれば、音声の合成や、ロボットの動作の生成が可能となる。

また、本発明はゲノム配列の認識（例えば、一乃至複数の癌患者からＤＮＡ情報を集めて、それらのデータを上記実施例で提案したalpha-ＨＭＭアルゴリズムに取り込んで癌になる確率を計算する等）、或いはロボットの動作認識などにも適用することができる。

Claims (23)

1. 観測データを時系列に格納する記憶手段と、
   前記観測データがどのような確率モデルであるのかを、隠れマルコフモデルの未知パラメータを算出することで推定する推定手段とを備えた隠れマルコフモデルの推定装置において、
   前記推定手段は、前記隠れマルコフモデルの高速化パラメータを設定する初期設定手段と、
   前記隠れマルコフモデルの未知パラメータとして、状態遷移する状態遷移確率ａ，出力が発生する出力確率ｂ，初期状態の初期状態確率πおよび尤度の各確率量と、
   状態遷移および出力の各期待値を更新設定する更新設定手段と、
   前記更新設定手段で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、
   前記記憶手段から読み出した観測データと、前記初期設定手段で設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算手段と、
   前記演算手段による計算の収束を判定し、
   計算が収束していなければ、前記演算手段で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、
   計算が収束していれば、前記演算手段で計算した新たな各確率量を最終的な値として出力させる判定手段とを備えた隠れマルコフモデルの推定装置。
2. t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂの集合の組合せを次の式で示し、
   
   ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
   前記演算手段は、前記新たな各確率量として、
   ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
   
   ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
   
   状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
   
   ことを特徴とする請求項１記載の隠れマルコフモデルの推定装置。
3. 前記更新設定手段は、前記高速化パラメータの値を１＜β＜３に設定することを特徴とする請求項２記載の隠れマルコフモデルの推定装置。
4. 記憶手段に観測データを時系列に格納し、
   前記観測データがどのような確率モデルであるのかを、推定手段が隠れマルコフモデルの未知パラメータを算出することで推定する隠れマルコフモデルの推定方法において、
   前記未知パラメータの算出は、前記隠れマルコフモデルの高速化パラメータを設定する初期設定ステップと、
   前記隠れマルコフモデルの未知パラメータとして、状態遷移する状態遷移確率ａ，出力が発生する出力確率ｂ，初期状態の初期状態確率πおよび尤度の各確率量と、
   状態遷移および出力の各期待値を更新設定する更新設定ステップと、
   前記更新設定ステップで更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定ステップで設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算ステップと、
   前記演算ステップによる計算の収束を判定し、計算が収束していなければ、前記演算ステップで計算した新たな各確率量および各期待値を前記更新設定ステップで設定更新させ、計算が収束していれば、前記演算ステップで計算した新たな各確率量を最終的な値として出力させる判定ステップとからなる隠れマルコフモデルの推定方法。
5. t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂの集合の組合せを次の式で示し、
   
   ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
   前記演算ステップは、前記新たな各確率量として、
   ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
   
   ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
   
   状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
   
   ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
6. t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂの集合の組合せを次の式で示し、
   
   ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
   前記演算ステップは、前記新たな各確率量として、
   ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
   
   ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
   
   状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
   
   ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
7. 前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、ｌ回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
   
   ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
   前記演算ステップは、前記新たな各確率量として、
   l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
   
   l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
   
   l+1回目の状態ｊで分岐ｋに移行する枝へ遷移する確率を特定する分岐確率cjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎcjk｜θlおよびＮcjk｜θl-1（jkはcの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
   
   前記分岐確率の確率密度関数の平均値である平均値ベクトルμjk｜θl+1（l+1はθの添え字）を次の式で計算し、
   
   l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σjk|θl+1（l+1はθの添え字）を次の式で計算する
   
   
   
   ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
8. 前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、ｌ回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
   
   ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
   前記演算ステップは、前記新たな各確率量として、
   l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
   
   l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
   
   
   
   l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率cjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎcjk｜θlおよびＮcjk｜θl-1（jkはcの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
   
   
   
   前記分岐確率の確率密度関数の平均値である平均値ベクトルμjk｜θl+1（l+1はθの添え字）を次の式で計算し、
   
   
   l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σjk|θl+1（l+1はθの添え字）を次の式で計算する
   
   
   ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
9. 前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、ｌ回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
   
   ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
   前記演算ステップは、前記新たな各確率量として、
   l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
   
   l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
   
   l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率cjk｜θl+（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎcjk｜θlおよびＮcjk｜θl-1（jkはcの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
   
   前記分岐確率の確率密度関数の平均値である平均値ベクトルμj｜θl+1（l+1はθの添え字）を次の式で計算し、
   
   l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σｊ|θl+1（l+1はθの添え字）を次の式で計算する
   
   
   
   ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
10. 前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、ｌ回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    
    ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    前記演算ステップは、前記新たな各確率量として、
    l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    
    l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    
    
    l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θl およびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    
    
    前記分岐確率の確率密度関数の平均値である平均値ベクトルμ-k｜θl+1（-はμの直上に付くものであり、l+1はθの添え字）を次の式で計算し、
    
    
    l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σ-ｋ|θl+1（-はΣの直上に付くものであり、l+1はθの添え字）を次の式で計算する
    
    
    ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
11. 前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、ｌ回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    
    ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    前記演算ステップは、前記新たな各確率量として、
    l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    
    l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    Dgを部分が重複しない集合として、l+1回目の状態jでグループ化gが行なわれるグループ化確率ｃjｇ｜θl+1（l+1はθの添え字）の値を次の式で計算し、
    
    前記分岐確率の確率密度関数の平均値である平均値ベクトルμjｇ｜θl+1（l+1はθの添え字）を次の式で計算し、
    
    l+1回目の状態jでグループ化gが行なわれるときの出力の共分散行列Σｊｇ|θl+1（l+1はθの添え字）を次の式で計算する
    
    
    
    ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
12. 前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、ｌ回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    
    ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    前記演算ステップは、前記新たな各確率量として、
    l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    
    l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    
    
    Dを部分が重複しない集合として、l+1回目の状態jでグループ化gが行なわれる確率の際に(g) g番目のグループへ遷移する確率を特定するグループ化確率ｃjｇ｜θl+1（l+1はθの添え字）の値を次の式で計算し、
    
    
    
    前記分岐確率の確率密度関数の平均値である平均値ベクトルμjｇ｜θl+1（l+1はθの添え字）を次の式で計算し、
    
    
    l+1回目の状態jでグループ化gが行なわれるときの出力の共分散行列Σｊｇ|θl+1（l+1はθの添え字）を次の式で計算する
    
    
    ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
13. 前記更新設定ステップで、前記高速化パラメータの値を１＜β＜３に設定することを特徴とする請求項５〜１２記載の隠れマルコフモデルの推定方法。
14. 記憶手段に時系列に格納された観測データがどのような確率モデルであるのかを、隠れマルコフモデルの未知パラメータを算出することで推定する推定手段として、コンピュータを機能させる隠れマルコフモデルの推定プログラムにおいて、
    前記推定手段を、前記隠れマルコフモデルの高速化パラメータを設定する初期設定手段と、
    前記隠れマルコフモデルの未知パラメータとして、状態遷移する状態遷移確率ａ，出力が発生する出力確率ｂ，初期状態の初期状態確率πおよび尤度の各確率量と、
    状態遷移および出力の各期待値を更新設定する更新設定手段と、
    前記更新設定手段で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定手段で設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算手段と、
    前記演算手段による計算の収束を判定し、計算が収束していなければ、前記演算手段による新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、前記演算手段による新たな各確率量を最終的な値として確定させる判定手段として機能させることを特徴とする隠れマルコフモデルの推定プログラム。
15. t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂの集合の組合せを次の式で示し、
    
    ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    前記演算手段は、前記新たな各確率量として、
    ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
    
    ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
    
    状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
    
    ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
16. t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂの集合の組合せを次の式で示し、
    
    ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    前記演算手段は、前記新たな各確率量として、
    ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
    
    ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
    
    状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
    
    ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
17. 前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、ｌ回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    
    ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    前記演算手段は、前記新たな各確率量として、
    l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    
    l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    前記分岐確率の確率密度関数の平均値である平均値ベクトルμjk｜θl+1（l+1はθの添え字）を次の式で計算し、
    
    l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σｊｋ|θl+1（l+1はθの添え字）を次の式で計算する
    
    
    
    ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
18. 前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、ｌ回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    
    ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    前記演算手段は、前記新たな各確率量として、
    l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    
    l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    
    
    l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    
    
    前記分岐確率の確率密度関数の平均値である平均値ベクトルμjk｜θl+1（l+1はθの添え字）を次の式で計算し、
    
    
    l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σｊｋ|θl+1（l+1はθの添え字）を次の式で計算する
    
    
    ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
19. 前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、ｌ回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    
    ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    前記演算手段は、前記新たな各確率量として、
    l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    
    l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    前記分岐確率の確率密度関数の平均値である平均値ベクトルμj｜θl+1（l+1はθの添え字）を次の式で計算し、
    
    l+1回目の状態jで分岐ｋに移行した時の出力確率の共分散行列Σｊ|θl+1（l+1はθの添え字）を次の式で計算する
    
    
    
    ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
20. 前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、ｌ回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    
    ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    前記演算手段は、前記新たな各確率量として、
    l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    
    l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    
    
    l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    
    
    前記分岐確率の確率密度関数の平均値である平均値ベクトルμ-k｜θl+1（-はμの直上に付くものであり、l+1はθの添え字）を次の式で計算し、
    
    
    l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率の共分散行列Σ-ｋ|θl+1（-はΣの直上に付くものであり、l+1はθの添え字）を次の式で計算する
    
    
    ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
21. 前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、ｌ回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    
    ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    前記演算手段は、前記新たな各確率量として、
    l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    
    l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    Dgを部分が重複しない集合として、l+1回目の状態jでグループ化gが行なわれるグループ化確率ｃjｇ｜θl+1（l+1はθの添え字）の値を次の式で計算し、
    
    前記グループ化確率の確率密度関数の平均値である平均値ベクトルμjｇ｜θl+1（l+1はθの添え字）を次の式で計算し、
    
    l+1回目の状態jでグループ化gが行なわれるときの出力の共分散行列Σｊｇ|θl+1（l+1はθの添え字）を次の式で計算する
    
    
    
    ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
22. 前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、ｌ回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    
    ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    前記演算手段は、前記新たな各確率量として、
    l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    
    l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    
    
    
    Dgを部分が重複しない集合として、l+1回目の状態jでグループ化gが行なわれる確率の際にg番目のグループへ遷移する確率を特定するグループ化確率ｃjｇ｜θl+1（l+1はθの添え字）の値を次の式で計算し、
    
    
    
    前記グループ化確率の確率密度関数の平均値である平均値ベクトルμjｇ｜θl+1（l+1はθの添え字）を次の式で計算し、
    
    
    l+1回目の状態jでグループ化gが行なわれるときの出力の共分散行列Σｊｇ|θl+1（l+1はθの添え字）を次の式で計算する
    
    
    ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
23. 前記更新設定手段は、前記高速化パラメータの値を１＜β＜３に設定することを特徴とする請求項１５〜２２記載の隠れマルコフモデルの推定プログラム。