(57)【要約】  隠れ状態数および観測確率の種類と共にモデルの候補数が指数的に増加しても高速にモデル選択を実現できる隠れ変数モデル推定装置を提供する。変分確率計算部７１は、周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値を最大化することによって変分確率を計算する。モデル推定部７２は、各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定する。収束判定部７３は、変分確率計算部７１が変分確率を計算する際に用いた基準値が収束したか否かを判定する。
2013179579.tif 000002
【特許請求の範囲】
【請求項１】
周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値を最大化することによって変分確率を計算する変分確率計算部と、
各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定するモデル推定部と、
変分確率計算部が変分確率を計算する際に用いた基準値が収束したか否かを判定する収束判定部とを備える
ことを特徴とする隠れ変数モデル推定装置。
【請求項２】
変分確率計算部が変分確率を計算し、モデル推定部が最適な隠れ変数モデルを推定し、収束判定部が、基準値が収束したか否かを判定するループ処理を繰り返し、
基準値が収束したときに、当該基準値と、その１つ前のループ処理における基準値とを比較し、大きい方の基準値に対応する隠れ変数モデルを最適な隠れ変数モデルとして選択する最適モデル選択部を備える
請求項１に記載の隠れ変数モデル推定装置。
【請求項３】
変分確率計算部が変分確率を計算し、モデル推定部が最適な隠れ変数モデルを推定し、収束判定部が、基準値が収束したか否かを判定するループ処理を繰り返し、
変分確率計算部の計算結果に応じて、所定の条件を満たす隠れ状態を除去する隠れ状態除去部を備える
請求項１に記載の隠れ変数モデル推定装置。
【請求項４】
モデル推定部は、隠れ変数モデルとして隠れマルコフモデルを推定する請求項１から請求項３のうちのいずれか１項に記載の隠れ変数モデル推定装置。
【請求項５】
周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値を最大化することによって変分確率を計算し、
各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定し、
変分確率を計算する際に用いた基準値が収束したか否かを判定する
ことを特徴とする隠れ変数モデル推定方法。
【請求項６】
変分確率を計算し、最適な隠れ変数モデルを推定し、基準値が収束したか否かを判定するループ処理を繰り返し、
基準値が収束したときに、当該基準値と、その１つ前のループ処理における基準値とを比較し、大きい方の基準値に対応する隠れ変数モデルを最適な隠れ変数モデルとして選択する
請求項５に記載の隠れ変数モデル推定方法。
【請求項７】
変分確率を計算し、最適な隠れ変数モデルを推定し、基準値が収束したか否かを判定するループ処理を繰り返し、
変分確率の計算結果に応じて、所定の条件を満たす隠れ状態を除去する
請求項５に記載の隠れ変数モデル推定方法。
【請求項８】
隠れ変数モデルとして隠れマルコフモデルを推定する請求項５から請求項７のうちのいずれか１項に記載の隠れ変数モデル推定方法。
【請求項９】
コンピュータに、
周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値を最大化することによって変分確率を計算する変分確率計算処理、
各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定するモデル推定処理、および、
変分確率計算処理で変分確率を計算する際に用いた基準値が収束したか否かを判定する収束判定処理
を実行させるための隠れ変数モデル推定プログラム。
【請求項１０】
コンピュータに、
変分確率計算処理、モデル推定処理、収束判定処理のループ処理を繰り返し実行させ、
基準値が収束したときに、当該基準値と、その１つ前のループ処理における基準値とを比較し、大きい方の基準値に対応する隠れ変数モデルを最適な隠れ変数モデルとして選択する最適モデル選択処理
を実行させる請求項９に記載の隠れ変数モデル推定プログラム。
【請求項１１】
コンピュータに、
変分確率計算処理、モデル推定処理、収束判定処理のループ処理を繰り返し実行させ、
変分確率計算処理の計算結果に応じて、所定の条件を満たす隠れ状態を除去する隠れ状態除去処理
を実行させる請求項９に記載の隠れ変数モデル推定プログラム。
【請求項１２】
コンピュータに、
モデル推定処理で、隠れ変数モデルとして隠れマルコフモデルを推定させる
請求項９から請求項１１のうちのいずれか１項に記載の隠れ変数モデル推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、系列依存性を持つ多変量データの隠れ変数モデル推定装置、隠れ変数モデル推定方法、および隠れ変数モデル推定プログラムに関し、特に、モデル事後確率を近似し、その下界を最大化する事によって系列依存性を持つ多変量データの隠れ変数モデルを推定する隠れ変数モデル推定装置、隠れ変数モデル推定方法、および隠れ変数モデル推定プログラムに関する。
【背景技術】
【０００２】
系列依存性を持つデータが種々存在する。系列依存性を持つデータの例として、例えば、時間依存性を有するデータや、文字系列に依存する文章、塩基系列に依存する遺伝子データ等が挙げられる。【０００３】
自動車から取得されるセンサデータ、健康診断の検査値履歴、電力需要履歴などに代表されるデータは、すべて「系列依存性（本例では時間的な依存性）」を持つ多変量データである。このようなデータの分析は、産業上重要な数多くの分野に適用される。例えば、自動車から取得されるセンサデータの分析によって、自動車の故障原因を解析して素早い修理を実現することが考えられる。また、健康診断の検査値履歴の分析によって、疾患のリスクの推定および疾患の予防を実現することが考えられる。また、電力需要履歴の分析によって、電力の需要を予測して過不足に備えられるようにすることが考えられる。【０００４】
このようなデータに対しては、系列依存性を持つ隠れ変数モデル（例えば、隠れマルコフモデル）を用いてモデル化する事が一般的である。例えば、隠れマルコフモデルを利用するためには、隠れ状態数、観測確率分布の種類、及び分布パラメータを決定する必要がある。隠れ状態数と観測確率分布の種類がわかっている場合には、Expectation Maximization法（例えば非特許文献１参照）を利用してパラメータを推定する事が可能である。【０００５】
隠れ状態数や観測確率の種類を決定する問題は、一般的に「モデル選択問題」や「システム同定問題」と呼ばれ、信頼性のあるモデルを構築するために極めて重要な問題である。そのための技術が種々提案されている。【０００６】
例えば、非特許文献２では、隠れ状態数を決定する方法として、変分ベイズ法によって変分自由エネルギーを最大化する方法が提案されている。また、例えば、非特許文献３では、隠れ状態数を決定する方法として、階層Dirichlet過程事前分布を用いたノンパラメトリックベイズ法が提案されている。【０００７】
また、非特許文献４では、時間依存性のない隠れ変数モデルの代表例である混合モデルに対して、完全周辺尤度関数を近似し、その下界を最大化することが述べられている。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】C. Bishop, Pattern Recognition and Machine Learning, Springer, 2007, pp.610-629
【非特許文献２】Beal, M. J. Variational Algorithms for Approximate Bayesian Inference. Chapter3, PhD thesis, University College London, 2003
【非特許文献３】van Gael, J., Saatci, Y, Teh, Y.-W., and Ghahramani, Z. Beam sampling for the infinite hidden Markov model. In ICML, 2008.
【非特許文献４】Ryohei Fujimaki, Satoshi Morinaga: Factorized Asymptotic Bayesian Inference for Mixture Modeling. Proceedings of the the fifteenth international conference on Artificial Intelligence and Statistics (AISTATS), 2012
【発明の概要】
【発明が解決しようとする課題】
【０００９】
非特許文献２に記載の方法では、周辺化尤度関数の下界を最大化する際に、変分分布上における隠れ状態と分布パラメータの独立性を仮定するため、周辺化尤度の近似精度が悪くなるという問題がある。【００１０】
非特許文献３に記載の方法では、モンテカルロ法に基づく最適化アルゴリズムが知られているが、計算量が非常に多くなるという問題がある。【００１１】
非特許文献２に記載の方法と非特許文献３に記載の方法で観測確率の種類を決定することは、計算量が極めて多くなるという問題のため、事実上困難である。【００１２】
この計算量の問題を、観測確率分布が混合多項式曲線となる場合を例に説明する。なお、隠れ状態については、以下の議論には影響しないため省略する。ある隠れ状態に対する観測が多項式曲線となる場合には、１次曲線（直線）、２次曲線、３次曲線など曲線の次数を正しく選ぶ必要がある。上記の方法では、例えば、隠れ状態数が３で直線と２次曲線が２つの場合、隠れ状態数が５で３次曲線が３つと４次曲線が２つの場合等の全てのモデルの候補に対して情報量基準を計算する必要があった。このモデルの候補の数は、例えば隠れ状態数が１０で、曲線の最大次数を１０とすると約十万通りとなり、隠れ状態数が２０で曲線の最大次数を２０とすると数百億通りとなる。このようにモデルの候補の数は、探索すべきモデルの複雑さと共に指数的に増加する。従って、非特許文献２および非特許文献３に記載の方法による計算は、事実上困難であった。【００１３】
また、非特許文献４に記載の技術では、隠れ変数間の独立性が必要となり、系列依存性のある隠れ変数モデルに適用できなかった。非特許文献４に記載の技術では、隠れ変数間の系列依存性を考慮していないため、隠れ変数の変分分布が、非特許文献４の式（１５）として算出される。しかし、隠れ変数間に系列依存性がある場合には、この式は適切でなく、適切なモデルが得られる保証はない。さらに、隠れ変数間の遷移確率が算出できないという問題がある。【００１４】
そこで、本発明は、多変量データに対する系列依存性を持つ隠れ変数モデルの学習問題において、隠れ状態数および観測確率の種類と共にモデルの候補数が指数的に増加しても高速にモデル選択を実現できるようにすることを目的とする。
【課題を解決するための手段】
【００１５】
本発明による隠れ変数モデル推定装置は、周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値を最大化することによって変分確率を計算する変分確率計算部と、各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定するモデル推定部と、変分確率計算部が変分確率を計算する際に用いた基準値が収束したか否かを判定する収束判定部とを備えることを特徴とする。【００１６】
また、本発明による隠れ変数モデル推定方法は、周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値を最大化することによって変分確率を計算し、各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定し、変分確率を計算する際に用いた基準値が収束したか否かを判定することを特徴とする。【００１７】
また、本発明による隠れ変数モデル推定プログラムは、コンピュータに、周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値を最大化することによって変分確率を計算する変分確率計算処理、各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定するモデル推定処理、および、変分確率計算処理で変分確率を計算する際に用いた基準値が収束したか否かを判定する収束判定処理を実行させる。
【発明の効果】
【００１８】
本発明によれば、隠れ状態数および観測確率の種類と共にモデルの候補数が指数的に増加しても高速にモデル選択を行える。
【図面の簡単な説明】
【００１９】
【図１】本発明の隠れ変数モデル推定装置の第１の実施形態を示すブロック図である。
【図２】隠れ変数変分確率計算処理部１０４の例を示すブロック図である。
【図３】本発明の第１の実施形態の処理経過の例を示すフローチャートである。
【図４】隠れ変数変分確率計算処理部１０４の動作を示すフローチャートである。
【図５】本発明の隠れ変数モデル推定装置の第２の実施形態を示すブロック図である。
【図６】本発明の隠れ変数モデル推定装置の概要を示すブロック図である。
【発明を実施するための形態】
【００２０】
以下、本発明の実施形態を図面を参照して説明する。なお、以下の説明では便宜的に、数式内の表記と、文章中の表記とが異なる場合がある。例えば、記号“~”を数式内では、変数の上部に記載するが、文章中では便宜的に右側に記載する等の違いがある。このような数式内の表記と文章中の表記との相違は、当業者が理解し得る範囲内の相違である。
【００２１】
本発明の隠れ変数モデル推定装置は、系列依存性を持つ隠れ変数モデルを推定する。以下の説明では、系列依存性を持つデータの例として時間依存性を持つデータを例にして説明するが、本発明におけるデータは系列依存性を持つデータであればよく、時間依存性を持つデータに限定されない。例えば、文字系列に依存するデータ、塩基系列に依存するデータ、あるいはその他の系列に依存するデータであってもよい。
【００２２】
また、以下の説明では、系列依存性を持つ隠れ変数モデルの最も代表的なモデルである隠れマルコフモデル（式（１）参照）を例にして具体的な説明を行う。
【００２３】
【数１】
2013179579.tif 000003
【００２４】
なお、以下の説明では時間依存するデータ列xn(n=1,...,N)が入力されると仮定する。ここで、各xnは、長さTnの多変量データ列（xn=xn1, ..., xnT，t=1,...,N)であるとする。次に、観測変数xntに対する隠れ変数znt=(znt1, ..., zntK)を定義する。zntk=1は、xntがk番目の隠れ状態から生成されたデータである事を意味し、zntk=0は、そうでない事を意味する。また、Σk=1K zntk = 1である。xとzの組は「完全変数」と呼ばれる。なお、その対比としてxを不完全変数と呼ぶ。完全変数に関する隠れマルコフモデルの同時分布は、式（１）中のP（x、z）として定義される。なお、隠れマルコフモデルの隠れ変数に対する変分分布は、時刻tにおけるk番目の隠れ状態zntkの分布q(zntk)、及び、時刻t-1から時刻tにおいてk番目の状態からj番目の状態へ遷移する分布q(znt-1k, zntj)と表される。
【００２５】
式（１）において、Kは隠れ状態数を表す。また、θ＝（α1, ...,αK, β1, ..., βK,φ1, ..., φK）は、モデルのパラメータを表す。ここで、αkはk番目の隠れ状態の初期確率を表し、βkはk番目の隠れ状態からの遷移確率を表し、φkはk番目の隠れ状態に対する観測パラメータを表す。また、S1, ..., SKは、φkに対応する観測確率の種類を表す。なお、S1からSKとなりうる候補は、例えば、多変量データの生成確率の場合には｛正規分布、対数正規分布、指数分布｝であったり、多項曲線出力の場合では、｛０次曲線、１次曲線、２次曲線、３次曲線｝であったりする。
【００２６】
なお、本明細書では具体的な例はすべて隠れマルコフを用いて説明をするが、その拡張モデル（例えば隠れセミマルコフモデル、因子化隠れマルコフモデルなど）の類似のモデルにも本発明を適用可能である。同様に、本明細書では、ターゲット変数をＸとした分布について説明しているが、回帰や判別のように、観測分布が条件付モデルＰ（Ｙ｜Ｘ）（Ｙはターゲットとなる確率変数）である場合に関しても本発明を適用可能である。
【００２７】
実施形態１．
図１は、本発明の隠れ変数モデル推定装置の第１の実施形態を示すブロック図である。隠れ変数モデル推定装置１００は、データ入力装置１０１と、隠れ状態数設定部１０２と、初期化処理部１０３と、隠れ変数変分確率計算処理部１０４と、モデル最適化処理部１０５と、最適性判定処理部１０６と、最適モデル選択処理部１０７と、モデル推定結果出力装置１０８とを備える。隠れ変数モデル推定装置１００には、入力データ１１１が入力され、入力データ１１１に対して隠れ状態数および観測確率の種類を最適化し、モデル推定結果１１２として出力する。
【００２８】
また、図２は、隠れ変数変分確率計算処理部１０４の例を示すブロック図である。隠れ変数変分確率計算処理部１０４は、前向き確率計算処理部１０４１と、正規化定数記憶部１０４２と、後ろ向き確率計算処理部１０４３と、前向き後ろ向き確率合算処理部１０４４を備える。隠れ変数変分確率計算処理部１０４には、入力データ１１１と、モデル最適化処理部１０５で推定された推定モデル１０４５とが入力され、隠れ変数変分確率１０４６と、前向き確率正規化定数１０４７とを出力する。
【００２９】
入力装置１０１は、入力データ１１１を入力するための入力インタフェース装置である。入力装置１０１には、入力データ１１１が入力される際に、観測確率の種類や、隠れ状態数の候補値など、モデルの推定に必要なパラメータも同時に入力される。
【００３０】
隠れ状態数設定部１０２は、入力された隠れ状態数の候補値からモデルの隠れ状態数を選択して設定する。以下、設定された隠れ状態数をＫと表記する。
【００３１】
初期化処理部１０３は、推定のための初期化処理を実施する。なお、初期化は任意の方法によって実施することが可能である。例としては、観測確率の種類を隠れ状態ごとにランダムに設定し、設定された種類に従って、各観測確率のパラメータをランダムに設定する方法や、隠れ変数の変分確率をランダムに設定する方法が挙げられる。
【００３２】
隠れ変数変分確率計算処理部１０４は、隠れ変数の変分確率を計算する。ここで、パラメータθは初期化処理部１０３あるいはモデル最適化処理部１０６で計算されているため、隠れ変数変分確率計算処理部１０４は、その値を利用する。隠れ変数変分確率計算処理部１０４は、次に定義する最適化基準Ａを最大化することによって変分確率を計算する。最適化基準Ａとは、周辺化対数尤度関数を完全変数に対する推定量（例えば最尤推定量や最大事後確率推定量）に関してラプラス近似した近似量の下界として定義される。なお、この下界は、完全変数に対する推定量の最適性と対数関数の凹性を用いることで導出することが可能である。
【００３３】
この手順を、隠れマルコフモデルを例に説明する。まず、周辺化対数尤度関数の下界を考える。この下界は、以下の式（２）で示される
【００３４】
【数２】
2013179579.tif 000004
【００３５】
なお、式（２）において、変分確率q（zN)を最大化する事で等号が成立する。ここで、完全変数に対する最尤推定量を用いて分子の完全変数の周辺化尤度をラプラス近似する事で、周辺化対数尤度関数の近似式として、以下に示す式（３）を得る。
【００３６】
【数３】
2013179579.tif 000005
【００３７】
ただし、上付きのバーは完全変数に対する最尤推定量を表す。また、D＊は下付きのパラメータ＊の次元を表す。
【００３８】
次に、式（３）に対して最尤推定量が対数尤度関数を最大化する性質と、対数関数が凹関数である事を利用して、式（３）の下界を以下の式（４）のように算出する。
【００３９】
【数４】
2013179579.tif 000006
【００４０】
隠れ変数の変分分布q(zntk)及びq(znt-1k, zntj)は、式（４）をqについて最大化する事によって算出される。ただし、上付き(i)を、隠れ変数変分確率計算処理部１０４、モデル最適化処理部１０５、最適判定処理部１０６の繰り返し計算における、(i)回目の繰り返しを表すとすると、q(i)はq~= q(i-1)、θ＝θ(i-1)と固定する。
【００４１】
なお、式（４）において下線を付した部分をＢとする。Ｂは、後述の式（８）で参照する。
【００４２】
図２を参照して、隠れ変数変分確率計算処理部１０４が備える要素について説明する。前向き確率計算処理部１０４１には、入力データ１１１と推定モデルが入力される。そして、前向き確率計算処理部１０４１は、時刻1から時刻tまでの観測（xn1, ..., xnt）が得られた場合のzntの確率を前向き確率として算出する。ただし、前向き確率は、最適化基準Aで算出されるモデル複雑度（例えば式（４）ではδtkに関する項）を考慮して算出される。また、前向き確率計算処理部１０４１は、zntの確率の隠れ状態に関する和を１とするための正規化定数を、正規化定数記憶部１０４２に記憶させる。
【００４３】
同様に、後ろ向き確率計算処理部１０４３は、時刻t+1からTまでの観測（xnt+1, ..., xnT）が得られた場合のxntの確率を後ろ向き確率として算出する。なお、後ろ向き確率の計算の際に、前向き確率算出と同時に得られる正規化定数を正規化定数記憶部１０４２から読み込む。ただし、後ろ向き確率は最適化基準Aで算出されるモデル複雑度（例えば式（４）ではδtkに関する項）を考慮して算出される。
【００４４】
最後に、前向き後ろ向き確率合算処理部１０４４は、前向き確率と後ろ向き確率から、変分分布を算出する。例えば、前向き後ろ向き確率合算処理部１０４４は、q(zntk)をxn1, ..., xnTが得られたときのzntkの確率として計算する。前向き後ろ向き確率合算処理部１０４４は、前向き確率と後ろ向き確率の積として、以下の式（５）の計算によって、q(zntk)を計算する。
【００４５】
【数５】
2013179579.tif 000007
【００４６】
また、前向き後ろ向き確率合算処理部１０４４は、xn1, ..., xnt-1が得られたときのznt-1jの確率と、隠れ状態jから隠れ状態kへ遷移する確率と、隠れ状態kにおいてxntが観測される確率と、（xnt+1, ..., xnT）が得られた場合のxntの確率の積として、q(znt-1j, zntk)を算出する。具体的には、前向き後ろ向き確率合算処理部１０４４は、以下の式（６）の計算によってq(znt-1j, zntk)を算出する（式（６）左辺のq~の定義は式（７）を参照）。
【００４７】
【数６】
2013179579.tif 000008
【００４８】
この手順を隠れマルコフモデルを例に説明すると、前向き確率及び後ろ向き確率は、以下の式（７）の計算によって算出される。
【００４９】
【数７】
2013179579.tif 000009
【００５０】
ただし、ftnk（式（７）の第１式）が前向き確率を表し、btnk（式（７）の第２式）が後ろ向き確率を表す。より具体的には、式（７）において、前向き確率と後ろ向き確率の両者とも、漸化式として記述されている。そして、前向き確率はt=1から順に算出することが可能であり、後ろ向き確率はt=Tから順に算出することが可能である。なお、正規化定数はζtnで算出される。後ろ向き確率計算処理部１０４３は、前向き確率計算処理部１０４１が前向き確率を算出する際に計算した正規化定数を利用して後ろ向き確率を算出すればよい。
【００５１】
また、式（５）の第３式において、δに関する乗算が含まれているが、これは、最適化基準Aで算出されるモデル複雑度を考慮していることを意味する。
【００５２】
モデル最適化処理部１０５は、式（４）に対してモデル（パラメータθ及びその種類S） を最適化する。具体的には、q及びq~を隠れ変数変分確率計算処理部１０４で算出された隠れ変数の変分分布（q(i)）に固定し、式（４）におけるGを最大化するモデルを算出する。この処理において重要な点は、式（４）によって定義されたGは、コンポーネントごとに最適化の関数を分解する事が可能なため、コンポーネント種類の組合せ(S1からSKのどの種類を指定するか)を考慮することなく、S1からSK及びパラメータφ1からφKを別々に最適化する事が可能な点である。これによって、コンポーネントの種類を最適化する際に、組み合わせ爆発を回避して最適化を実行する事が可能となる。
【００５３】
最適性判定処理部１０６は、式（４）で計算される最適化基準Aの収束を判定する。最適化基準Aが収束したと判定していない場合には、隠れ変数変分確率計算処理部１０４から最適性判定処理部１０６の処理を繰り返す。なお、式（４）で計算される最適化基準Aの算出は、隠れ状態が独立ではないためΣzn q(zn) log q(zn)の計算に指数時間の計算量を必要とするが、正規化定数記憶部１０４２に記憶されている正規化定数を利用して効率的に計算をする事が可能である。例えば隠れマルコフモデルの場合には、以下の式（８）のように計算される。
【００５４】
【数８】
2013179579.tif 000010
【００５５】
式（８）に示すＢは、式（４）において下線を付した部分である。
【００５６】
隠れ変数変分確率計算処理部１０４から最適性判定処理部１０６の処理を繰り返し、変分分布とモデルを更新する事で、適切なモデルを選択する事が可能となる。なお、この繰り返しによって最適化基準Aが単調に増加する事が保証される。
【００５７】
最適化基準Aが収束している場合、隠れ状態数設定部１０２で設定された隠れ状態数Kに対して、隠れ変数変分確率計算処理部１０４から最適性判定処理部１０６のループ処理で算出される最適化基準Aと、その１つ前のループ処理で算出された最適化基準Aのうち、大きい方の最適化基準Aに対応するモデルを最適なモデルとして設定する。全ての候補値についてモデルの最適化が完了した場合には、処理がモデル推定結果出力装置１０８へ移り、まだ最適化の済んでいない候補が存在する場合には、隠れ状態数設定部１０２へ処理が移る。
【００５８】
モデル推定結果出力装置１０８は、最適な隠れ状態数、観測確率の種類、パラメータ、変分分布などをモデル推定結果出力結果１１２として出力する。
【００５９】
隠れ状態数設定部１０２、初期化処理部１０３、隠れ変数変分確率計算処理部１０４（前向き確率計算処理部１０４１、正規化定数記憶部１０４２、後ろ向き確率計算処理部１０４３、前向き後ろ向き確率合算処理部１０４４）、モデル最適化処理部１０５、最適性判定処理部１０６、最適モデル選択処理部１０７およびモデル推定結果出力装置１０８は、例えば、隠れ変数モデル推定プログラムに従って動作するコンピュータのＣＰＵによって実現される。ＣＰＵが、隠れ変数モデル推定プログラムを記録したコンピュータ読み取り可能な記録媒体から隠れ変数モデル推定プログラムを読み込み、上記の各要素として動作すればよい。
【００６０】
また、隠れ状態数設定部１０２、初期化処理部１０３、隠れ変数変分確率計算処理部１０４、モデル最適化処理部１０５、最適性判定処理部１０６、最適モデル選択処理部１０７およびモデル推定結果出力装置１０８が別々のハードウェアで実現されていてもよい。また、隠れ変数変分確率計算処理部１０４においても、前向き確率計算処理部１０４１、正規化定数記憶部１０４２、後ろ向き確率計算処理部１０４３、前向き後ろ向き確率合算処理部１０４４が別々のハードウェアで実現されていてもよい。
【００６１】
図３は、本発明の第１の実施形態の処理経過の例を示すフローチャートである。データ入力装置１０１を介して入力データ１１１が入力される（ステップＳ１００）。
【００６２】
次に、隠れ状態数設定部１０２は、入力された隠れ状態数の候補値のうち、まだ最適化の行なわれていない隠れ状態数を選択し設定する（ステップＳ１０１）。
【００６３】
次に、初期化処理部１０３は、設定された隠れ状態数に対して、推定のため、パラメータや隠れ変数変分確率の初期化処理を実施する（ステップＳ１０２）。
【００６４】
次に、隠れ変数変分確率計算処理部１０４は、隠れ変数の変分確率を計算する（ステップＳ１０３）。
【００６５】
次に、モデル最適化処理部１０５は、各隠れ状態に対して観測確率の種類とパラメータの推定を実施する（ステップＳ１０４）。この処理は、各隠れ状態のモデルの最適化であると言うことができる。
【００６６】
次に、最適性判定処理部１０６は、最適化基準A が収束したかを判定する。（例えば、Ｓ１０５）。最適性判定処理部１０６は、直近のステップＳ１０３〜Ｓ１０５のループ処理において得られた最適化基準Aと、その１つ前のステップＳ１０３〜Ｓ１０５のループ処理で得られた最適化基準Aの差を計算し、その差の絶対値が予め定められた閾値以下になっていれば、最適化基準A が収束したと判定してよい。また、その差の絶対値が閾値より大きければ、最適性判定処理部１０６は、最適化基準A が収束していないと判定してよい。なお、絶対値による最適化基準Aの差の算出は一例であり、例えば相対的な差によって収束を判定する等の方法を採用してもよい。
【００６７】
ステップＳ１０５において、最適化基準A が収束していないと判定された場合には、ステップＳ１０３〜Ｓ１０５の処理を繰り返す。
【００６８】
ステップＳ１０５において、最適化基準A が収束したと判定された場合、最適モデル選択処理部１０７は、直近のステップＳ１０３〜Ｓ１０５のループ処理において最適化されたモデル（隠れ状態数、観測確率の種類、パラメータ）の最適化基準Aと、その１つ前のステップＳ１０３〜Ｓ１０５のループ処理で最適化されたモデルの最適化基準Aとを比較し、値の大きい方の最適化基準Aに対応するモデルを、最適なモデルとして設定する（ステップＳ１０６）。
【００６９】
次に、隠れ状態数設定部１０２は、推定されていない隠れ状態数の候補が残っているか否かを判定する（ステップＳ１０７）。隠れ状態数の候補が残っている場合には、ステップＳ１０２〜Ｓ１０７の処理を繰り返す。一方、隠れ状態数の候補が残っていない場合には、モデル推定結果出力装置１０８がモデル推定結果を出力し（ステップＳ１０８）、処理を完了する。
【００７０】
図４は、本実施形態における隠れ変数変分確率計算処理部１０４の動作（換言すれば、ステップＳ１０３の処理経過）を示すフローチャートである。
【００７１】
前向き確率計算処理部１０４１は、n番目のデータのt番目の時刻に対する前向き確率ft(i)nkを算出する（ステップＳ１１１）。このとき、前向き確率計算処理部１０４１は、正規化定数も算出し、正規化定数記憶部１０４２に記憶させる（ステップＳ１１２）。
【００７２】
続いて、前向き確率計算処理部１０４１は、すべての時刻tに対して前向き確率の算出が完了したかを確認し（ステップＳ１１３）、未完了の場合にはステップＳ１１１，Ｓ１１２の処理を繰り返す。完了した場合には、ステップＳ１１４の処理に移る。
【００７３】
後ろ向き確率計算処理部１０４３は、n番目のデータのt番目の時刻に対する後ろ向き確率bt(i)nkを算出する（ステップＳ１１４）。そして、後ろ向き確率計算処理部１０４３は、すべての時刻tに対して後ろ向き確率の算出が完了したかを確認し（ステップＳ１１５）、未完了の場合にはステップＳ１１４の処理を繰り返す。完了した場合には、ステップＳ１１６の処理に移る。
【００７４】
前向き後ろ向き確率合算処理部１０４４は、n番目のデータのすべての時刻に対して前向き確率、後ろ向き確率の合算処理を行い、変分分布の計算を行う（ステップＳ１１６）。
【００７５】
続いて、前向き後ろ向き確率合算処理部１０４４は、nに関してすべてのデータに対して変分分布の算出処理が完了しているかを確認する（ステップＳ１１７）。未完了の場合には、ステップＳ１１１以降の処理を繰り返し、完了した場合には、処理を終了する。
【００７６】
隠れ状態数および観測確率の種類と共にモデルの候補数が指数的に増加する場合であっても、上記のような本発明の動作（特に、隠れ変数変分確率計算処理部１０４の動作）によって、高速にモデル選択を実現することができる。
【００７７】
また、非特許文献４に記載の技術と本発明を比較すると、既に説明したように、非特許文献４に記載の技術では、隠れ変数間の独立性が必要となり、系列依存性のある隠れ変数モデルに適用できなかった。これに対し本発明は、系列依存性を持つ多変量データの隠れ変数モデルを推定することができる。
【００７８】
実施形態２．
図５は、本発明の隠れ変数モデル推定装置の第２の実施形態を示すブロック図である。第２の実施形態の隠れ変数モデル推定装置２００は、第１の実施形態の隠れ変数モデル推定装置１００（図１参照）と比較して、最適モデル選択処理部１０７を備えず、隠れ状態数選択処理部２０１を備える。
【００７９】
データ入力装置１０１、隠れ状態数設定部１０２、初期化処理部１０３、隠れ変数変分確率計算処理部１０４、モデル最適化処理部１０５、最適性判定処理部１０６およびモデル推定結果出力装置１０８に関しては、第１の実施形態と同様である。
【００８０】
また、第１の実施形態の隠れ変数モデル推定装置１００は、隠れ状態数の候補に対してモデル最適化を行い、最適化基準Aを最大化するモデルを選択する。これに対して、第２の実施形態の隠れ変数モデル推定装置２００では、隠れ変数変分確率計算処理部１０４の処理の後で、隠れ状態数選択処理部２０１が、小さくなった隠れ状態をモデルから除去する。
【００８１】
具体的には、隠れ状態数選択処理部２０１は、隠れ変数変分確率計算処理部１０４で算出されたq(zntk)に対して、以下の式（９）の状態を満たす隠れ状態を除去する。
【００８２】
【数９】
2013179579.tif 000011
【００８３】
式（９）の右辺に示すεは、入力データ１１１と同時に入力される閾値である。すなわち、隠れ状態数選択処理部２０１は、閾値ε以下である隠れ状態を除去する。
【００８４】
式（９）による隠れ状態の除去は、以下の理由によって正しく動く事が説明される。まず、式（７）の前向き確率を観察すると、小さい隠れ状態（すなわち小さいδtkに対応する隠れ状態）に対する前向き確率は小さくなる。また、後ろ向き確率では、小さい隠れ状態は前の状態に対する寄与が小さい。従って、前向き確率と後ろ向き確率から算出される変分分布は、小さい隠れ状態に対する確率が繰り返し最適化を通して徐々に小さくなっていく（前の更新ステップにおいて小さい隠れ状態ほど、次の更新ステップでも小さくなりやすくなるため）。このような構成とする事で、隠れ変数モデル推定装置１００のように複数の隠れ状態数の候補に対して最適化をする必要がなく、隠れ状態数、観測確率の種類とパラメータ、変分分布を同時に推定し、計算コストを抑える事が可能という利点がある。
【００８５】
第２の実施形態において、隠れ状態数設定部１０２、初期化処理部１０３、隠れ変数変分確率計算処理部１０４、隠れ状態数選択処理部２０１、モデル最適化処理部１０５、最適性判定処理部１０６およびモデル推定結果出力装置１０８は、例えば、隠れ変数モデル推定プログラムに従って動作するコンピュータのＣＰＵによって実現される。ＣＰＵが、隠れ変数モデル推定プログラムを記録したコンピュータ読み取り可能な記録媒体から隠れ変数モデル推定プログラムを読み込み、上記の各要素として動作すればよい。また、第２の実施形態におけるこれらの各要素が別々のハードウェアで実現されていてもよい。
【実施例１】
【００８６】
以下、本発明の第１の実施形態の応用例を、自動車のセンサデータに対する走行モード分析を例にして説明する。以下の例では理解を容易にするため１次元の例を説明するが、多次元でも同様に適用可能である。【００８７】
第１の実施形態の隠れ変数モデル推定装置を利用することで、自動車に設置された複数のセンサから取得される多次元時系列データに対して、例えば「走行モード」のように、時系列を、複数の異なる性質へ分解することができる。センサデータからの故障診断や異常挙動の検出を考えた場合に、走行モードによってセンサの振る舞いは大きく異なる。そのため、モードへ分解し分析する事が必要であり、これを自動化する事は非常に重要である。【００８８】
例えば、エンジン回転数をXとし、速度をYとした、多項回帰出力の隠れマルコフモデルを考える。このとき、推定すべきモデルは、隠れ状態数、各隠れ状態に対する回帰次数（Sk）、回帰パラメータ（φk）、初期確率（αk）、遷移確率（βk）、変分分布（q）である。【００８９】
まず、エンジン回転数と速度の時系列データとともに、隠れ状態数の候補値としてK=1から10までを隠れ変数モデル推定装置１００に入力する。隠れ状態数設定部１０２は、K=1から10まで順に隠れ状態数を設定する。次に、初期化処理部１０３は、初期化処理として、K個の隠れ状態に対して、回帰次数及びその他のパラメータをランダムに設定する。次に、隠れ変数変分確率計算処理部１０４から最適性判定処理部１０６によってモデルの推定を行う。この処理を通じて、例えば、エンジン回転数が一定で速度が増加している状態（等加速）に対応するXからYへの０次多項式、エンジン回転数も速度も減少している状態（減速中）に対応するXからYへの１次多項式、エンジン回転数が急激に増加し速度が徐々に増加する状態（急加速）に対応するXからYへの2次多項式など、異なる走行状態が、異なる次数と係数を持つ回帰モデルとして自動的に分離される。さらに、最適モデル選択処理部１０７が、最もよい隠れ状態数を自動選択するため、例えばドライバーに応じて異なる運転特性（モード）の数を、自動的に検出し、適切な数の走行モードに分離する事が可能である。
【実施例２】
【００９０】
以下、本発明の第２の実施形態の応用例を、診療ログ（レセプトデータ）からの疾病パタン分析を例にして説明する。例えば、心筋梗塞などは高血圧や糖尿病といった生活習慣病を事前に併発している事が多く、また生活習慣病は一度治ったとしても、再発する事が多い。そのような疾病のパタンを分析する事で、疾病のリスク低減への方策検討や、生活習慣指導へ活用する事が可能である。【００９１】
本例では、高血圧にかかっているかどうかの論理値（かかっていれば１、かかっていなければ０）を複数の疾病について並べた、多次元の論理値ベクトル時系列を入力データとする。推定するモデルは、多次元のベルヌーイ観測の隠れマルコフモデルとする。【００９２】
まず、入力データとともに、隠れ状態数としてKmax、選択のしきい値εを入力する。隠れ状態の候補値はKmaxに設定され、ベルヌーイ分布のパラメータがランダムに初期化される。次に、隠れ変数変分確率計算処理部１０４から最適性判定処理部１０６によってモデルの推定を行う。この処理を通じて、高血圧と糖尿病が合併しているパタン、高脂血症が治ったり再発したりを繰り返すパタン（薬による治療中）、生活習慣病がほとんど発生しないパタンなどに分離するとともに、典型的ではないパタンに対応する隠れ状態は小さくなり、隠れ状態数選択処理部２０１によって除去され、最終的な推定結果として典型的なパタンのみを抽出する事が可能である。【００９３】
図６は、本発明の隠れ変数モデル推定装置の概要を示すブロック図である。本発明の隠れ変数モデル推定装置は、変分確率計算部７１と、モデル推定部７２と、収束判定部７３とを備える。【００９４】
変分確率計算部７１（例えば、隠れ変数変分確率計算処理部１０４）は、周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値（例えば、最適化基準Ａ）を最大化することによって変分確率を計算する。【００９５】
モデル推定部７２（例えば、モデル最適化処理部１０５）は、各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定する。【００９６】
収束判定部７３（例えば、最適化判定処理部１０６）は、変分確率計算部７１が変分確率を計算する際に用いた基準値が収束したか否かを判定する。【００９７】
また、変分確率計算部７１が変分確率を計算し、モデル推定部７２が最適な隠れ変数モデルを推定し、収束判定部７３が、基準値が収束したか否かを判定するループ処理を繰り返し、基準値が収束したときに、当該基準値と、その１つ前のループ処理における基準値とを比較し、大きい方の基準値に対応する隠れ変数モデルを最適な隠れ変数モデルとして選択する最適モデル選択部（例えば、最適モデル選択処理部１０７）を備える構成であってもよい。【００９８】
また、変分確率計算部７１が変分確率を計算し、モデル推定部７２が最適な隠れ変数モデルを推定し、収束判定部７３が、基準値が収束したか否かを判定するループ処理を繰り返し、変分確率計算部の計算結果に応じて、所定の条件を満たす隠れ状態を除去する隠れ状態除去部（例えば、隠れ状態数選択処理部２０１）を備える構成であってもよい。【００９９】
また、モデル推定部７２が、隠れ変数モデルとして隠れマルコフモデルを推定する構成であってもよい。【０１００】
この出願は、２０１２年５月３１に出願された米国仮出願61/653,855を基礎とする優先権を主張し、その開示の全てをここに取り込む。【０１０１】
以上、実施形態を参照して本願発明を説明したが、本願発明は上記の実施形態に限定されるものではない。本願発明の構成や詳細には、本願発明のスコープ内で当業者が理解し得る様々な変更をすることができる。
【産業上の利用の可能性】
【０１０２】
本発明は、系列依存性を持つ多変量データの隠れ変数モデル推定装置に好適に適用される。
【符号の説明】
【０１０３】
１０１  データ入力装置
１０２  隠れ状態数設定部
１０３  初期化処理部
１０４  隠れ変数変分確率計算処理部
１０５  モデル最適化処理部
１０６  最適性判定処理部
１０７  最適モデル選択処理部
１０８  モデル推定結果出力装置
２０１  隠れ状態数選択処理部
【図１】
2013179579.tif 000012
【図２】
2013179579.tif 000013
【図３】
2013179579.tif 000014
【図４】
2013179579.tif 000015
【図５】
2013179579.tif 000016
【図６】
2013179579.tif 000017
【国際調査報告】
2013179579.tif 000018
2013179579.tif 000019
2013179579.tif 000020
2013179579.tif 000021