(57)【要約】  特殊な条件下以外であっても未知パラメータの計算が可能な隠れマルコフモデル推定アルゴリズムを得る。  推定手段１４は、隠れマルコフモデルの未知パラメータとして、状態遷移確率ａ，出力確率ｂ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび出力の期待値Ｎｂを更新設定する更新設定手段２４と、更新設定手段２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記レジスタ１から読み出した観測データと、初期設定手段２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する演算手段２６とを備える。
2012008184.tif 000002
【特許請求の範囲】
【請求項１】
観測データを時系列に格納する記憶手段と、
前記観測データがどのような確率モデルであるのかを、隠れマルコフモデルの未知パラメータを算出することで推定する推定手段とを備えた隠れマルコフモデルの推定装置において、
前記推定手段は、前記隠れマルコフモデルの高速化パラメータを設定する初期設定手段と、
前記隠れマルコフモデルの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を更新設定する更新設定手段と、
前記更新設定手段で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定手段で設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算手段と、
前記演算手段による計算の収束を判定し、計算が収束していなければ、前記演算手段で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、前記演算手段で計算した新たな各確率量を最終的な値として出力させる判定手段とを備えた隠れマルコフモデルの推定装置。
【請求項２】
t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂとの集合の組合せを次の式で示し、
【数１】
2012008184.tif 000214
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
前記演算手段は、前記新たな各確率量として、
ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
【数２】
2012008184.tif 000215
ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
【数３】
2012008184.tif 000216
状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
【数４】
2012008184.tif 000217
ことを特徴とする請求項１記載の隠れマルコフモデルの推定装置。
【請求項３】
前記更新設定手段は、前記高速化パラメータの値を１＜β＜３に設定することを特徴とする請求項２記載の隠れマルコフモデルの推定装置。
【請求項４】
記憶手段に観測データを時系列に格納し、
前記観測データがどのような確率モデルであるのかを、推定手段が隠れマルコフモデルの未知パラメータを算出することで推定する隠れマルコフモデルの推定方法において、
前記未知パラメータの算出は、前記隠れマルコフモデルの高速化パラメータを設定する初期設定ステップと、
前記隠れマルコフモデルの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を更新設定する更新設定ステップと、
前記更新設定ステップで更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定ステップで設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算ステップと、
前記演算ステップによる計算の収束を判定し、計算が収束していなければ、前記演算ステップで計算した新たな各確率量および各期待値を前記更新設定ステップで設定更新させ、計算が収束していれば、前記演算ステップで計算した新たな各確率量を最終的な値として出力させる判定ステップとからなる隠れマルコフモデルの推定方法。
【請求項５】
t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂとの集合の組合せを次の式で示し、
【数５】
2012008184.tif 000218
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
前記演算ステップは、前記新たな各確率量として、
ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
【数６】
2012008184.tif 000219
ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
【数７】
2012008184.tif 000220
状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
【数８】
2012008184.tif 000221
ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
【請求項６】
t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂとの集合の組合せを次の式で示し、
【数９】
2012008184.tif 000222
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
前記演算ステップは、前記新たな各確率量として、
ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
【数１０】
2012008184.tif 000223
ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
【数１１】
2012008184.tif 000224
状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
【数１２】
2012008184.tif 000225
ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
【請求項７】
前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σとの集合の組合せを次の式で示し、
【数１３】
2012008184.tif 000226
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
前記演算ステップは、前記新たな各確率量として、
l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
【数１４】
2012008184.tif 000227
l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数１５】
2012008184.tif 000228
l+1回目の状態ｊで分岐ｋに移行する枝へ遷移する確率を特定する分岐確率cjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎcjk｜θlおよびＮcjk｜θl-1（jkはcの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数１６】
2012008184.tif 000229
前記分岐確率の確率密度関数の平均値である平均値ベクトルμjk｜θl+1（l+1はθの添え字）を次の式で計算し、
【数１７】
2012008184.tif 000230
l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σjk|θl+1（l+1はθの添え字）を次の式で計算する
【数１８】
2012008184.tif 000231
【数１９】
2012008184.tif 000232
【数２０】
2012008184.tif 000233
ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
【請求項８】
前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
【数２１】
2012008184.tif 000234
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
前記演算ステップは、前記新たな各確率量として、
l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
【数２２】
2012008184.tif 000235
l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数２３】
2012008184.tif 000236
【数２４】
2012008184.tif 000237
【数２５】
2012008184.tif 000238
l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率cjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎcjk｜θlおよびＮcjk｜θl-1（jkはcの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数２６】
2012008184.tif 000239
【数２７】
2012008184.tif 000240
【数２８】
2012008184.tif 000241
前記分岐確率の確率密度関数の平均値である平均値ベクトルμjk｜θl+1（l+1はθの添え字）を次の式で計算し、
【数２９】
2012008184.tif 000242
【数３０】
2012008184.tif 000243
l+1回目の状態jで状態jで分岐ｋに移行した時の出力の共分散行列Σjk|θl+1（l+1はθの添え字）を次の式で計算する
【数３１】
2012008184.tif 000244
【数３２】
2012008184.tif 000245
ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
【請求項９】
前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
【数３３】
2012008184.tif 000246
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
前記演算ステップは、前記新たな各確率量として、
l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
【数３４】
2012008184.tif 000247
l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数３５】
2012008184.tif 000248
l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率cjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎcjk｜θlおよびＮcjk｜θl-1（jkはcの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数３６】
2012008184.tif 000249
前記分岐確率の確率密度関数の平均値である平均値ベクトルμj｜θl+1（l+1はθの添え字）を次の式で計算し、
【数３７】
2012008184.tif 000250
l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σｊ|θl+1（l+1はθの添え字）を次の式で計算する
【数３８】
2012008184.tif 000251
【数３９】
2012008184.tif 000252
【数４０】
2012008184.tif 000253
ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
【請求項１０】
前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
【数４１】
2012008184.tif 000254
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
前記演算ステップは、前記新たな各確率量として、
l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
【数４２】
2012008184.tif 000255
l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数４３】
2012008184.tif 000256
【数４４】
2012008184.tif 000257
【数４５】
2012008184.tif 000258
l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数４６】
2012008184.tif 000259
【数４７】
2012008184.tif 000260
【数４８】
2012008184.tif 000261
前記分岐確率の確率密度関数の平均値である平均値ベクトルμ-k｜θl+1（-はμの直上に付くものであり、l+1はθの添え字）を次の式で計算し、
【数４９】
2012008184.tif 000262
【数５０】
2012008184.tif 000263
l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σ-ｋ|θl+1（-はΣの直上に付くものであり、l+1はθの添え字）を次の式で計算する
【数５１】
2012008184.tif 000264
【数５２】
2012008184.tif 000265
ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
【請求項１１】
前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
【数５３】
2012008184.tif 000266
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
前記演算ステップは、前記新たな各確率量として、
l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
【数５４】
2012008184.tif 000267
l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数５５】
2012008184.tif 000268
Dgを部分が重複しない集合として、l+1回目の状態jでグループ化gが行なわれるグループ化確率ｃjｇ｜θl+1（l+1はθの添え字）の値を次の式で計算し、
【数５６】
2012008184.tif 000269
前記分岐確率の確率密度関数の平均値である平均値ベクトルμjｇ｜θl+1（l+1はθの添え字）を次の式で計算し、
【数５７】
2012008184.tif 000270
l+1回目の状態jでグループ化gが行なわれるときの出力の共分散行列Σｊｇ|θl+1（l+1はθの添え字）を次の式で計算する
【数５８】
2012008184.tif 000271
【数５９】
2012008184.tif 000272
【数６０】
2012008184.tif 000273
ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
【請求項１２】
前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
【数６１】
2012008184.tif 000274
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
前記演算ステップは、前記新たな各確率量として、
l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
【数６２】
2012008184.tif 000275
l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数６３】
2012008184.tif 000276
【数６４】
2012008184.tif 000277
【数６５】
2012008184.tif 000278
Dgを部分が重複しない集合として、l+1回目の状態jでグループ化gが行なわれる確率の際にg番目のグループへ遷移する確率を特定するグループ化確率ｃjｇ｜θl+1（l+1はθの添え字）の値を次の式で計算し、
【数６６】
2012008184.tif 000279
【数６７】
2012008184.tif 000280
【数６８】
2012008184.tif 000281
前記分岐確率の確率密度関数の平均値である平均値ベクトルμjｇ｜θl+1（l+1はθの添え字）を次の式で計算し、
【数６９】
2012008184.tif 000282
【数７０】
2012008184.tif 000283
l+1回目の状態jでグループ化gが行なわれるときの出力の共分散行列Σｊｇ|θl+1（l+1はθの添え字）を次の式で計算する
【数７１】
2012008184.tif 000284
【数７２】
2012008184.tif 000285
ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
【請求項１３】
前記更新設定ステップで、前記高速化パラメータの値を１＜β＜３に設定することを特徴とする請求項５〜１２記載の隠れマルコフモデルの推定方法。
【請求項１４】
記憶手段に時系列に格納された観測データがどのような確率モデルであるのかを、隠れマルコフモデルの未知パラメータを算出することで推定する推定手段として、コンピュータを機能させる隠れマルコフモデルの推定プログラムにおいて、
前記推定手段を、前記隠れマルコフモデルの高速化パラメータを設定する初期設定手段と、
前記隠れマルコフモデルの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を更新設定する更新設定手段と、
前記更新設定手段で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定手段で設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算手段と、
前記演算手段による計算の収束を判定し、計算が収束していなければ、前記演算手段による新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、前記演算手段による新たな各確率量を最終的な値として確定させる判定手段として機能させることを特徴とする隠れマルコフモデルの推定プログラム。
【請求項１５】
t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂとの集合の組合せを次の式で示し、
【数７３】
2012008184.tif 000286
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
前記演算手段は、前記新たな各確率量として、
ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
【数７４】
2012008184.tif 000287
ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
【数７５】
2012008184.tif 000288
状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
【数７６】
2012008184.tif 000289
ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
【請求項１６】
t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂとの集合の組合せを次の式で示し、
【数７７】
2012008184.tif 000290
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
前記演算手段は、前記新たな各確率量として、
ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
【数７８】
2012008184.tif 000291
ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
【数７９】
2012008184.tif 000292
状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
【数８０】
2012008184.tif 000293
ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
【請求項１７】
前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σとの集合の組合せを次の式で示し、
【数８１】
2012008184.tif 000294
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
前記演算手段は、前記新たな各確率量として、
l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
【数８２】
2012008184.tif 000295
l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数８３】
2012008184.tif 000296
l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数８４】
2012008184.tif 000297
前記分岐確率の確率密度関数の平均値である平均値ベクトルμjk｜θl+1（l+1はθの添え字）を次の式で計算し、
【数８５】
2012008184.tif 000298
l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σｊｋ|θl+1（l+1はθの添え字）を次の式で計算する
【数８６】
2012008184.tif 000299
【数８７】
2012008184.tif 000300
【数８８】
2012008184.tif 000301
ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
【請求項１８】
前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
【数８９】
2012008184.tif 000302
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
前記演算手段は、前記新たな各確率量として、
l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
【数９０】
2012008184.tif 000303
l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数９１】
2012008184.tif 000304
【数９２】
2012008184.tif 000305
【数９３】
2012008184.tif 000306
l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数９４】
2012008184.tif 000307
【数９５】
2012008184.tif 000308
【数９６】
2012008184.tif 000309
前記分岐確率の確率密度関数の平均値である平均値ベクトルμjk｜θl+1（l+1はθの添え字）を次の式で計算し、
【数９７】
2012008184.tif 000310
【数９８】
2012008184.tif 000311
l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σｊｋ|θl+1（l+1はθの添え字）を次の式で計算する
【数９９】
2012008184.tif 000312
【数１００】
2012008184.tif 000313
ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
【請求項１９】
前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
【数１０１】
2012008184.tif 000314
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
前記演算手段は、前記新たな各確率量として、
l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
【数１０２】
2012008184.tif 000315
l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数１０３】
2012008184.tif 000316
l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数１０４】
2012008184.tif 000317
前記分岐確率の確率密度関数の平均値である平均値ベクトルμj｜θl+1（l+1はθの添え字）を次の式で計算し、
【数１０５】
2012008184.tif 000318
l+1回目の状態jで分岐ｋに移行した時の出力確率の共分散行列Σｊ|θl+1（l+1はθの添え字）を次の式で計算する
【数１０６】
2012008184.tif 000319
【数１０７】
2012008184.tif 000320
【数１０８】
2012008184.tif 000321
ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
【請求項２０】
前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
【数１０９】
2012008184.tif 000322
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
前記演算手段は、前記新たな各確率量として、
l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
【数１１０】
2012008184.tif 000323
l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数１１１】
2012008184.tif 000324
【数１１２】
2012008184.tif 000325
【数１１３】
2012008184.tif 000326
l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数１１４】
2012008184.tif 000327
【数１１５】
2012008184.tif 000328
【数１１６】
2012008184.tif 000329
前記分岐確率の確率密度関数の平均値である平均値ベクトルμ-k｜θl+1（-はμの直上に付くものであり、l+1はθの添え字）を次の式で計算し、
【数１１７】
2012008184.tif 000330
【数１１８】
2012008184.tif 000331
l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率の共分散行列Σ-ｋ|θl+1（-はΣの直上に付くものであり、l+1はθの添え字）を次の式で計算する
【数１１９】
2012008184.tif 000332
【数１２０】
2012008184.tif 000333
ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
【請求項２１】
前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
【数１２１】
2012008184.tif 000334
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
前記演算手段は、前記新たな各確率量として、
l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
【数１２２】
2012008184.tif 000335
l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数１２３】
2012008184.tif 000336
Dgを部分が重複しない集合として、l+1回目の状態jでグループ化gが行なわれるグループ化確率ｃjｇ｜θl+1（l+1はθの添え字）の値を次の式で計算し、
【数１２４】
2012008184.tif 000337
前記グループ化確率の確率密度関数の平均値である平均値ベクトルμjｇ｜θl+1（l+1はθの添え字）を次の式で計算し、
【数１２５】
2012008184.tif 000338
l+1回目の状態jでグループ化gが行なわれるときの出力の共分散行列Σｊｇ|θl+1（l+1はθの添え字）を次の式で計算する
【数１２６】
2012008184.tif 000339
【数１２７】
2012008184.tif 000340
【数１２８】
2012008184.tif 000341
ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
【請求項２２】
前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
【数１２９】
2012008184.tif 000342
ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
前記演算手段は、前記新たな各確率量として、
l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
【数１３０】
2012008184.tif 000343
l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
【数１３１】
2012008184.tif 000344
【数１３２】
2012008184.tif 000345
【数１３３】
2012008184.tif 000346
Dgを部分が重複しない集合として、l+1回目の状態jでグループ化gが行なわれる確率の際にg番目のグループへ遷移する確率を特定するグループ化確率ｃjｇ｜θl+1（l+1はθの添え字）の値を次の式で計算し、
【数１３４】
2012008184.tif 000347
【数１３５】
2012008184.tif 000348
【数１３６】
2012008184.tif 000349
前記グループ化確率の確率密度関数の平均値である平均値ベクトルμjｇ｜θl+1（l+1はθの添え字）を次の式で計算し、
【数１３７】
2012008184.tif 000350
【数１３８】
2012008184.tif 000351
l+1回目の状態jでグループ化gが行なわれるときの出力の共分散行列Σｊｇ|θl+1（l+1はθの添え字）を次の式で計算する
【数１３９】
2012008184.tif 000352
【数１４０】
2012008184.tif 000353
ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
【請求項２３】
前記更新設定手段は、前記高速化パラメータの値を１＜β＜３に設定することを特徴とする請求項１５〜２２記載の隠れマルコフモデルの推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、あるデータ系列が与えられたときに、それを生成する隠れマルコフモデル（ＨＭＭ：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の間の凸ダイバージェンスについて考える。
【００２８】
【数１】
2012008184.tif 000003
【００２９】
ここで、ＹはＫ次元ユークリッド空間である。関数ｆ(r）は、(０，∞）において凸である。関数f（r）の双対関数g（r）は、次式を満たす必要がある。
【００３０】
【数２】
2012008184.tif 000004
【００３１】
上記数１における「≧」という記号は、ｐ＝ｑを殆どどこでも保持する場合に限り保持する。ｆ（1）の正規化が任意であるため、ここでは、Ｄｆ（ｐ・ｑ）とＤｇ（ｑ・ｐ）がｐとｑとの間の擬似距離として作用するように、ｆ（1）を次のように選定する。
【００３２】
【数３】
2012008184.tif 000005
【００３３】
留意すべきなのは、次の２つの式が、カルバック・ライブラー・ダイバージェンス（Kullback-Leibler divergence）を発生させる、ということである。もしｐとｑが確率量の関数であれば、数１における積分は総和に置き換えられる。これは、以後説明する積分の全てにおいて同様である。
【００３４】
【数４】
2012008184.tif 000006
【００３５】
【数５】
2012008184.tif 000007
【００３６】
上記数１〜数５の関係は、対数より有能であるものの、対数と類似した関数であることを意味している。連続して２回微分可能な関数のクラスについて考えると、数２は以下の等式を与える。
【００３７】
【数６】
2012008184.tif 000008
【００３８】
また、ｒ＝１付近では以下の式がそれぞれ成り立つ。
【００３９】
【数７】
2012008184.tif 000009
【００４０】
【数８】
2012008184.tif 000010
【００４１】
ここで、o（1）は高次無限小の項を意味する。もし、乗算的に分離可能な凸関数を使用するのであれば、o（1）は不要であることに注目されたい。
【００４２】
【数９】
2012008184.tif 000011
【００４３】
以下の関数は、その一例である。
【００４４】
【数１０】
2012008184.tif 000012
【００４５】
もし、上記の関数を、ｆ（1）＝０およびｆ”（1）＝１であるようにシフトして正規化すれば、ｋ＝−ｃ（１−ｃ）となる。前記数７のコアとなる項は以下の式を有している。
【００４６】
【数１１】
2012008184.tif 000013
【００４７】
上式において、次の式は興味深い関数であり、これはパラメータｃによって凸性が調整できる単調関数である。
【００４８】
【数１２】
2012008184.tif 000014
【００４９】
これは、以下の式の関係が保たれることを指し示すのに重要である。
【００５０】
【数１３】
2012008184.tif 000015
【００５１】
したがって、Ｌ(c)（r）はパラメータｃを有する法則化された対数と見なされる。つまり、ｃ-対数と呼ぶことができる。しかしながら、パラメータｃの増加はその数値を減少させる。そこで、ここでは以下の関係によって、法則化された対数をパラメータ化する。
【００５２】
【数１４】
2012008184.tif 000016
【００５３】
このように、α-対数（alpha-log）は以下の式に示すような凸ダイバージェンスＤｆ（ｐ・ｑ）から始まることによって得られる。
【００５４】
【数１５】
2012008184.tif 000017
【００５５】
ここでの対数は次式のようになる。
【００５６】
【数１６】
2012008184.tif 000018
【００５７】
注目すべきは、次の式に示す凸関数が、ｒ=ｐ/ｑとして上記数１に適用されるとαダイバージェンスが発生し、これはalpha-ＥＭアルゴリズムにおいて重要な役割を果たす。
【００５８】
【数１７】
2012008184.tif 000019
【００５９】
αダイバージェンスは、次の式であらわせる。
【００６０】
【数１８】
2012008184.tif 000020
【００６１】
上式において、ｐ＝ｑが殆どどこでも適用される場合に限り、等式のゼロが達成される。
【００６２】
次に、alpha-ＨＭＭ再推定アルゴリズムの原点となるalpha-ＥＭアルゴリズムの期待値と最大化の処理ステップについて説明する。ここではＰＹ｜ψ（ｙ｜ψ）を、ψによってパラメータ化されたＹにおける観測データｙの確率密度または確率量とする。集合Ｙは、不完全データの集合と見なされる。ｘ∈Ｘ（太字）は、未知または消失したデータを含む理想的観測の完全データあるいは増補データとする。そして、不完全データのpdf（probability density function：確率密度関数）あるいはpmf（probability mass function：確率量関数）は次式に示される。
【００６３】
【数１９】
2012008184.tif 000021
【００６４】
ここでの積分範囲は、次式のようになり、pdfの積分はpmfの総和となる。
【００６５】
【数２０】
2012008184.tif 000022
【００６６】
そして、条件付きのpdfまたはpmfは次式の通りとなる。
【００６７】
【数２１】
2012008184.tif 000023
【００６８】
pdfまたはpmfの不完全なデータのalpha-log尤度比は、次式の通りである。
【００６９】
【数２２】
2012008184.tif 000024
【００７０】
ここで、次の式におけるαダイバージェンスの計算は、alpha-ＥＭアルゴリズムの基本的な関係を示している。
【００７１】
【数２３】
2012008184.tif 000025
【００７２】
その式は、以下の通りである。
【００７３】
【数２４】
2012008184.tif 000026
【００７４】
したがって、最大化のための目的関数は以下の式のような定量化（quantity）によって表される。
【００７５】
【数２５】
2012008184.tif 000027
【００７６】
上記説明によって、alpha-ＥＭアルゴリズムがalpha-ＨＭＭ再推定アルゴリズムの開始点と位置付けられることがわかる。
【００７７】
alpha-ＥＭアルゴリズムでは、次のステップで処理が行なわれる。
【００７８】
初期設定：最初のサイクル用のパラメータφを設定する。
【００７９】
Ｅ-ステップ：φが与えられたときのψの関数として、前記数２５を計算する。
【００８０】
Ｍ-ステップ：次の式にしたがって、更新パラメータを計算する。
【００８１】
【数２６】
2012008184.tif 000028
【００８２】
Ｕ-ステップ：φをψ*に置換して、収束するまではＥ-ステップに戻る。 
もし、閉じた部分で正確な最大化が得られなければ、数２６に示す最大化が改善した計算に置き換わる。
【００８３】
また、別なalpha-ＧＥＭアルゴリズムでは、次のステップで処理が行なわれる。
【００８４】
初期設定：最初のサイクル用のパラメータφを設定する。
【００８５】
Ｅ-ステップ：φが与えられたときのψの関数として、前記数２５を計算する。
【００８６】
Ｍ-ステップ：次の式にしたがって、更新パラメータを計算する。
【００８７】
【数２７】
2012008184.tif 000029
【００８８】
Ｕ-ステップ：φをψ＋に置換して、収束するまではＥ-ステップに戻る。alpha-ＨＭＭの近似されたバージョンが、ＧＥＭアルゴリズムと見なすことができる点に注目されたい。
【００８９】
隠れマルコフモデル（ＨＭＭ）の公式化の問題は、一連のランダムな観測変数データを発生させる最良のマルコフモデルを推定することにある。観測データの変数列は、次の式であらわせる。この式におけるtはマルコフモデルを生成するためのデータの順序を表すものであり、繰り返しのインデックスとは意味が異なる。
【００９０】
【数２８】
2012008184.tif 000030
【００９１】
ここでの「最良」とは、隠れマルコフモデルの尤度が最大化されることを意味する。隠れマルコフモデルを公式化するために、次に定義する各確率を推定する必要がある。
【００９２】
（ａ）状態遷移確率
状態がｉからｊに遷移する確率のことであり、次式のように定義される。
【００９３】
【数２９】
2012008184.tif 000031
【００９４】
上記式において、ｓi∈Σは有限な範囲での総和Σの中の１つの状態を示すアルファベットの記号である。状態ｓkから状態ｓlへのつながりが無い場合、状態遷移確率ａkl＝０となる。これは、状態遷移における従前のトポロジーを反映している。
【００９５】
（ｂ）出力確率
状態ｓjにおいて出力ｙkが発生する確率のことであり、次式のように定義される。
【００９６】
【数３０】
2012008184.tif 000032
【００９７】
ここで、上記式における各要素は以下の通りである。
【００９８】
【数３１】
2012008184.tif 000033
【００９９】
（ｃ）初期状態確率
最初の状態がｉである確率のことであり、次式のように定義される。
【０１００】
【数３２】
2012008184.tif 000034
【０１０１】
推定される上記（ａ）〜（ｃ）の各要素は、以下のように示される。
【０１０２】
【数３３】
2012008184.tif 000035
【０１０３】
つまり、ここで用いる記号θは、状態遷移確率，出力確率および初期状態確率を含ませたものである。
【０１０４】
一方、隠れマルコフモデル推定の問題に関し、前記数２８で示したランダムな観測データの変数列は、その値が次の式であらわされる。ここで、tはデータの順序を表すインデックスである。
【０１０５】
【数３４】
2012008184.tif 000036
【０１０６】
与えられた変数列は、その変数列を発生できる最大尤度モデルθを推定する。この隠れマルコフモデル推定の問題は、以下の解釈をもたらす。
（ａ）数２８で示す観測データが、前述した不完全データＹとして与えられる。
（ｂ）次の式に示す状態遷移は、消失データＺである。
【０１０７】
【数３５】
2012008184.tif 000037
【０１０８】
（ｃ）完全データは、次式のように示される。
【０１０９】
【数３６】
2012008184.tif 000038
【０１１０】
この確率は、次式で示される。
【０１１１】
【数３７】
2012008184.tif 000039
【０１１２】
（ｄ）前述の三重項θは、数３６に示す完全データの確率構造を提供する。
【０１１３】
この解釈は期待値最大化によって、実際にはalpha-ＥＭアルゴリズムによって、前記ＨＭＭ問題の公式化を可能にする。そして、以下の重要な定理が得られる。
【０１１４】
定理：ψとφを、２つのマルコフモデルの各パラメータセットとする。そして、α≦１に対して、次の数３８の式は、数３９という意味を持つ。
【０１１５】
【数３８】
2012008184.tif 000040
【０１１６】
【数３９】
2012008184.tif 000041
【０１１７】
ここで数３９の不等式は、数３８の不等式が成立する場合に成立するものである。
【０１１８】
証明：基本の方程式である数２４において、次式に示す不完全データのalpha-log尤度比を設定する。
【０１１９】
【数４０】
2012008184.tif 000042
【０１２０】
数３５に示す状態Ｓ（太字）が消失データであるとすると、Ｘ＝（Ｓ，Ｙ）（各記号は何れも太字）であるので、前記数２５のＱ-関数は次式で算出できる。
【０１２１】
【数４１】
2012008184.tif 000043
【０１２２】
前記数１８の特性と数２４の基本的な方程式によって、この定理が保たれる。
【０１２３】
数４１の特性は、数４０の形態で推定されるマルコフモデルψが、マルコフモデルφよりも高い確率で数列ｙ（太字）を生成することを意味する。すなわち、 数３３に示すマルコフモデルクラスでのalpha-ＥＭアルゴリズムの繰り返しは、必然的な局所最適性のコストで最良の隠れマルコフモデルを生成する。したがって、数３３に示すマルコフモデルのパラメータで数３８を繰り返すことが、ＨＭＭアルゴリズムの理論版となる。ここでは、α＝１の場合（後述する高速化パラメータの値が、β＝３である場合)は、収束限界であるために実用的ではない。
【０１２４】
alpha-ＥＭアルゴリズムとしての主な問題は、ソフトウェアとして実現できる具体的なコンピュータアルゴリズムをいかに提供するのかにある。数３８のalpha-ＨＭＭは理論上の形態であり、これはlog-ＥＭアルゴリズムと同様に一般的なＨＭＭである。ソフトウェアで実行可能なalpha-ＨＭＭとして、ここでは２つのバージョンを提示するが、第１のバージョンは正確であるものの、未来の情報の計算を必要とする。すなわち、得られたアルゴリズムはnon-causalである。一方、本発明で特に提案する第２のバージョンは、時間シフトと確率の近似を用いることによって、未来の情報を含めることなく計算が可能となる。なお、ここで提示するalpha-ＨＭＭアルゴリズムの生成において、全てはpmfであるので、積分が総和となる。以下、それぞれのバージョンについて説明する。
【０１２５】
（Ａ）第１の正確なバージョン
・non-causalな更新
正確なバージョンは、条件付き微分によって得られる。このバージョンで得られる更新がnon-causalであることは、既に述べている。この形態は次のバージョンのcausalな形態を得るための原点となるものである。
【０１２６】
次式に示すように、全ての更新式は、ｘ＝（ｓ，ｙ）（各記号は何れも太字）に対する前記Ｑ-関数の最大化から得られる。
【０１２７】
【数４２】
2012008184.tif 000044
【０１２８】
ここで、ｔはalpha-ＨＭＭに対する繰り返しの指標（インデックス）である。そして、パラメータ集合は次式に示すようになる。なお、明細書中、小文字のπ、a、b等の記号は集合として表記する場合は大文字で示すようにしている。
【０１２９】
【数４３】
2012008184.tif 000045
【０１３０】
・Πθt+1（Πは太字、t+1はθの添え字）の推定
初期状態確率πｉ｜θt+1（πは太字、t+1はθの添え字）の更新は、次の修正されたＱ-関数における微分の最大化によって得られる。
【０１３１】
【数４４】
2012008184.tif 000046
【０１３２】
この計算は、次式の微分を用いることによって達成できる。
【０１３３】
【数４５】
2012008184.tif 000047
【０１３４】
そして、次式を用いることで、数４４のλを除去する。
【０１３５】
【数４６】
2012008184.tif 000048
【０１３６】
すると、次の更新式を与える。
【０１３７】
【数４７】
2012008184.tif 000049
【０１３８】
上記数４７は観測データｙ（太字）に対して計算できるように見えるが、左辺に含まれるθt+1が右辺に含まれている。これは、non-causalityの存在による自己撞着性を示している。しかしながら、これは次のバージョンの近似によって解決されることとなる。したがって、確率のための２つの更なる更新式を得ることができた。
【０１３９】
・Ａθt+1（Ａは太字、t+1はθの添え字）の推定
状態遷移確率ａij｜θt+1（t+1はθの添え字）の更新は、同様に微分によって得られる。
【０１４０】
【数４８】
2012008184.tif 000050
【０１４１】
ここで次式を用いることで、λを除去する。
【０１４２】
【数４９】
2012008184.tif 000051
【０１４３】
すると、次の更新式を与える。
【０１４４】
【数５０】
2012008184.tif 000052
【０１４５】
ここでのＮij（ｓ）（ｓは太字）は、状態ｓにおいてｉからｊに移行する状態遷移の数である。数５０も右辺にθt+1を含んでいる。これも次のバージョンで解決されることとなる。
【０１４６】
・Ｂθt+1（Ｂは太字、t+1はθの添え字）の推定
出力確率ｂjk｜θt+1（t+1はθの添え字）の更新は、再度の微分によって得られる。
【０１４７】
【数５１】
2012008184.tif 000053
【０１４８】
ここで次式を用いることで、λを除去する。
【０１４９】
【数５２】
2012008184.tif 000054
【０１５０】
すると、次の更新式を与える。
【０１５１】
【数５３】
2012008184.tif 000055
【０１５２】
ここでのＮｂjk（ｓ）（ｓは太字、jkはｂの添え字）は、数列ｓ（太字）によって生じる出力確率ｂjkに対する事象の数である。数５３も右辺にθt+1を含んでいる。これは次のバージョンで解決されることとなる。
【０１５３】
（Ｂ）第２のCausalな近似バージョン
前述したように、第１のバージョンに示す一連の更新式（数４７，数５０，数５３）は、自己撞着性を解決するために近似を必要とする。
【０１５４】
・時間のシフト
alpha-ＥＭアルゴリズムの収束は、適切な収束判定基準によって、次の式のようになることを意味する。
【０１５５】
【数５４】
2012008184.tif 000056
【０１５６】
これは、以下の近似式を与える。ここで、o(1)は高位の無限小を意味する。
【０１５７】
【数５５】
2012008184.tif 000057
【０１５８】
これは、繰り返し指標tをシフトすることによって、数４７，数５０，数５３に示す更新式の自己撞着性を解決することになる。
【０１５９】
・αのシフト
前記時間のシフトは、尤度比の期待値のための確率環境が反転されることを意味する。したがって、元のパラメータαはαcausalに変換される。αとαcausalの関係は、以下の誘導によって得られる。ここで、数４７，数５０，数５３における尤度比の期待値の核となる部分を考慮すると、数４７の場合は以下の関係が成り立つ。
【０１６０】
【数５６】
2012008184.tif 000058
【０１６１】
この誘導は尤度比の時系列を維持しながら、期待値が現在の環境Ｐ（ｓ｜ｙ，θt）（ｓ，ｙは太字）から、未来の環境Ｐ（ｓ｜ｙ，θt+1）（ｓ，ｙは太字）に変化することを意味する。上記数５６を見るとθｔ＋１の初期確率πiを左辺として計算するのに、第一式の右辺では同じ時系列のθｔ＋１に依存する関数（項）をしなければならないという問題点があり、このままでは計算できないと言う自己撞着性を含んでいた。そこで、まず第二式では第一式の分母分子を｛｝内の項のα乗でそれぞれ除算する。さらに第三式ではテーラー展開によって微小近似すると共に時間シフトの概念を導入しており、これによって自己撞着矛盾を解決している。言い換えると、右辺には左辺と同時系列のθｔ＋１に依存する関数（項）を用いない形になっている。この詳細な考え方を以下に示す。
【０１６２】
ａ）未来に関する尤度比を現在の環境で計算すると、次式のようになる。
【０１６３】
【数５７】
2012008184.tif 000059
【０１６４】
ｂ）未来に関する尤度比を未来の環境で計算すると、次式のようになる。
【０１６５】
【数５８】
2012008184.tif 000060
【０１６６】
したがって、時間シフト前後の関係は以下の通りである。
【０１６７】
【数５９】
2012008184.tif 000061
【０１６８】
そこから得られる関係式は、次式のようになる。
【０１６９】
【数６０】
2012008184.tif 000062
【０１７０】
したがって、α＝−１のlog-ＥＭとlog-ＨＭＭは、ここでのαcausal＝１の場合に相当する。なお、符号の使用を簡単化するために、以後、高速化パラメータの値となるαcausalを、次式のようにβとしてあらわす。すると、数６０の右側の条件式はβ≦３を示すことになり、従来技術であるlog-ＥＭとlog-ＨＭＭはβ＝１の場合に対応することがわかる。（ただし、β＝３だと後述の実験結果より発散してしまうことがわかっている。）
【０１７１】
【数６１】
2012008184.tif 000063
【０１７２】
・テーラー展開
前記全ての更新式（数４７，数５０，数５３）は、尤度比の計算に能力を必要とする。これは計算に時間がかかるものであり、Ｔが増加するに従って扱いにくくなる。テーラー展開は、alpha-logアルゴリズムにおける尤度比の長所を失うことなく、計算を簡単化することができるもので、尤度比は次式のように近似化される。
【０１７３】
【数６２】
2012008184.tif 000064
【０１７４】
数５５と数６２を数４７，数５０，数５３に適用すると、過去の情報を利用するcausalで、かつ計算上効率的なalpha-ＨＭＭアルゴリズムが得られる。注目すべきは、これが数２７のalpha-ＧＥＭアルゴリズムに対応するということである。
【０１７５】
（Ｃ）ソフトウェアで実行可能なalpha-ＨＭＭアルゴリズム
・Πθt+1（Πは太字、t+1はθの添え字）の推定
高次の項ｏ（1）を放棄することによって数５５と数６２の近似式を結合し、数４７に適用すると以下の更新式が与えられる。
【０１７６】
【数６３】
2012008184.tif 000065
【０１７７】
・Ａθt+1（Ａは太字、t+1はθの添え字）の推定
数５５と数６２の近似式を適用すると、数５０は以下の更新式を与える。
【０１７８】
【数６４】
2012008184.tif 000066
【０１７９】
ここで、Ｎａij｜θt（ijはａの添え字、tはθの添え字）は、条件付き確率Ｐ（ｓ｜ｙ，θt）（ｓ，ｙは太字）の下でのＮij（ｓ）（ｓは太字）の期待値であり、不完全データｙ（ｙは太字）によって生じる状態遷移のカウント数により推定される。
【０１８０】
・Ｂθt+1（Ｂは太字、t+1はθの添え字）の推定
数５５と数６２の近似式を適用すると、数５３は以下の更新式を与える。
【０１８１】
【数６５】
2012008184.tif 000067
【０１８２】
ここで、ＮｂjK｜θt（jkはｂの添え字、tはθの添え字）は、条件付き確率Ｐ（ｓ｜ｙ，θt）（ｓ，ｙは太字）の下でのＮｂjK（ｓ）（ｓは太字）の期待値であり、不完全データｙ（ｙは太字）によって生じる状態遷移のカウント数によって推定される。
【０１８３】
・従来のlog-ＨＭＭとの比較
上記の数６３，数６４，数６５は、ソフトウェアとして実行可能である。これらの更新式は、以下に示す共通の特性を有する。
（１）β＝１すなわちα＝−１の場合は、log-ＥＭアルゴリズムから得られた従来のＨＭＭアルゴリズムと一致する。
（２）数６０と数６１との間の関係により、αが行なうのと同様に、パラメータβは更新用の現在と過去の情報の各要素を調整する。
（３）βに依存する、すなわちαに依存するすべての過去の更新値は、格納データを参照することによってのみ得ることができる。
（４）現在の項は、log-ＨＭＭアルゴリズムに対する前方−後方アルゴリズムを用いることで効率的に計算できる。
【０１８４】
上記特性（１）〜（４）は、alpha-ＨＭＭアルゴリズムの各更新サイクルがlog-ＨＭＭアルゴリズムのオーバーヘッドより僅かに増えてだけであることを示唆している。したがって、設計パラメータα、あるいは高速化パラメータβが適切に選択されれば、本発明のalpha-ＨＭＭアルゴリズムは、従来のlog-ＨＭＭアルゴリズムよりも高速に収束することが予測される。この違いが生じる理由は、数６３、数６４、数６５の式において、従来のlog-ＥＭから得られたＨＭＭアルゴリズムに相当するβ＝１を代入すると分母分子の第二項が消去され、分母分子の第一項が約分されて尤度比が残らなかったのに対し、β＞1で設定すると分母分子の第二項が消去されず、単純に分母分子を約分できないために尤度比の影響を残すことができるようになったことだと考えられる。
【０１８５】
以下、上記理論に基づく好ましい実施例について説明する。
【０１８６】
実施例１、２では、観測データが離散的な場合において、観測データが単一配列（一本の観測データ）を有する場合、複数配列（Ｍ本の観測データ）を有する場合について、それぞれ説明する。
【０１８７】
実施例３、４では、観測データが連続的な場合において、観測データが単一配列（一本の観測データ）を有する場合、複数配列（Ｍ本の観測データ）を有する場合について、それぞれ説明する。
【０１８８】
実施例５、６では、観測データが半連続的な場合において、観測データが単一配列（一本の観測データ）を有する場合、複数配列（Ｍ本の観測データ）を有する場合について、それぞれ説明する。
【０１８９】
実施例７では、観測データが離散的かつ連続的な場合、すなわち離散系列と連続系列が混在している場合において、単一配列を有する場合と複数配列を有する場合をまとめて説明する。
【実施例１】
【０１９０】
図１は、本発明の第１実施例において、上述したalpha-ＨＭＭ再推定アルゴリズムを実行可能にするプログラムの処理手順をあらわしたものである。【０１９１】
同図において、１は合計でＴ個の観測データを格納する記憶手段としてのレジスタで、各観測データは時刻τが１から順に時系列に並んで格納される。ここでｙτは個々のデータ値を示し、ｙはＴ個のデータ値の集合を示しており、本実施例では配列が１本の観測データをレジスタ１に格納している。推定装置１０は、レジスタ１に記憶される一列のデータ値が、どのような確率構造（モデル）を有しているのかを、ＨＭＭの未知パラメータを算出する以下のステップＳ１〜ステップＳ１０の各手順に従って推定解析するものである。【０１９２】
ステップＳ１は、前述した高速化パラメータの値βを設定する部分である。従来のＨＭＭ推定アルゴリズムは、β＝１の場合に相当する。つまり、β＝１という特殊な条件下であれば、従来のalpha-ＨＭＭ再推定アルゴリズムであっても、上記確率モデルの計算は可能である。しかし、本実施例で提案する新たなalpha-ＨＭＭ再推定アルゴリズムでは、確率モデルの計算が可能な高速化パラメータの値βを１≦β＜３の範囲に拡張することができる。これは、従来のβ＝１を特例として含むものである。なお、βが大きな値であるほど、推定装置１０としての処理の高速性は増すが、収束性を保持するにはβ＜３でなければならない。【０１９３】
ステップＳ２は、初期確率と収束判定値を決める部分である。これは、後述するステップＳ３〜ステップＳ８の手順を繰り返す前に行なわれる。推定装置１０が最終的に算出しようとするＨＭＭの確率構造は、次式のようにあらわせる。これは、前述の数３７と等しい。【０１９４】
【数６６】
2012008184.tif 000068
【０１９５】
上記式において、πS0（０はｓの添え字）は最初の状態ｓ０における確率であり、ａSτ-1Sτ（τ−１およびτはｓの添え字）は時刻τ−１の状態ｓτ-1が時刻τの状態ｓτに移る確率であり、ｂSτ（ｙτ）（τはｓの添え字）は時刻τの状態ｓτに移行したときにｙτが出力される確率である。推定装置１０は、レジスタ１に格納された観測データを読み出して、その観測データが最も出現しやすくなる初期状態確率πと、状態遷移確率ａと、出力確率ｂとを推定するが、このステップＳ２では、次式に示すように、それらの確率π，ａ，ｂの初期値を決定する。【０１９６】
【数６７】
2012008184.tif 000069
【０１９７】
上式において、θ０は数４３にあるように、０回目の繰り返しにおける初期の確率π，ａ，ｂの組合せを示し、推定装置１０は、その条件で状態ｉから状態ｊに移る状態遷移確率ａij｜θ０（０はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂjk｜θ0（０はθの添え字）の値と、最初に状態ｉとなる初期状態確率πｉ｜θ0（０はθの添え字）の値とをそれぞれ決定する。【０１９８】
またステップＳ２では、対数尤度に基づく収束範囲を決めることで、ステップＳ３〜ステップＳ９の繰り返しを終了させるための収束判定値を決定する。この収束判定値はステップＳ８で用いられ、具体的には後述する数８３で示される。【０１９９】
これらの値をステップＳ１，Ｓ２で決定すると、ステップＳ３の手順に移行して、確率量とカウント値の設定が行なわれる。ステップＳ３は、最初にステップＳ２で求めた初期確率を利用し、それ以降はステップＳ９で更新された確率量とカウント値を利用して、数６８に示す各確率値と、数６９に示す各カウント値を、実際に推定装置１０のメモリ（図示せず）に設定する部分である。【０２００】
【数６８】
2012008184.tif 000070
【０２０１】
【数６９】
2012008184.tif 000071
【０２０２】
次のステップＳ４は、レジスタ１から観測データを読み出して、ステップＳ３で設定された上式の確率量から、数７０と数７１に示す前向き確率を計算する部分である。なお、ここに示すαは確率値であり、前述したパラメータの値αとは異なる。ここでは、数７０に示す確率値αが計算され、そこから数７１に示す確率値（尤度）Ｐ（ｙ｜θt）が計算される。【０２０３】
【数７０】
2012008184.tif 000072
【０２０４】
【数７１】
2012008184.tif 000073
【０２０５】
次のステップＳ５は、レジスタ１から観測データを読み出して、ステップＳ３で設定された上式の確率量から、数７２に示す前向き確率を計算する部分である。なお、ここに示すβは確率値であり、前述したパラメータの値βとは異なる。ここでは、数７２に示す確率値βが計算される。【０２０６】
【数７２】
2012008184.tif 000074
【０２０７】
上記前向き確率と後向き確率は、計算回数を減らすために、既存のＨＭＭ推定アルゴリズムにも組み込まれていたものである。【０２０８】
続くステップＳ６は、ステップＳ４で計算された確率値αと、ステップＳ５で計算された確率値βを用い、レジスタ１から観測データを読み出して、数７３に示す状態遷移のカウント値と、数７４に示す出力のカウント値をそれぞれ計算する部分である。【０２０９】
【数７３】
2012008184.tif 000075
【０２１０】
【数７４】
2012008184.tif 000076
【０２１１】
次のステップＳ７は、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の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したことが、その効果として挙げられる。【０２１３】
本実施例では、配列が１本の観測データをレジスタ１に格納しているが、その場合のステップＳ７における確率用の更新式は、次式のようにあらわせる。【０２１４】
【数７５】
2012008184.tif 000077
【０２１５】
【数７６】
2012008184.tif 000078
【０２１６】
【数７７】
2012008184.tif 000079
【０２１７】
ステップＳ８は、計算の収束を判定する部分である。ここでは、直前のステップＳ７で計算された新たな確率量に基づく尤度Ｐを用い、その尤度Ｐと前記ステップＳ２における収束判定値との比較により、ステップＳ３〜ステップＳ９の計算が収束したか否かを判定する。近似的には、繰り返しの回数を指定して、指定回数に達したら、計算が収束したと判定してもよい。【０２１８】
ステップＳ８において計算が収束していないと判定された場合、ステップＳ９に移行して手順が繰り返される。ステップＳ９は、計算された確率量と事象のカウント値を、繰り返し計算のために更新する部分である。本実施例では、次式に示すように、２つの過去値をシフトさせる。【０２１９】
【数７８】
2012008184.tif 000080
【０２２０】
ここでは、繰り返しの回数が１つ増えることにより、t回目の確率量およびカウント値がt-1回目の確率量およびカウント値に更新され、t+1回目の確率量およびカウント値がt回目の確率量およびカウント値に更新される。なお、ステップＳ３とステップＳ９を一纏めにして、更新された確率量およびカウント値をそのまま次の繰り返しのために、推定装置１０のメモリに設定記憶させてもよい。【０２２１】
一方、ステップＳ８において計算が収束していると判定された場合、ステップＳ１０に移行して、次式に示す計算された三組の確率量を用いたＨＭＭを採用する。【０２２２】
【数７９】
2012008184.tif 000081
【０２２３】
上記数７９で計算した各値を用いて、推定装置１０は数６６に示すＨＭＭの確率構造を出力することができる。【０２２４】
図２は、図１に示すアルゴリズムを実現する推定装置１０の構成を示している。同図において、推定装置１０は前記レジスタ１を内蔵する入力手段１２と、レジスタ１に格納した観測データを適宜読み出し、上記ステップＳ１〜ステップＳ１０の手順を実行して数６６に示すＨＭＭの確率構造を推定する推定手段１４と、この推定手段１４で得たＨＭＭの確率構造を出力する出力手段１６とにより概ね構成される。本実施例の入力手段１２は、１本の配列を有する観測データをレジスタ１に格納できる構造となっているが、複数本の配列を有する観測データを取り扱う推定装置１００については、次の実施例で詳しく説明する。【０２２５】
推定装置１０のハードウェア構成は、例えば演算処理部としてＣＰＵを備えたコンピュータで実現することができる。その場合、図１に示す推定装置１０の処理手順を実行するプログラムが、メモリなどの記録媒体に記憶される。当該プログラムをどこに記憶するのかは限定せず、例えば通信手段を介してプログラムがコンピュータにダウンロードされる構成であってもよい。【０２２６】
推定手段１４は、前記ステップＳ１における高速化パラメータの値βや、ステップＳ２における初期の確率値および収束判定の条件を設定する初期設定手段２２と、ステップＳ３における確率量およびカウント値の設定や、ステップＳ９における確率量およびカウント値の更新を行なう更新設定手段２４と、ステップＳ４における前向き確率の計算や、ステップＳ５における後ろ向き確率の計算を行なうと共に、そこからステップＳ６における状態遷移のカウント値と出力のカウント値をそれぞれ計算し、さらにステップＳ７における新たな確率量の計算を行なう演算手段２６と、ステップＳ８における計算の収束判定を行ない、計算が収束していなければ、前記更新設定手段２４による確率量およびカウント値の更新を行なわせる一方で、計算が収束していれば、ステップＳ１０において、演算手段２６が直前に計算した新たな確率量を、ＨＭＭの最終的なパラメータ値として確定させる判定手段２８とを備えている。【０２２７】
また推定装置１０には、前記高速化パラメータの値や、初期の確率値および収束判定の条件や、更新設定される確率量およびカウント値の値を読み書き可能に格納するメモリ３０の他に、必要に応じて高速化パラメータの値や収束判定の条件を操作入力するキーボードやマウスなどの操作手段３２が、推定手段１４に接続して設けられる。これらの各装置構成により、上述したステップＳ１〜ステップＳ１０の手順が実行される。【０２２８】
次に、本実施例における試験用の観測データでの実験結果を説明する。【０２２９】
（Ａ）生成されたデータを用いた速度の評価
観測データの集合が強い凸性でなければ（そして強い凹面でもなければ）、全ての最適化アルゴリズムは局所的に最適となる。存在するデータの殆ど全ては、そのような性質を持っている。従来のあるいはlog-ＨＭＭアルゴリズムは例外的ではないし、alpha-ＨＭＭアルゴリズムも例外的ではない。したがって、そうしたＨＭＭアルゴリズムを適用する前に、発生のメカニズムが分かっている人工的なデータ集合を生成することが必要になってくる。【０２３０】
ここでは、数３３に示す数列をレジスタ１に格納する入力データとし、そのパラメータを次式のように指定する。【０２３１】
【数８０】
2012008184.tif 000082
【０２３２】
このマルコフ連鎖を使用し、[0, 1]上の一様な乱数によって、1000個のサンプルを有するデータ集合を生成した。完全な推定には、以下の困難があることを先に述べる。
（ａ）各データ集合の殆どが、局所的最適を与える対象である。
（ｂ）たとえ局所的最適が避けられたとしても、正確な推定には無限の長さの数列が必要である。
（ｃ）コンピュータが生成した乱数は単なる擬似的な乱数にすぎない。【０２３３】
したがって、以下が満たされれば満足するものとする。
（１）数３３に示すパラメータ集合の数値が、数８０からずれないこと
（２）alpha-ＨＭＭアルゴリズムの収束したパラメータ集合がlog-ＨＭＭのそれと近いこと
ステップＳ２における初期値は、次のように設定する。これにより、初期状態を除いて事前の情報は必要としないものとする。【０２３４】
【数８１】
2012008184.tif 000083
【０２３５】
以下の理由のため初期状態は固定した。
（ａ）分布された初期確率は、不安定な局所最適性を得やすい。これは、alpha-ＨＭＭアルゴリズムの各収束速度を比較するのに不要な曖昧さをもたらす。
（ｂ）ＨＭＭ推定の後の認識問題において、一つの初期状態はビタビ（Viterbi）アルゴリズムによって選択される。【０２３６】
図３は、対数尤度（log-likelihood）、すなわち繰り返し回数を考慮した尤度Ｐ（ｙ｜θt）に関する収束の傾向を示している。ここでは、数６０および数６１を満たす多様な値で実験した。【０２３７】
この図から以下のことがわかる。
（１）β=αcausalが増加するにつれて対数尤度の立ち上がりが早くなる。すなわち、収束するまでの繰り返し回数が少なくなる。
（２）alpha-ＥＭアルゴリズムの限界(α＝１)は、alpha-ＨＭＭアルゴリズムではβ=αcausal＝３に相当するが、この場合、発散してしまう。【０２３８】
さらに２つの視点で、alpha-ＨＭＭアルゴリズムの性能を比較する必要がある。一つは局所最適の比較であり、もう一つはＣＰＵ時間である。局所最適の比較では、次式に示すように、θに対する全パラメータと共に、それに付随した対数尤度の値を一覧にする必要がある。【０２３９】
【数８２】
2012008184.tif 000084
【０２４０】
上記の実験結果の一覧を見るとβを1〜2.75の範囲で変化させても対数尤度（ＬＬ）の値がほとんど変化せずに正確さを保っていることがわかる。この実験ではβ≧３に設定した実験は発散する結果になるために省略しているが、３に近い2.75に設定しても実験結果にそれほど影響なく高速化できることがわかる。ただし、このような実験結果の収束限界となるβの値は入力データによって若干異なる。【０２４１】
次は、収束するまでのＣＰＵ時間の比較である。この比較はソフトウェアの実行手段とハードウェアの手法に依存する。ここでは、alpha-ＨＭＭアルゴリズムを標準的なPOSIX環境におけるＣ‐コードとして実行する。これらのコードは標準的なシングルコアプロセッサで走らせている。前記ステップＳ８で行われる収束は、以下の判定基準で測定した。【０２４２】
【数８３】
2012008184.tif 000085
【０２４３】
上式において、Ｐnewは今回の計算で得た確率Ｐの値であり、Ｐoldは前回の計算で得た確率Ｐの値である。この収束判定基準は更新前後で確率の値がほとんど変化していなければ収束したと判断するものである。この収束判定式結果は次の表に示す通りであり、ここでは繰り返し回数とＣＰＵ時間による収束の比較を示している。【０２４４】
【表１】
2012008184.tif 000086
【０２４５】
局所最適性とＣＰＵ時間とを合せて比較すると、以下のことがわかる。
（ａ）収束する各対数尤度は非常に接近しているものの、何れも局所最大値が存在する。
（ｂ）対数のまたは従来のＨＭＭアルゴリズムは、常に最良の対数尤度をもたらす訳ではない。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な形態を使用したことにある。更新式の一例は、次式のように表せる。【０２５０】
【数８４】
2012008184.tif 000087
【０２５１】
ここで、θ＾ｔ＋１とθ＾ｔ（＾はθの上部に記される）は、より少ないサンプル（例えば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はθの添え字）を利用して、前記数６５の式で計算する構成となっている。【０２５５】
そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。【０２５６】
さらに、ここでの更新設定手段２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。【０２５７】
なお上述した本実施例の作用効果は、初期設定手段２２としての動作を実行する初期設定ステップと、更新設定手段２４としての動作を実行する更新設定ステップと、演算手段２６としての動作を実行する演算ステップと、判定手段２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。
【実施例２】
【０２５８】
図４は、本発明の第２実施例におけるプログラムの処理手順をあらわしたものである。【０２５９】
同図において、１０１は合計でＴ個の観測データを格納する記憶装置としてのレジスタで、各観測データのそれぞれは、時刻τが１から順に時系列に並んで格納される。ここでｙτは個々のデータ値を示し、ｙはＴ個のデータ値の集合を示し、ｎは何本目の配列かをあらわすインデックスを示しており、本実施例では複数の配列を有するＭ（＝２以上の整数）本の観測データをレジスタ１０１に格納している。推定装置１００は、レジスタ１０１に記憶されるＭ列のデータ値が、どのような確率構造（モデル）を有しているのかを、以下のステップＳ１１〜ステップＳ２０の各手順に従って推定解析するものである。【０２６０】
ステップＳ１１は、前述した高速化パラメータの値βを設定する部分である。これは上記第１実施例のステップＳ１に相当するもので、高速化パラメータの値βは、従来のβ＝１を特例として含む１≦β＜３の範囲に拡張される。なお、βが大きな値であるほど、推定装置１００としての処理の高速性は増すが、収束性を保持するにはβ＜３でなければならない。【０２６１】
ステップＳ１２は、初期確率と収束判定値を決める部分である。これは、前記ステップＳ２に相当するもので、後述するステップＳ１３〜ステップＳ１８の手順を繰り返す前に行なわれる。推定装置１００は、レジスタ１０１に格納された観測データを読み出して、その観測データが最も出現しやすくなる初期状態確率πと、状態遷移確率ａと、出力確率ｂとを推定するが、このステップＳ１２では、次式に示すように、それらの確率π，ａ，ｂの初期値を決定する。【０２６２】
【数８５】
2012008184.tif 000088
【０２６３】
上式において、θ０は数４３にあるように、０回目の繰り返しにおける初期の確率π，ａ，ｂの組合せを示し、推定装置１００は、その条件で状態ｉから状態ｊに移る状態遷移確率ａij｜θ０（０はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂjk｜θ0（０はθの添え字）の値と、最初に状態ｉとなる初期状態確率πｉ｜θ0（０はθの添え字）の値とをそれぞれ決定する。【０２６４】
またステップＳ１２では、対数尤度に基づく収束範囲を決めることで、ステップＳ１３〜ステップＳ１９の繰り返しを終了させるための収束判定値を決定する。この収束判定値はステップＳ１８で用いられ、具体的には前記数８３で示される。【０２６５】
これらの値をステップＳ１１，Ｓ１２で決定すると、ステップＳ１３の手順に移行して、確率量とカウント値の設定が行なわれる。ステップＳ１３は、最初にステップＳ１２で求めた初期確率を利用し、それ以降はステップＳ１９で更新された確率量とカウント値を利用して、数８６に示す各確率値と、数８７に示す各カウント値を、実際に推定装置１００のメモリ（図示せず）に設定する部分である。【０２６６】
【数８６】
2012008184.tif 000089
【０２６７】
【数８７】
2012008184.tif 000090
【０２６８】
次のステップＳ１４は、レジスタ１０１から観測データを読み出して、ステップＳ３で設定された上式の確率量から、数８８と数８９に示す前向き確率を計算する部分である。なお、ここに示すαは確率値であり、前述したパラメータの値αとは異なる。ここでは、数８８に示す確率値αが計算され、そこから数８９に示す確率値（尤度）Ｐ（ｙ(n)｜θt）が計算される。【０２６９】
【数８８】
2012008184.tif 000091
【０２７０】
【数８９】
2012008184.tif 000092
【０２７１】
次のステップＳ１５は、レジスタ１０１から観測データを読み出して、ステップＳ１３で設定された上式の確率量から、数９０に示す前向き確率を計算する部分である。なお、ここに示すβは確率値であり、前述したパラメータの値βとは異なる。ここでは、数９０に示す確率値βが計算される。【０２７２】
【数９０】
2012008184.tif 000093
【０２７３】
上記前向き確率と後向き確率は、計算回数を減らすために、既存のＨＭＭ推定アルゴリズムにも組み込まれていたものである。【０２７４】
続くステップＳ１６は、ステップＳ１４で計算された確率値αと、ステップＳ５で計算された確率値βを用い、レジスタ１０１から観測データを読み出して、数９１に示す状態遷移のカウント値と、数９２に示す出力のカウント値をそれぞれ計算する部分である。【０２７５】
【数９１】
2012008184.tif 000094
【０２７６】
【数９２】
2012008184.tif 000095
【０２７７】
次のステップＳ１７は、t+1回目に繰り返される確率π，ａ，ｂの組合せを条件として、状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値と、状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値とを、ステップＳ１４で求めたt回目の確率Ｐ（ｙ(n)｜θt）と、それよりも１回前の確率Ｐ（ｙ(n)｜θt-1）を利用し、且つステップＳ５で求めたt回目のカウント値Ｎ(n)ａjk｜θt，Ｎ(n)ｂjk｜θt（jkはａの添え字、tはθの添え字）と、それよりも１回前のカウント値Ｎ(n)ａjk｜θt-1，Ｎ(n)ｂjk｜θt-1（jkはａの添え字、t-1はθの添え字）を利用して計算する部分であり、これは前記数６３，数６４，６５に対応している。【０２７８】
このステップＳ１７では、前記ステップＳ１１で設定したパラメータの値βが用いられ、本アルゴリズムの主要な特徴部分となっている。特に、添え字ｔ+1の部分の確率を計算する上で、時間シフトした添え字tや添え字t-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したことが、その効果として挙げられる。【０２７９】
本実施例では、配列がＭ本の観測データをレジスタ１０１に格納しているが、その場合のステップＳ１７における確率用の更新式は、次式のようにあらわせる。【０２８０】
【数９３】
2012008184.tif 000096
【０２８１】
【数９４】
2012008184.tif 000097
【０２８２】
【数９５】
2012008184.tif 000098
【０２８３】
ステップＳ１８は、計算の収束を判定する部分である。ここでは、直前のステップＳ１７で計算された新たな確率量に基づく尤度Ｐを用い、その尤度Ｐと前記ステップＳ１２における収束判定値との比較により、ステップＳ１３〜ステップＳ１９の計算が収束したか否かを判定する。近似的には、繰り返しの回数を指定して、指定回数に達したら、計算が収束したと判定してもよい。【０２８４】
ステップＳ１８において計算が収束していないと判定された場合、ステップＳ１９に移行して手順が繰り返される。ステップＳ１９は、計算された確率量と事象のカウント値を、繰り返し計算のために更新する部分である。本実施例では、次式に示すように、２つの過去値をシフトさせる。【０２８５】
【数９６】
2012008184.tif 000099
【０２８６】
ここでは、繰り返しの回数が１つ増えることにより、t回目の確率量およびカウント値がt-1回目の確率量およびカウント値に更新され、t+1回目の確率量およびカウント値がt回目の確率量およびカウント値に更新される。【０２８７】
一方、ステップＳ１８において計算が収束していると判定された場合、ステップＳ２０に移行して、次式に示す計算された三組の確率量を用いたＨＭＭを採用する。【０２８８】
【数９７】
2012008184.tif 000100
【０２８９】
上記数９７で計算した各値を用いて、推定装置１０は数６６に示すＨＭＭの確率構造を出力することができる。【０２９０】
図５は、図４に示すアルゴリズムを実現する推定装置１００の構成を示している。同図において、推定装置１００は前記レジスタ１０１を内蔵する入力手段１１２と、レジスタ１０１に格納した観測データを適宜読み出し、上記ステップＳ１１〜ステップＳ２０の手順を実行して数６６に示すＨＭＭの確率構造を推定する推定手段１１４と、この推定手段１１４で得たＨＭＭの確率構造を出力する出力手段１１６とにより概ね構成される。本実施例の入力手段１１２は、Ｍ本の配列を有する観測データをレジスタ１０１に格納できる構造となっている。推定装置１００のハードウェア構成については、前記第１実施例の推定装置１０と同様であるため、ここでは説明を省略する。【０２９１】
推定手段１１４は、初期設定手段１２２と、更新設定手段１２４と、演算手段１２６と、判定手段１２８とを備えている。これらは前記第１実施例の初期設定手段２２，更新設定手段２４，演算手段２６および判定手段２８にそれぞれ対応するもので、取り扱う観測データがＭ本の配列になった以外は、第１実施例と同様に機能する。【０２９２】
また推定装置１１０には、前記高速化パラメータの値や、初期の確率値および収束判定の条件や、更新設定される確率量およびカウント値の値を読み書き可能に格納するメモリ１３０の他に、必要に応じて高速化パラメータの値や収束判定の条件を操作入力するキーボードやマウスなどの操作手段１３２が、推定手段１１４に接続して設けられる。これらの各装置構成により、上述したステップＳ１１〜ステップＳ２０の手順が実行される。【０２９３】
以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１０１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１１４とを備えたＨＭＭの推定装置１００において、推定手段１１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段１２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，出力確率ｂ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび出力の期待値Ｎｂを更新設定する更新設定手段１２４と、更新設定手段１２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数６８，数６９を参照）と共に、前記レジスタ１０１から読み出した観測データと、初期設定手段１２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数６３，数６４，数６５，数９３，数９４，数９５を参照）演算手段１２６と、演算手段１２６による計算の収束を判定し、計算が収束していなければ、演算手段１２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段１２６で計算した新たな各確率量を最終的な値として出力させる判定手段１２８とを備えている。【０２９４】
このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，出力確率ｂ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび出力の期待値Ｎｂを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。【０２９５】
また、ここでの演算手段１２６は、ｔ+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はデータの順序を表すインデックスである。【０３０５】
【数９８】
2012008184.tif 000101
【０３０６】
各ｙtは、スカラーあるいはベクトルである。上式は、単一配列の場合であるが、M本の複数配列で与えられるならば、ソースデータ列は次式の通りである。記号｛｝は集合を示している。【０３０７】
【数９９】
2012008184.tif 000102
【０３０８】
まず、単一配列では、数９８の観測データｙ（ｙは太字）が与えられる。この場合のHMMの課題は、次式の最尤推定法（Maximum Likelihood Estimation：MLE）の認識において、最良なモデルを見つけることである。【０３０９】
【数１００】
2012008184.tif 000103
【０３１０】
ここで、S（Sは太字）は、次式のような状態遷移系列を意味している。【０３１１】
【数１０１】
2012008184.tif 000104
【０３１２】
ランダムな変数とそれらの値のために、小文字のｙ（ｙは太字）およびｓ（ｓは太字）は同様に用いられる（通常、ランダム変数は大文字で示される）。【０３１３】
各確率は以下の数１０２〜数１０５の通りである。【０３１４】
・初期状態確率【０３１５】
【数１０２】
2012008184.tif 000105
【０３１６】
・状態遷移確率【０３１７】
【数１０３】
2012008184.tif 000106
【０３１８】
ここで、もし、状態ｓｉから状態ｓｊへのつながりが無ければ、状態ｉから状態ｊに移る状態遷移確率はａij＝０となる。これは、状態遷移のための前のトポロジーを反映している。【０３１９】
・出力確率【０３２０】
【数１０４】
2012008184.tif 000107
【０３２１】
ここで、上式の各要素は次式の通りである。【０３２２】
【数１０５】
2012008184.tif 000108
【０３２３】
上記数１０２〜数１０５に示した確率の集合をまとめてθとして次式に示す。【０３２４】
【数１０６】
2012008184.tif 000109
【０３２５】
そして、確率的なデータ構造について以下のように解釈している。不完全なデータはｙ（ｙは太字）であり、推定される消失データはｓ（ｓは太字）であり、完全なデータは、数１００の確率量と数１０６のパラメータを有しており、次式に示されるｘ（ｘは太字）として定義する。【０３２６】
【数１０７】
2012008184.tif 000110
【０３２７】
・Alpha-EMアルゴリズム
不完全なデータ、消失データ、完全なデータによるHMMの解釈は、EMアルゴリズムに匹敵する。本紙の目的は、新たなHMM推定アルゴリズムを見つけることであるため、alpha-EMアルゴリズムから始まる道筋を示している。【０３２８】
観測データｙ（ｙは太字）は推定すべきパラメータの全体を表すψによってパラメータ化されるため、Ｐy|ψ（ｙ｜ψ）（ｙは太字）を確率密度あるいは確率量とする。ｘ∈Ｘ（ｘとＸは太字）を、消失データを含む理想的な観測結果である完全なデータ、あるいは拡張されたデータとする。そして、不完全なデータ確率密度関数(pdf)、あるいは確率量関数(pmf)は次式で示される。【０３２９】
【数１０８】
2012008184.tif 000111
【０３３０】
ここで、積分する範囲を次式とする。【０３３１】
【数１０９】
2012008184.tif 000112
【０３３２】
pdfの積分は、pmfの総和となる。そして、条件付きのpdfあるいはpmfは次式である。【０３３３】
【数１１０】
2012008184.tif 000113
【０３３４】
alpha-EMアルゴリズムでは、次式のようなalpha対数が用いられる。【０３３５】
【数１１１】
2012008184.tif 000114
【０３３６】
ここで、α＝−１の場合は対数、すなわち、次式である。【０３３７】
【数１１２】
2012008184.tif 000115
【０３３８】
alpha-EMアルゴリズムは、alpha対数に関して、不完全なデータの尤度比を考慮する必要があるため、次式のようになる。【０３３９】
【数１１３】
2012008184.tif 000116
【０３４０】
ここで、φとψは、繰返し最大化ステップにおける数１０６に対して、古いモデルのパラメータと新しいモデルのパラメータを示している。それから、alpha-EMアルゴリズムの基本的な方程式が次式として得られる。【０３４１】
【数１１４】
2012008184.tif 000117
【０３４２】
ここで、D(α)は、常に負とならない２つの条件付きの確率Pｘ|ｙ,φ（x|y,φ）（ｘとｙは太字）とPｘ|ｙ,ψ（x|y,ψ）（ｘとｙは太字）の間のalphaダイバージェンスである。上式において、次式のＱ関数が重要である。Ｅは定量化（quantity）によって表される最大化のための目的関数を示している。【０３４３】
【数１１５】
2012008184.tif 000118
【０３４４】
数１１４のために、もし、このＱ関数が正ならば、a＜1の範囲で、数１１４の左辺である不完全なデータのalpha対数尤度比も正となる。したがって、alpha-EMアルゴリズムとその変化形であるalpha-GEMアルゴリズムについて以下に示す。【０３４５】
[Alpha-EMアルゴリズム]
初期化: 数１０６の初期値を選択し、φとして用いる。
E-Step: 数１１５の計算を実行する。
M-step: 更新パラメータを次式で算出する。【０３４６】
なお、"arg max"とは最大値を与える変数を意味する記号である。すなわち、次式は、Ｑ（α）（ｘ|ｙ、φ）（ψ｜φ）（ｘとｙは太字であり、（ｘ|ｙ、φ）はＱの下付き）を変数ψの関数と考えたとき、その最大値を与える変数の値をψ*とする、という意味である。【０３４７】
【数１１６】
2012008184.tif 000119
【０３４８】
U-step: φをψ*によって置き換え、収束することを確認する。収束しない場合は、E-stepに戻って繰り返しが繰返される。【０３４９】
[Alpha-GEMアルゴリズム]
これは、上記M-stepを次式のようにＱ関数を非負とするψ+の算出に置き換えたアルゴリズムである。【０３５０】
【数１１７】
2012008184.tif 000120
【０３５１】
まず、alpha-HMMアルゴリズムの近似バージョンがalpha-GEMアルゴリズムであることに着目すべきである。数１１４を有するalpha-EMアルゴリズムの長所によるalpha-HMM推定アルゴリズムの基本的な特性について説明する。【０３５２】
完全なデータをｘ＝（s，ｙ）（ｘ、ｓ、ｙは太字）とし、ｓ（ｓは太字）を消失データとし、ｙ（ｙは太字）を不完全なデータとし、数１１５は次式に相当するものとする。【０３５３】
【数１１８】
2012008184.tif 000121
【０３５４】
もし、上式の値が負にならなければ、以下の不等式が成り立つ。【０３５５】
【数１１９】
2012008184.tif 000122
【０３５６】
ここで、alpha-HMMの抽象バージョンに留意すべきである。【０３５７】
・単一配列のalpha-HMM 
Non-Causalな更新方程式
離散型のアルファベット系列ｙ（ｙは太字）の場合において、最大化されるQ関数は次式で示される。なお、ｌはアルゴリズムの繰返しの指標（インデックス）である。一方、上記文章中に出てくるｔはデータの順番に対応する指標（インデックス）である。したがって、ｌとｔは区別されるものである。【０３５８】
【数１２０】
2012008184.tif 000123
【０３５９】
まず、状態遷移確率aijの更新方程式について説明する。状態遷移確率aij は、その更新後に確率量とする必要があるので、ラグランジュの未定係数法（Lagrange multiplier）を使用する必要がある。したがって、次式のように最大値を算出するための微分を行なう。【０３６０】
【数１２１】
2012008184.tif 000124
【０３６１】
そうすると、上式から次式が得られる。【０３６２】
【数１２２】
2012008184.tif 000125
【０３６３】
ここで、Ｎij(ｓ)（ｓは太字）は、iからjへの状態遷移の数と位置である。次のステップでは、ソフトウェアで実行可能なアルゴリズムを目指すため、以下の問題を解決する必要がある。
(a)Non-causalityの存在: 右辺にθl+1(θは太字）が含まれている。このままでは、自己撞着性のため計算できない。
(b)右辺の算出には、状態数がＮでデータ数がＴのとき、Ｏ(ＮT)すなわちＮのＴ乗のオーダーの演算が必要である。【０３６４】
上記の問題は、次のセクションで解決されることとなるが、その前に、出力確率bjk|θl+1（l+1はθの添え字）と初期状態確率πi|θl+1（l+1はθの添え字）の２つの更なる更新式について次式に示す。なお、Ｎjk(ｓ)（ｓは太字）は、状態ｓt＝ｊでの出力yt＝kの発生回数である。【０３６５】
【数１２３】
2012008184.tif 000126
【０３６６】
【数１２４】
2012008184.tif 000127
【０３６７】
Causal近似と拡張系列: 離散出力の場合
更新方程式である数１２２〜数１２４のコアとなる部分は、繰返し指標（インデックス）のシフトによるCausal近似によって、次式のように変換することができる。【０３６８】
【数１２５】
2012008184.tif 000128
【０３６９】
したがって、P（ｙ|θl）＝P（ｙ|θl-1）＋o(1)(yは太字）の領域で次の等式を得る。【０３７０】
【数１２６】
2012008184.tif 000129
【０３７１】
ここで、αcausalという表記は、以後、βと表記する。【０３７２】
いま、数１２５が算出可能であるが、演算の複雑度のために、もう一つの近似が必要となる。このため、我々は次式の系列拡張を用いている。【０３７３】
【数１２７】
2012008184.tif 000130
【０３７４】
そして、causal近似式である数１２５の適用と数１２７の系列拡張は、以下に示す遷移確率の更新方程式を与える。なお、これらの確率の更新式である上記数１２２〜数１２４は、それぞれ上記実施例１の数７５〜７７に対応するものである。【０３７５】
単一配列の離散型alpha-HMMの遷移確率aij|θl+1(l+1は、θの添え字)は、次式である。【０３７６】
【数１２８】
2012008184.tif 000131
【０３７７】
ここで、以下の各特性について理解しておくことが重要である。
[特性１]
β＝αcausal＝１の場合は、従来のlog-HMMの方法に帰趨する。
[特性２]
数１２８の分子は、現在と過去の更新項の重み付けられた総和である。分母も同様である。
[特性３]
数１２８の第２行及び第３行は、これらの確率計算に伴う複雑度を抑えた従来の前方-後方アルゴリズムの方法に匹敵する。
[特性４]
alpha-HMMに唯一、追加する必要があるものは、θl-1での更新項を記憶することである。これは、CPU時間に直接現れているように繰返しを減らすことを意味する。実際の実験結果でも、この予想通りになっている。【０３７８】
次に、残りの２つの更新方程式を示す。【０３７９】
単一配列の離散型alpha-HMMの出力確率ｂjk|θl+1(l+1は、θの添え字)は、次式である。【０３８０】
【数１２９】
2012008184.tif 000132
【０３８１】
単一配列の離散型alpha-HMMの初期状態確率πi|θl+1(l+1は、θの添え字)は、次式である。【０３８２】
【数１３０】
2012008184.tif 000133
【０３８３】
（変形例）
以下、近似式の変形例について説明する。なお、以下の変形例は本発明の他の実施例でも同様に適用可能なものである。上述した隠れマルコフモデルの推定法では、有限個の過去情報を用いて各確率量をすべて次のような計算に基づいて計算している。【０３８４】
【数１３１】
2012008184.tif 000134
【０３８５】
この高速化パラメータβは、利用者が選ぶデザインパラメータであり、現在値による繰り返し値に対し、1回前の繰り返し値を重みづけする効果を有している。この過去の値に重みづけを行って加算することは、有限の過去にまでさかのぼって行うこともできる。すなわち、ετ≧０をデザインパラメータとして設定し、次式のように置換してもよい。【０３８６】
【数１３２】
2012008184.tif 000135
【０３８７】
【数１３３】
2012008184.tif 000136
【０３８８】
そうすると、数１３１は次式となる。【０３８９】
【数１３４】
2012008184.tif 000137
【０３９０】
次に展開による近似式を用いる場合について説明すると、上記数１３１や数１３４は、次式のように展開できる（例えばテーラー展開）。【０３９１】
【数１３５】
2012008184.tif 000138
【０３９２】
ここで、o(1)は高位の無限小を意味する記号であり、上式の右辺第一項は，従来の隠れマルコフモデル推定アルゴリズムに相当している。そして、モーメンタム項に相当する右辺第二項は、展開により近似される高速化項の意味をもつ。【０３９３】
なお、Ｔ＝１かつε１＝Ｐ（ｙ｜θl-1）／Ｐ（ｙ｜θl）（ｙはそれぞれ太字），すなわち重みを用いた場合は、重みを用いていない場合に比べて、収束速度はわずかに遅くなるものの、β＝３付近での収束性能の安定性が良くなる。換言すると、図３のような性能曲線において、計算が収束する立ち上がり部分では、高速化に伴い、波状の上下の変動が生じやすくなるが、その変動が少なくなる。【０３９４】
・複数配列のalpha-HMM
もし、HMMの予め設計されたトポロジーは、エルゴードなものであれば、単一の長いトレーニング系列ｙ（ｙは太字）は、十分である。もし、選択されたトポロジーが、吸収状態を有しているならば、複数のトレーニング系列を使うことが望まれる。このような複数系列に対するalpha-HMM推定の更新式は、上記の単一系列に対する方法を利用する形で得ることができる。【０３９５】
離散的なシンボルの場合において、Ｓ（Ｓは太字）を次式のようにＭ個の状態遷移系列の集合とする。【０３９６】
【数１３６】
2012008184.tif 000139
【０３９７】
すると、複数配列のＱ関数は次式である。【０３９８】
【数１３７】
2012008184.tif 000140
【０３９９】
ここで、Ｐは、次式に示されるMarkov過程の確率である。【０４００】
【数１３８】
2012008184.tif 000141
【０４０１】
上式において、初期状態確率πs0(n)(0はｓの下付きの添え字であり、（ｎ）は上付きの添え字である)、状態遷移確率ast-1(n)st(n) (t-1とtは、ｓの下付きの添え字であり、（ｎ）は上付きの添え字である)、及び出力確率ｂst(n)(yt(n)) (tは、ｓの下付きの添え字であり、（ｎ）は上付きの添え字である)の形態が、系列指標（インデックス）ｎから独立している点に着目すべきである。【０４０２】
初期状態確率πi|θl+1(l+1は、θの添え字)によって、初期状態の更新方程式の誘導は、Ｑ関数の数１３７の微分から始まる。上述の単一配列との違いは、上式からも明らかなように微分がｎ回現れることである。そして、次式のようにnon-causalな方程式が得られる。【０４０３】
【数１３９】
2012008184.tif 000142
【０４０４】
ここで、【０４０５】
【数１４０】
2012008184.tif 000143
【０４０６】
であり、上式の分母の関数ｆ(ｎ)は次式である。【０４０７】
【数１４１】
2012008184.tif 000144
【０４０８】
そして、causalの繰返し指標（インデックス）のシフトと系列の拡張は、次の更新方程式を与える。【０４０９】
【数１４２】
2012008184.tif 000145
【０４１０】
同様に、状態遷移と出力の更新方程式が得られる。【０４１１】
【数１４３】
2012008184.tif 000146
【０４１２】
ここで、【０４１３】
【数１４４】
2012008184.tif 000147
【０４１４】
及び【０４１５】
【数１４５】
2012008184.tif 000148
【０４１６】
また、次式の出力確率が得られる。【０４１７】
【数１４６】
2012008184.tif 000149
【０４１８】
ここで、上式で省略された各項はそれぞれ次式となる。【０４１９】
【数１４７】
2012008184.tif 000150
【０４２０】
【数１４８】
2012008184.tif 000151
【０４２１】
上記数１４７のΣの下にある記号「ｔ：ｙt＝ｋ」の意味は、ｙt＝ｋが成り立つｔのみを対象にする、という意味である。すなわち、上記数１４７の場合、ｙ＝｛ｙｔ｝T（ｔ＝１）（左辺のyは太字であり、ｔ＝１は｛ｙｔ｝の下付きの添え字）において、ｙt＝ｋが成り立つｔのみについて加算する、という意味である。【０４２２】
（Ｂ）連続系列の単一配列
ここから、連続的な出力系列の場合の本実施例の具体的な説明を行なう。【０４２３】
もし、出力系列ｙ＝{ｙt}Ｔ（ｔ＝１）（出力系列ｙは太字であり、ｔ＝１は｛ｙt｝の下付き）が、連続的な多変数の観測結果として現れれば、数１２８〜数１３０と似ているが、少しだけ異なる更新方程式が得られる。この場合、ｙtは太字の文字で示されていないが、適切な次元のユークリッド空間におけるベクトルである。【０４２４】
このような連続的なアルファベットの場合、最尤推定法（MLE）の問題は、以下の尤度を最大化することである。【０４２５】
【数１４９】
2012008184.tif 000152
【０４２６】
ここで、状態stで、出力確率ｂｓtｋt（ｙt）（ｓtｋtのtは、それぞれｓ，ｋの下付きの添え字）において、kt番目の枝へ遷移する確率を特定する分岐確率Ｃｓtｋt（ｓtｋtのtは、それぞれｓ，ｋの下付きの添え字）は、ｙtの確率密度関数である。我々は、これをガウス密度（Gaussian density）とみなしている。そして、状態ｊで状態ｋが出力される出力確率ｂｊｋ（ｙt）は次式となる。【０４２７】
【数１５０】
2012008184.tif 000153
【０４２８】
ここで、μｊｋは平均値ベクトルであり、Σｊｋは共分散行列であり、総和と混同してはいけない。平均値ベクトルとは、分岐確率の確率密度関数の平均値を示している（以後、同様）。そして、状態jでの出力確率密度関数(pdf)は、ｂｊ（ｙt）がpdfとなるように次式となる。【０４２９】
【数１５１】
2012008184.tif 000154
【０４３０】
ここで、上式の記号について説明する。【０４３１】
まず、離散系列の場合、ｋを出力記号の種類を表すインデックスとすると、出力ykのとる値は有限種類であるため、ｂj(yk)はｂｊｋと略記できる。しかし、連続系列の場合は、離散値ではなく、連続値であるため、そのような表記はできない。そのため、時間を表すインデックスをｔとして、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は、それぞれｓ，ｋの下付きの添え字）を除いた部分であるため、図６、図７の対応関係との矛盾は生じない。【０４３５】
混合確率の場合のため、消失データはｓ（ｓは太字）とｃ（ｃは太字）である。そのため、Ｑ関数は次式である。【０４３６】
【数１５２】
2012008184.tif 000155
【０４３７】
数１４９、数１５２から明らかなように、初期確率と状態遷移確率の更新方程式は、それぞれ数１２８、数１３０と同様である。分岐確率ｃｊｋの更新方程式は、状態遷移確率aijの場合と同様にラグランジュの未定係数法（Lagrange multiplier）によって得られる。【０４３８】
【数１５３】
2012008184.tif 000156
【０４３９】
そして、右辺は、状態遷移確率の場合と同様にcausalで計算可能となり、次式のように変形される。【０４４０】
【数１５４】
2012008184.tif 000157
【０４４１】
次式は、平均値ベクトルμｊｋの更新方程式である。μｊｋについて数１５２の直接微分として、次のnon-causalの方程式が得られる。【０４４２】
【数１５５】
2012008184.tif 000158
【０４４３】
そして、繰返し指標（インデックス）のシフト、系列の拡張、総和の変更により、次式が得られる。【０４４４】
【数１５６】
2012008184.tif 000159
【０４４５】
数１５６が、θl-1の過去情報が十分に利用できる点を示していることに着目すべきである。【０４４６】
共分散行列の更新には、行列微分が必要である。Ｑ関数の数１５２を共分散行列の逆行列Σ-1ｊｋ（ｊｋはΣの添え字）について微分することにより、次のnon-causalな方程式が得られる。【０４４７】
【数１５７】
2012008184.tif 000160
【０４４８】
そして、繰返し指標（インデックス）のシフト、系列の拡張、および総和の変更により、次の更新式が得られる。【０４４９】
【数１５８】
2012008184.tif 000161
【０４５０】
ここで、上式の各項は以下の通りである。【０４５１】
【数１５９】
2012008184.tif 000162
【０４５２】
【数１６０】
2012008184.tif 000163
【０４５３】
共分散行列の更新式である数１５８には、効果的に十分に利用できる過去情報の形態を有していることに着目すべきである。なお、記号Σ-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と同等である。【０４５６】
図８は、本発明の第３実施例におけるプログラムの処理手順をあらわしたものである。このフローチャート全体の流れは、前記第１実施例の場合とほぼ同様である。ただし、本実施例では、観測データｙが、連続的な多変数の観測結果として現れる単一配列の連続系列データである点で異なる。また、初期値として設定され、繰返し計算される未知パラメータには、初期状態確率、状態遷移確率、分岐確率、平均値ベクトル、共分散行列があるため、各ステップにおける計算式も異なる。【０４５７】
そのため、置き換わる数式の対応関係について説明する。まず、繰返し更新される未知パラメータの集合の組み合わせを次式とする。【０４５８】
【数１６１】
2012008184.tif 000164
【０４５９】
推定装置１０が最終的に算出しようとするHMMの確率構造は、数６６の代わりに数１４９となる。ステップＳ３２で決定する未知パラメータの初期値は、数６７の代わりに次式となる。【０４６０】
【数１６２】
2012008184.tif 000165
【０４６１】
ステップＳ３２における収束判定値の決定方法は、数８３と同じである。【０４６２】
ステップＳ３３では、最初にステップＳ３２で決定した未知パラメータの初期値を利用し、それ以降はステップＳ３９で更新された未知パラメータとカウント値を利用して、次の数１６３に示す各確率値と、数１６４に示す各カウント値を、実際に推定装置１０のメモリ（図示せず）に設定する。なお、前記実施例１、２と異なり、指標（インデックス）は、データの順番のｔではなく、アルゴリズムの繰り返しのｌを付して設定する。【０４６３】
【数１６３】
2012008184.tif 000166
【０４６４】
【数１６４】
2012008184.tif 000167
【０４６５】
その後のステップＳ３４〜Ｓ３６の演算動作は、前記実施例１の数７０〜数７４の一部の符号を置き換えて適用することによって行う。具体的には、出力確率ｂｊ．ｙτ+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の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したこと等は、前記実施例１と同様にその効果として挙げられる。【０４６８】
本実施例では、配列が１本の観測データをレジスタ１に格納しているが、その場合のステップＳ３７における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数１３０、状態遷移確率の更新式は数１２８、分岐確率の更新式は数１５４、平均値ベクトルの更新式は数１５６、共分散行列の更新式は数１５８〜数１６０である。【０４６９】
ステップＳ３８では、前記実施例１と同様に、前記ステップ３７で計算された新たな確率量に基づく尤度Ｐを用いて、その尤度ＰとステップＳ３２における収束判定値との比較により、ステップＳ３３〜Ｓ３９の計算が収束したか否かを判定する。収束してなければ、ステップＳ３９に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように２つの過去値をシフトさせる。【０４７０】
【数１６５】
2012008184.tif 000168
【０４７１】
ここでは、繰り返しの回数が１つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。なお、ステップＳ３３とステップＳ３９を一纏めにして、更新された確率量およびカウント値をそのまま次の繰り返しのために、推定装置１０のメモリに設定記憶させてもよい。【０４７２】
一方、ステップＳ３８において計算が収束していると判定された場合、ステップＳ４０に移行して、次式に示す計算された５組の確率量を用いたＨＭＭを採用する。【０４７３】
【数１６６】
2012008184.tif 000169
【０４７４】
上記数１６６で計算した各値を用いて、推定装置１０は数１４９に示すＨＭＭの確率構造を出力することができる。【０４７５】
図８に示すアルゴリズムを実現する推定装置の構成は、図２に示した前記第１実施例の推定装置１０と同様であるため、ここでは説明を省略する。【０４７６】
以上のように本実施例では、入力される観測データを時系列に格納する記憶手段としてのレジスタ１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１４とを備えたＨＭＭの推定装置１０において、推定手段１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを更新設定する更新設定手段２４と、更新設定手段２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数１６３，数１６４を参照）と共に、前記レジスタ１から読み出した観測データと、初期設定手段２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１２８，数１３０，数１５４，数１５６，数１５８〜数１６０を参照）演算手段２６と、演算手段２６による計算の収束を判定し、計算が収束していなければ、演算手段２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段２６で計算した新たな各確率量を最終的な値として出力させる判定手段２８とを備えている。【０４７７】
このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。【０４７８】
また、ここでの演算手段２６は、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の更新方程式が得られる。初期状態確率と状態遷移確率の更新方程式は、それぞれ数１４２〜数１４５と同様である。しかし、出力の更新方程式の集合は、離散的なシンボルの場合と異なる。我々は、分岐確率、平均値ベクトル、及び共分散行列の更新方程式を必要としている。【０４８５】
分岐確率については、初期状態確率と状態遷移確率の場合と同様にラグランジュの未定係数法を用いて算出可能である。そして、次式の更新方程式が得られる。【０４８６】
【数１６７】
2012008184.tif 000170
【０４８７】
ここで、各要素は次式である。【０４８８】
【数１６８】
2012008184.tif 000171
【０４８９】
【数１６９】
2012008184.tif 000172
【０４９０】
Ｑ関数の数１３７での直接ベクトル微分であるμｊｋ｜θl+1（l+1はθの下付きの添え字）が適用される。その更新方程式は次式である。【０４９１】
【数１７０】
2012008184.tif 000173
【０４９２】
ここで、【０４９３】
【数１７１】
2012008184.tif 000174
【０４９４】
数１７０、数１７１と同様に、共分散行列の更新方程式は、Σ-1ｊｋ｜θl（ｊｋ｜θlはΣの下付き添え字であり、さらにlはθの下付き添え字）について行列微分を用いることによって得られる。【０４９５】
【数１７２】
2012008184.tif 000175
【０４９６】
ここで、Ｇ(n)θl（θlはＧの下付き添え字）は次式である。【０４９７】
【数１７３】
2012008184.tif 000176
【０４９８】
図９は、本発明の第４実施例におけるプログラムの処理手順をあらわしたものである。このフローチャート全体の流れは、前記第３実施例の場合とほぼ同様である。ただし、本実施例では、観測データｙが、連続的な多変数の観測結果として現れる複数配列の連続系列データである点で異なる。【０４９９】
置き換わる数式の対応関係について説明する。まず、推定装置１００が最終的に算出しようとするHMMの確率構造は、数６６の代わりに数１４９となる。ステップＳ４２で決定する未知パラメータの初期値は、数１６２と同じである。【０５００】
ステップＳ４２における収束判定値の決定方法は、数８３と同じである。【０５０１】
ステップＳ４３では、最初にステップＳ４２で決定した未知パラメータの初期値を利用し、それ以降はステップＳ４９で更新された未知パラメータとカウント値を利用して、数１７４に示す各確率値と、数１７５に示す各カウント値を、実際に推定装置１００のメモリ（図示せず）に設定する。なお、指標（インデックス）は、データの順番のｔではなく、アルゴリズムの繰り返しのｌを付して設定する。【０５０２】
【数１７４】
2012008184.tif 000177
【０５０３】
【数１７５】
2012008184.tif 000178
【０５０４】
その後のステップＳ４４〜Ｓ４６の演算動作は、前記実施例２の数８８〜数９２の一部の符号を置き換えて適用することによって行う。具体的には、出力確率ｂｊ．ｙτ+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の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したこと等は、前記実施例２と同様にその効果として挙げられる。【０５０７】
本実施例では、配列がＭ本の観測データをレジスタ１０１に格納しているが、その場合のステップＳ４７における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数１４２、状態遷移確率の更新式は数１４３〜数１４５、分岐確率の更新式は数１６７〜数１６９、平均値ベクトルの更新式は数１７０、数１７１、共分散行列の更新式は数１７２、数１７３である。【０５０８】
ステップＳ４８では、前記実施例２と同様に、前記ステップＳ４７で計算された新たな確率量に基づく尤度Ｐを用いて、その尤度ＰとステップＳ４２における収束判定値との比較により、ステップＳ４３〜Ｓ４９の計算が収束したか否かを判定する。収束してなければ、ステップＳ４９に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように２つの過去値をシフトさせる。【０５０９】
【数１７６】
2012008184.tif 000179
【０５１０】
ここでは、繰り返しの回数が１つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。【０５１１】
一方、ステップＳ４８において計算が収束していると判定された場合、ステップＳ５０に移行して、次式に示す計算された５組の確率量を用いたＨＭＭを採用する。【０５１２】
【数１７７】
2012008184.tif 000180
【０５１３】
上記数１７７で計算した各値を用いて、推定装置１０１は数１４９に示すＨＭＭの確率構造を出力することができる。【０５１４】
図９に示すアルゴリズムを実現する推定装置の構成は、図５に示した前記第２実施例の推定装置１０１と同様であるため、ここでは説明を省略する。【０５１５】
以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１０１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１１４とを備えたＨＭＭの推定装置１００において、推定手段１１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段１２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを更新設定する更新設定手段１２４と、更新設定手段１２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数１７４，数１７５を参照）と共に、前記レジスタ１０１から読み出した観測データと、初期設定手段１２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１４２〜数１４５，数１６７〜数１７３を参照）演算手段１２６と、演算手段１２６による計算の収束を判定し、計算が収束していなければ、演算手段１２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段１２６で計算した新たな各確率量を最終的な値として出力させる判定手段１２８とを備えている。【０５１６】
このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。【０５１７】
また、ここでの演算手段１２６は、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)離散型の場合のbjkの役割を連続モードのcjkに割り当てる。また、平均値ベクトルと共分散行列は遷移状態jに依存しないという場合を考慮する。そして、図６の離散型の場合をまっすぐ延長した構造となる。これを半連続HMMと呼んでいる。この構造のもう一つの解釈は、MLE-VQ HMM (Maximum Likelihood Vector Quantization HMM)である。【０５２３】
半連続のalpha-HMMモデルは、数１４９を変形し、次式となる。【０５２４】
【数１７８】
2012008184.tif 000181
【０５２５】
したがって、半連続のalpha-HMMの更新方程式は以下のようになる。
[単一配列の半連続のalpha-HMMの初期状態確率]
更新式は、数１３０である。
[単一配列の半連続のalpha-HMMの状態遷移確率]
更新式は、数１２８である。
[単一配列の半連続のalpha-HMMの分岐確率]
更新式は、数１５４である。
[単一配列の半連続のalpha-HMMの平均値ベクトル]
次式のように、数１５６において、μｊｋ｜θl+1（l+1はθの下付きの添え字）をμｊ｜θl+1（l+1はθの下付きの添え字）とし、右辺のkt=kの項を除去したものである。これは、状態ｊに依存しないためである（以後、同様である）。【０５２６】
【数１７９】
2012008184.tif 000182
【０５２７】
[単一配列の半連続のalpha-HMMの共分散行列]
次式のように、数１５８〜数１６０において、Σｊｋ｜θl+1（l+1はθの添え字），μｊｋ｜θl（lはθの添え字），μｊｋ｜θl-1（l-1はθの添え字）を、それぞれΣｊ｜θl+1（l+1はθの添え字），μｊ｜θl（lはθの添え字），μｊ｜θl-1（l-1はθの添え字）とし、右辺のkt=kの項を除去したものである。【０５２８】
【数１８０】
2012008184.tif 000183
【０５２９】
ここで、各要素は次式である。【０５３０】
【数１８１】
2012008184.tif 000184
【０５３１】
【数１８２】
2012008184.tif 000185
【０５３２】
本発明の第５実施例におけるプログラムの処理手順は、前記第３実施例の場合とほぼ同様であるため、図８を用いて説明する。ただし、本実施例では、観測データｙが、連続的な多変数の観測結果として現れる単一配列の半連続系列データである点で異なる。【０５３３】
本実施例は、前記第３実施例の場合の特例と考えることもできるので、置き換わる数式の対応関係に絞って説明する。まず、繰返し更新される未知パラメータの集合の組み合わせを数１６１とし、推定装置１０が最終的に算出しようとするHMMの確率構造は、数６６の代わりに数１７８とする。ステップＳ３２で決定する未知パラメータの初期値は、数６７の代わりに次式となる。【０５３４】
【数１８３】
2012008184.tif 000186
【０５３５】
ステップＳ３２における収束判定値の決定方法は、数８３と同じである。【０５３６】
ステップＳ３３では、最初にステップＳ３２で決定した未知パラメータの初期値を利用し、それ以降はステップＳ３９で更新された未知パラメータとカウント値を利用して、次の数１８４に示す各確率値と、数１８５に示す各カウント値を、実際に推定装置１０のメモリ（図示せず）に設定する。【０５３７】
【数１８４】
2012008184.tif 000187
【０５３８】
【数１８５】
2012008184.tif 000188
【０５３９】
その後のステップＳ３４〜Ｓ３６の演算動作は、前記実施例３と同様に前記実施例１の数７０〜数７４の一部の符号を置き換えて適用することによって行なう。【０５４０】
ステップＳ３７では、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の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したこと等は、前記実施例３と同様にその効果として挙げられる。【０５４２】
本実施例では、配列が１本の観測データをレジスタ１に格納しているが、その場合のステップＳ３７における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数１３０、状態遷移確率の更新式は数１２８、分岐確率の更新式は数１５４、平均値ベクトルの更新式は数１７９、共分散行列の更新式は数１８０〜数１８２である。【０５４３】
ステップＳ３８では、前記実施例３と同様に、前記ステップ３７で計算された新たな確率量に基づく尤度Ｐを用いて、その尤度ＰとステップＳ３２における収束判定値との比較により、ステップＳ３３〜Ｓ３９の計算が収束したか否かを判定する。収束してなければ、ステップＳ３９に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように２つの過去値をシフトさせる。【０５４４】
【数１８６】
2012008184.tif 000189
【０５４５】
ここでは、繰り返しの回数が１つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。なお、ステップＳ３３とステップＳ３９を一纏めにして、更新された確率量およびカウント値をそのまま次の繰り返しのために、推定装置１０のメモリに設定記憶させてもよい。【０５４６】
一方、ステップＳ３８において計算が収束していると判定された場合、ステップＳ４０に移行して、次式に示す計算された５組の確率量を用いたＨＭＭを採用する。【０５４７】
【数１８７】
2012008184.tif 000190
【０５４８】
上記数１８７で計算した各値を用いて、推定装置１０は数１７８に示すＨＭＭの確率構造を出力することができる。【０５４９】
図８に示すアルゴリズムを実現する推定装置の構成は、図２に示した前記第３実施例の推定装置１０と同様であるため、ここでは説明を省略する。【０５５０】
以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１４とを備えたＨＭＭの推定装置１０において、推定手段１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを更新設定する更新設定手段２４と、更新設定手段２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数１８４，数１８５を参照）と共に、前記レジスタ１から読み出した観測データと、初期設定手段２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１３０，数１２８，数１５４，数１７９〜数１８２を参照）演算手段２６と、演算手段２６による計算の収束を判定し、計算が収束していなければ、演算手段２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段２６で計算した新たな各確率量を最終的な値として出力させる判定手段２８とを備えている。【０５５１】
このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。【０５５２】
また、ここでの演算手段２６は、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はθの添え字）を利用して、前記数１５４、数１７９〜数１８２の式で計算する構成となっている。【０５５３】
そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。【０５５４】
さらに、ここでの更新設定手段２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。【０５５５】
なお上述した本実施例の作用効果は、初期設定手段２２としての動作を実行する初期設定ステップと、更新設定手段２４としての動作を実行する更新設定ステップと、演算手段２６としての動作を実行する演算ステップと、判定手段２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。
【実施例６】
【０５５６】
本実施例では、半連続的なシンボルで複数配列のデータ系列の場合の実施形態について説明する。【０５５７】
この場合の複数配列の更新方程式は、平均ベクトルと共分散の状態依存性を制限することによって得られる。初期状態確率と状態遷移確率と分岐確率の更新方程式は、数１４２〜数１４５、および数１６７〜数１６９と同様である。したがって、次式のようになる。ここで、各記号の上にバー(-)を付し、バー(-）のない場合と記号を区別している。これは、例えば、数１７１と以下の数１８９では数式が異なっており、同じ記号が使えないためである。すなわち、バー(-）それ自体に特別な意味はない。【０５５８】
また、本実施例では、前記第５実施例において添え字がｊだった記号は、添え字がｋになっている。これは、kt=kの項の代わりにｓｔ＝ｊの項が除去されたためである。（以後、同様である）。【０５５９】
【数１８８】
2012008184.tif 000191
【０５６０】
【数１８９】
2012008184.tif 000192
【０５６１】
上記数１８８と同様に、共分散行列の更新方程式は、状態依存性を除去することによって次式のように得られる。【０５６２】
【数１９０】
2012008184.tif 000193
【０５６３】
【数１９１】
2012008184.tif 000194
【０５６４】
本発明の第６実施例におけるプログラムの処理手順は、前記第４実施例の場合とほぼ同様であるため、図９を用いて説明する。ただし、本実施例では、観測データｙが、連続的な多変数の観測結果として現れる複数配列の半連続系列データである点で異なる。【０５６５】
本実施例は、前記第４実施例の場合の特例と考えることもできるので、置き換わる数式の対応関係に絞って説明する。まず、繰返し更新される未知パラメータの集合の組み合わせを数１６１とし、推定装置１００が最終的に算出しようとするHMMの確率構造は、数６６の代わりに数１７８とする。ステップＳ４２で決定する未知パラメータの初期値は、数１８２と同じである。【０５６６】
ステップＳ４２における収束判定値の決定方法は、数８３と同じである。【０５６７】
ステップＳ４３では、最初にステップＳ４２で決定した未知パラメータの初期値を利用し、それ以降はステップＳ４９で更新された未知パラメータとカウント値を利用して、次の数１９２に示す各確率値と、数１９３に示す各カウント値を、実際に推定装置１００のメモリ（図示せず）に設定する。【０５６８】
【数１９２】
2012008184.tif 000195
【０５６９】
【数１９３】
2012008184.tif 000196
【０５７０】
その後のステップＳ４４〜Ｓ４６の演算動作は、前記実施例４と同様に前記実施例２の数８８〜数９２の一部の符号を置き換えて適用することによって行なう。【０５７１】
ステップＳ４７では、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の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したこと等は、前記実施例４と同様にその効果として挙げられる。【０５７３】
本実施例では、配列がＭ本の観測データをレジスタ１０１に格納しているが、その場合のステップＳ４７における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数１４２、状態遷移確率の更新式は数１４３〜数１４５、分岐確率の更新式は数１６７〜数１６９、平均値ベクトルの更新式は数１８８、数１８９、共分散行列の更新式は数１９０、数１９１である。【０５７４】
ステップＳ４８では、前記実施例４と同様に、前記ステップＳ４７で計算された新たな確率量に基づく尤度Ｐを用いて、その尤度ＰとステップＳ４２における収束判定値との比較により、ステップＳ４３〜Ｓ４９の計算が収束したか否かを判定する。収束してなければ、ステップＳ４９に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように２つの過去値をシフトさせる。【０５７５】
【数１９４】
2012008184.tif 000197
【０５７６】
ここでは、繰り返しの回数が１つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。【０５７７】
一方、ステップＳ４８において計算が収束していると判定された場合、ステップＳ５０に移行して、次式に示す計算された５組の確率量を用いたＨＭＭを採用する。【０５７８】
【数１９５】
2012008184.tif 000198
【０５７９】
上記数１９５で計算した各値を用いて、推定装置１０１は数１７８に示すＨＭＭの確率構造を出力することができる。【０５８０】
図９に示すアルゴリズムを実現する推定装置の構成は、図５に示した前記第２実施例の推定装置１０１と同様であるため、ここでは説明を省略する。【０５８１】
以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１０１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１１４とを備えたＨＭＭの推定装置１００において、推定手段１１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段１２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ-(-はμの上に付く），共分散行列Σ-(-はΣの上に付く），初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを更新設定する更新設定手段１２４と、更新設定手段１２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数１９２，数１９３を参照）と共に、前記レジスタ１０１から読み出した観測データと、初期設定手段１２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１４２〜数１４５，数１６７〜数１６９，数１８８〜数１９１を参照）演算手段１２６と、演算手段１２６による計算の収束を判定し、計算が収束していなければ、演算手段１２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段１２６で計算した新たな各確率量を最終的な値として出力させる判定手段１２８とを備えている。【０５８２】
このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ-(-はμの上に付く），共分散行列Σ-(-はμの上に付く），初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
前記数１６７〜数１７３の式で計算する構成となっている。【０５８３】
また、ここでの演算手段１２６は、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の６つのタイプ、すなわち、{離散的,連続的,半連続的}×{単一配列,複数配列}について述べた。本実施例は、区切りのある連続アルファベットの場合である。これは、離散的な文字と連続的な文字が混在している場合の解釈である。この離散と連続が混在している更新方程式は、数１５４、数１６９にある総和を２段の総和に変形し、ゼロ出力の許容によって得られる。例えば、数１５４の分母の総和は、次式のような２段の総和に分割される。【０５８８】
【数１９６】
2012008184.tif 000199
【０５８９】
ここで、Dgは、部分が重複しない集合である。Ｇは、その部分の濃度である。もし、Dgが連続的なアルファベットの特定のサブクラスに相当するならば、このことは、そのアルファベットが離散的なシンボルを伴うものと見なすことができる。換言すると、データが連続している各区間において、それぞれ総和を取り、さらに当該総和どうしについて総和を取っている。この離散的で連続的なアルファベットの場合にも、上記実施例のように単一配列の場合と複数配列の場合がある。したがって、この申請は、全部で８つのタイプのalpha-HMMを有している。【０５９０】
なお、上式の（・）という記号は、Σの中を省略し、総和の範囲を示すΣだけについて考察するために用いたものである。【０５９１】
この場合には、離散値は連続値を出す集団をグループ化した際のラベルと考えればよい。すなわち、データ系列は次式のように示される。【０５９２】
【数１９７】
2012008184.tif 000200
【０５９３】
ここで、ラベルは、次式として考える。【０５９４】
【数１９８】
2012008184.tif 000201
【０５９５】
上式のＧは、Ｇ≦Ｋ（分岐枝の数）である。そして、新たにｙt-（-はｙtの上に付く）をｙtと書き直せば、ΣTt=1（・）（TはΣの直上に付き、ｔ＝１は直下に付く）がグループ化に対応する部分を含んでいる場合には，これをΣt∈Dg（・）（t∈DgがΣの直下に付く）に変更するだけでよい．ただし、Dg＝｛t|ct＝ｇ｝である。【０５９６】
以下、各パラメータの更新式について述べる。
・離散連続混在系列の単数配列の場合
[単一配列の離散連続混在系列のalpha-HMMの初期状態確率]
初期状態確率は数１３０と同じである。
[単一配列の離散連続混在系列のalpha-HMMの状態遷移確率]
状態遷移確率は数１２８と同じである。
[単一配列の離散連続混在系列のalpha-HMMの分岐確率]
分岐確率はグループ化情報を併せもち、数１５４に代わって次式となる。以後、これを単にグループ化確率とよぶ。【０５９７】
【数１９９】
2012008184.tif 000202
【０５９８】
[単一配列の離散連続混在系列のalpha-HMMの平均値ベクトル]
各グループに対する平均値ベクトルは数１５６に代わって次式となる。【０５９９】
【数２００】
2012008184.tif 000203
【０６００】
[単一配列の離散連続混在系列のalpha-HMMの共分散行列]
各グループに対する共分散行列は、数１５８〜数１６０に代わって、次式となる。【０６０１】
【数２０１】
2012008184.tif 000204
【０６０２】
各要素は次式である。【０６０３】
【数２０２】
2012008184.tif 000205
【０６０４】
【数２０３】
2012008184.tif 000206
【０６０５】
・離散連続混在系列の複数配列の場合
[複数配列の離散連続混在系列のalpha-HMMの初期状態確率]
初期状態確率は数１４２と同じである。
[複数配列の離散連続混在系列のalpha-HMMの状態遷移確率]
状態遷移確率は数１４３〜数１４５と同じである。
[複数配列の離散連続混在系列のalpha-HMMの分岐確率]
分岐確率すなわちグループ化確率は、数１６７〜数１６９の変形として次式のようになる。【０６０６】
【数２０４】
2012008184.tif 000207
【０６０７】
各要素は次式である。【０６０８】
【数２０５】
2012008184.tif 000208
【０６０９】
【数２０６】
2012008184.tif 000209
【０６１０】
[複数配列の離散連続混在系列のalpha-HMMの平均値ベクトル]
平均値ベクトルは、数１７０、数１７１の変形として次式となる。【０６１１】
【数２０７】
2012008184.tif 000210
【０６１２】
各要素は次式である。【０６１３】
【数２０８】
2012008184.tif 000211
【０６１４】
[複数配列の離散連続混在系列のalpha-HMMの共分散行列]
共分散行列は、数１７２、数１７３の変形として次式となる。【０６１５】
【数２０９】
2012008184.tif 000212
【０６１６】
各要素は次式である。【０６１７】
【数２１０】
2012008184.tif 000213
【０６１８】
本実施例も上記第３〜６実施例と同様に、符号や未知パラメータの更新式が異なるだけであり、それらの数式を、単一配列の場合には、図２に示すハードウェア構成と図８に示すフローチャートに適用し、複数配列の場合には、図５に示すハードウェア構成と図９に示すフローチャートに適用するものであるため、ここでは、詳細な説明は省略する。なお、繰返し更新される未知パラメータの集合の組み合わせは数１６１とし、各更新式は、上述したように単一配列の場合は、数１３０、数１２８、数１９９〜２０３であり、複数配列の場合は、数１４２〜数１４５、数２０４〜数２１０である。【０６１９】
なお、本発明は上記実施例に限定されるものではなく、本発明の趣旨を逸脱しない範囲で変更可能である。例えば、上記各実施例の近似計算において、第３実施例で示した変形例のように遡る過去値を増やしたり、高次の近似にしたりすることも可能である。【０６２０】
以上のように本実施例（単一配列の場合もほとんど同様であるため、ここでは省略する）においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１０１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１１４とを備えたＨＭＭの推定装置１００において、推定手段１１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段１２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを更新設定する更新設定手段１２４と、更新設定手段１２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記レジスタ１０１から読み出した観測データと、初期設定手段１２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１４２〜数１４５，数２０４〜数２１０を参照）演算手段１２６と、演算手段１２６による計算の収束を判定し、計算が収束していなければ、演算手段１２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段１２６で計算した新たな各確率量を最終的な値として出力させる判定手段１２８とを備えている。【０６２１】
このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。【０６２２】
また、ここでの演算手段１２６は、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-ＨＭＭアルゴリズムに取り込んで癌になる確率を計算する等）、或いはロボットの動作認識などにも適用することができる。
【図１】
2012008184.tif 000354
【図２】
2012008184.tif 000355
【図３】
2012008184.tif 000356
【図４】
2012008184.tif 000357
【図５】
2012008184.tif 000358
【図６】
2012008184.tif 000359
【図７】
2012008184.tif 000360
【図８】
2012008184.tif 000361
【図９】
2012008184.tif 000362
【国際調査報告】
2012008184.tif 000363
2012008184.tif 000364
2012008184.tif 000365
2012008184.tif 000366