JP6020561B2 - 隠れ変数モデル推定装置および方法 - Google Patents

Description

本発明は、系列依存性を持つ多変量データの隠れ変数モデル推定装置、隠れ変数モデル推定方法、および隠れ変数モデル推定プログラムに関し、特に、モデル事後確率を近似し、その下界を最大化する事によって系列依存性を持つ多変量データの隠れ変数モデルを推定する隠れ変数モデル推定装置、隠れ変数モデル推定方法、および隠れ変数モデル推定プログラムに関する。

系列依存性を持つデータが種々存在する。系列依存性を持つデータの例として、例えば、時間依存性を有するデータや、文字系列に依存する文章、塩基系列に依存する遺伝子データ等が挙げられる。

自動車から取得されるセンサデータ、健康診断の検査値履歴、電力需要履歴などに代表されるデータは、すべて「系列依存性（本例では時間的な依存性）」を持つ多変量データである。このようなデータの分析は、産業上重要な数多くの分野に適用される。例えば、自動車から取得されるセンサデータの分析によって、自動車の故障原因を解析して素早い修理を実現することが考えられる。また、健康診断の検査値履歴の分析によって、疾患のリスクの推定および疾患の予防を実現することが考えられる。また、電力需要履歴の分析によって、電力の需要を予測して過不足に備えられるようにすることが考えられる。

このようなデータに対しては、系列依存性を持つ隠れ変数モデル（例えば、隠れマルコフモデル）を用いてモデル化する事が一般的である。例えば、隠れマルコフモデルを利用するためには、隠れ状態数、観測確率分布の種類、及び分布パラメータを決定する必要がある。隠れ状態数と観測確率分布の種類がわかっている場合には、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

非特許文献２に記載の方法では、周辺化尤度関数の下界を最大化する際に、変分分布上における隠れ状態と分布パラメータの独立性を仮定するため、周辺化尤度の近似精度が悪くなるという問題がある。

非特許文献３に記載の方法では、モンテカルロ法に基づく最適化アルゴリズムが知られているが、計算量が非常に多くなるという問題がある。

非特許文献２に記載の方法と非特許文献３に記載の方法で観測確率の種類を決定することは、計算量が極めて多くなるという問題のため、事実上困難である。

この計算量の問題を、観測確率分布が混合多項式曲線となる場合を例に説明する。なお、隠れ状態については、以下の議論には影響しないため省略する。ある隠れ状態に対する観測が多項式曲線となる場合には、１次曲線（直線）、２次曲線、３次曲線など曲線の次数を正しく選ぶ必要がある。上記の方法では、例えば、隠れ状態数が３で直線と２次曲線が２つの場合、隠れ状態数が５で３次曲線が３つと４次曲線が２つの場合等の全てのモデルの候補に対して情報量基準を計算する必要があった。このモデルの候補の数は、例えば隠れ状態数が１０で、曲線の最大次数を１０とすると約十万通りとなり、隠れ状態数が２０で曲線の最大次数を２０とすると数百億通りとなる。このようにモデルの候補の数は、探索すべきモデルの複雑さと共に指数的に増加する。従って、非特許文献２および非特許文献３に記載の方法による計算は、事実上困難であった。

また、非特許文献４に記載の技術では、隠れ変数間の独立性が必要となり、系列依存性のある隠れ変数モデルに適用できなかった。非特許文献４に記載の技術では、隠れ変数間の系列依存性を考慮していないため、隠れ変数の変分分布が、非特許文献４の式（１５）として算出される。しかし、隠れ変数間に系列依存性がある場合には、この式は適切でなく、適切なモデルが得られる保証はない。さらに、隠れ変数間の遷移確率が算出できないという問題がある。

そこで、本発明は、多変量データに対する系列依存性を持つ隠れ変数モデルの学習問題において、隠れ状態数および観測確率の種類と共にモデルの候補数が指数的に増加しても高速にモデル選択を実現できるようにすることを目的とする。

本発明による隠れ変数モデル推定装置は、周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値を最大化することによって変分確率を計算する変分確率計算部と、各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定するモデル推定部と、変分確率計算部が変分確率を計算する際に用いた基準値が収束したか否かを判定する収束判定部とを備えることを特徴とする。

また、本発明による隠れ変数モデル推定方法は、周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値を最大化することによって変分確率を計算し、各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定し、変分確率を計算する際に用いた基準値が収束したか否かを判定することを特徴とする。

また、本発明による隠れ変数モデル推定プログラムは、コンピュータに、周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値を最大化することによって変分確率を計算する変分確率計算処理、各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定するモデル推定処理、および、変分確率計算処理で変分確率を計算する際に用いた基準値が収束したか否かを判定する収束判定処理を実行させる。

本発明によれば、隠れ状態数および観測確率の種類と共にモデルの候補数が指数的に増加しても高速にモデル選択を行える。

本発明の隠れ変数モデル推定装置の第１の実施形態を示すブロック図である。 隠れ変数変分確率計算処理部１０４の例を示すブロック図である。 本発明の第１の実施形態の処理経過の例を示すフローチャートである。 隠れ変数変分確率計算処理部１０４の動作を示すフローチャートである。 本発明の隠れ変数モデル推定装置の第２の実施形態を示すブロック図である。 本発明の隠れ変数モデル推定装置の概要を示すブロック図である。

以下、本発明の実施形態を図面を参照して説明する。なお、以下の説明では便宜的に、数式内の表記と、文章中の表記とが異なる場合がある。例えば、記号“~”を数式内では、変数の上部に記載するが、文章中では便宜的に右側に記載する等の違いがある。このような数式内の表記と文章中の表記との相違は、当業者が理解し得る範囲内の相違である。

本発明の隠れ変数モデル推定装置は、系列依存性を持つ隠れ変数モデルを推定する。以下の説明では、系列依存性を持つデータの例として時間依存性を持つデータを例にして説明するが、本発明におけるデータは系列依存性を持つデータであればよく、時間依存性を持つデータに限定されない。例えば、文字系列に依存するデータ、塩基系列に依存するデータ、あるいはその他の系列に依存するデータであってもよい。

また、以下の説明では、系列依存性を持つ隠れ変数モデルの最も代表的なモデルである隠れマルコフモデル（式（１）参照）を例にして具体的な説明を行う。

なお、以下の説明では時間依存するデータ列xⁿ(n=1,...,N)が入力されると仮定する。ここで、各xⁿは、長さTnの多変量データ列（xⁿ=xⁿ1, ..., xⁿT，t=1,...,N)であるとする。次に、観測変数xⁿtに対する隠れ変数zⁿt=(zⁿt1, ..., zⁿtK)を定義する。zⁿtk=1は、xⁿtがk番目の隠れ状態から生成されたデータである事を意味し、zⁿtk=0は、そうでない事を意味する。また、Σ_k=1 ^K zⁿtk = 1である。xとzの組は「完全変数」と呼ばれる。なお、その対比としてxを不完全変数と呼ぶ。完全変数に関する隠れマルコフモデルの同時分布は、式（１）中のP（x、z）として定義される。なお、隠れマルコフモデルの隠れ変数に対する変分分布は、時刻tにおけるk番目の隠れ状態zⁿtkの分布q(zⁿtk)、及び、時刻t-1から時刻tにおいてk番目の状態からj番目の状態へ遷移する分布q(zⁿt-1k, zⁿtj)と表される。

式（１）において、Kは隠れ状態数を表す。また、θ＝（α1, ...,αK, β1, ..., βK,φ1, ..., φK）は、モデルのパラメータを表す。ここで、αkはk番目の隠れ状態の初期確率を表し、βkはk番目の隠れ状態からの遷移確率を表し、φkはk番目の隠れ状態に対する観測パラメータを表す。また、S1, ..., SKは、φkに対応する観測確率の種類を表す。なお、S1からSKとなりうる候補は、例えば、多変量データの生成確率の場合には｛正規分布、対数正規分布、指数分布｝であったり、多項曲線出力の場合では、｛０次曲線、１次曲線、２次曲線、３次曲線｝であったりする。

なお、本明細書では具体的な例はすべて隠れマルコフを用いて説明をするが、その拡張モデル（例えば隠れセミマルコフモデル、因子化隠れマルコフモデルなど）の類似のモデルにも本発明を適用可能である。同様に、本明細書では、ターゲット変数をＸとした分布について説明しているが、回帰や判別のように、観測分布が条件付モデルＰ（Ｙ｜Ｘ）（Ｙはターゲットとなる確率変数）である場合に関しても本発明を適用可能である。

実施形態１．
図１は、本発明の隠れ変数モデル推定装置の第１の実施形態を示すブロック図である。隠れ変数モデル推定装置１００は、データ入力装置１０１と、隠れ状態数設定部１０２と、初期化処理部１０３と、隠れ変数変分確率計算処理部１０４と、モデル最適化処理部１０５と、最適性判定処理部１０６と、最適モデル選択処理部１０７と、モデル推定結果出力装置１０８とを備える。隠れ変数モデル推定装置１００には、入力データ１１１が入力され、入力データ１１１に対して隠れ状態数および観測確率の種類を最適化し、モデル推定結果１１２として出力する。

また、図２は、隠れ変数変分確率計算処理部１０４の例を示すブロック図である。隠れ変数変分確率計算処理部１０４は、前向き確率計算処理部１０４１と、正規化定数記憶部１０４２と、後ろ向き確率計算処理部１０４３と、前向き後ろ向き確率合算処理部１０４４を備える。隠れ変数変分確率計算処理部１０４には、入力データ１１１と、モデル最適化処理部１０５で推定された推定モデル１０４５とが入力され、隠れ変数変分確率１０４６と、前向き確率正規化定数１０４７とを出力する。

入力装置１０１は、入力データ１１１を入力するための入力インタフェース装置である。入力装置１０１には、入力データ１１１が入力される際に、観測確率の種類や、隠れ状態数の候補値など、モデルの推定に必要なパラメータも同時に入力される。

隠れ状態数設定部１０２は、入力された隠れ状態数の候補値からモデルの隠れ状態数を選択して設定する。以下、設定された隠れ状態数をＫと表記する。

初期化処理部１０３は、推定のための初期化処理を実施する。なお、初期化は任意の方法によって実施することが可能である。例としては、観測確率の種類を隠れ状態ごとにランダムに設定し、設定された種類に従って、各観測確率のパラメータをランダムに設定する方法や、隠れ変数の変分確率をランダムに設定する方法が挙げられる。

隠れ変数変分確率計算処理部１０４は、隠れ変数の変分確率を計算する。ここで、パラメータθは初期化処理部１０３あるいはモデル最適化処理部１０６で計算されているため、隠れ変数変分確率計算処理部１０４は、その値を利用する。隠れ変数変分確率計算処理部１０４は、次に定義する最適化基準Ａを最大化することによって変分確率を計算する。最適化基準Ａとは、周辺化対数尤度関数を完全変数に対する推定量（例えば最尤推定量や最大事後確率推定量）に関してラプラス近似した近似量の下界として定義される。なお、この下界は、完全変数に対する推定量の最適性と対数関数の凹性を用いることで導出することが可能である。

この手順を、隠れマルコフモデルを例に説明する。まず、周辺化対数尤度関数の下界を考える。この下界は、以下の式（２）で示される

なお、式（２）において、変分確率q（zN)を最大化する事で等号が成立する。ここで、完全変数に対する最尤推定量を用いて分子の完全変数の周辺化尤度をラプラス近似する事で、周辺化対数尤度関数の近似式として、以下に示す式（３）を得る。

ただし、上付きのバーは完全変数に対する最尤推定量を表す。また、D_＊は下付きのパラメータ＊の次元を表す。

次に、式（３）に対して最尤推定量が対数尤度関数を最大化する性質と、対数関数が凹関数である事を利用して、式（３）の下界を以下の式（４）のように算出する。

隠れ変数の変分分布q(zⁿtk)及びq(zⁿt-1k, zⁿtj)は、式（４）をqについて最大化する事によって算出される。ただし、上付き(i)を、隠れ変数変分確率計算処理部１０４、モデル最適化処理部１０５、最適判定処理部１０６の繰り返し計算における、(i)回目の繰り返しを表すとすると、q⁽ⁱ⁾はq~= q^(i-1)、θ＝θ^(i-1)と固定する。

なお、式（４）において下線を付した部分をＢとする。Ｂは、後述の式（８）で参照する。

図２を参照して、隠れ変数変分確率計算処理部１０４が備える要素について説明する。前向き確率計算処理部１０４１には、入力データ１１１と推定モデルが入力される。そして、前向き確率計算処理部１０４１は、時刻1から時刻tまでの観測（xⁿ1, ..., xⁿt）が得られた場合のzⁿtの確率を前向き確率として算出する。ただし、前向き確率は、最適化基準Aで算出されるモデル複雑度（例えば式（４）ではδtkに関する項）を考慮して算出される。また、前向き確率計算処理部１０４１は、zⁿtの確率の隠れ状態に関する和を１とするための正規化定数を、正規化定数記憶部１０４２に記憶させる。

同様に、後ろ向き確率計算処理部１０４３は、時刻t+1からTまでの観測（xⁿt+1, ..., xⁿT）が得られた場合のxⁿtの確率を後ろ向き確率として算出する。なお、後ろ向き確率の計算の際に、前向き確率算出と同時に得られる正規化定数を正規化定数記憶部１０４２から読み込む。ただし、後ろ向き確率は最適化基準Aで算出されるモデル複雑度（例えば式（４）ではδtkに関する項）を考慮して算出される。

最後に、前向き後ろ向き確率合算処理部１０４４は、前向き確率と後ろ向き確率から、変分分布を算出する。例えば、前向き後ろ向き確率合算処理部１０４４は、q(zⁿtk)をxⁿ1, ..., xⁿTが得られたときのzⁿtkの確率として計算する。前向き後ろ向き確率合算処理部１０４４は、前向き確率と後ろ向き確率の積として、以下の式（５）の計算によって、q(zⁿtk)を計算する。

また、前向き後ろ向き確率合算処理部１０４４は、xⁿ1, ..., xⁿt-1が得られたときのzⁿt-1jの確率と、隠れ状態jから隠れ状態kへ遷移する確率と、隠れ状態kにおいてxntが観測される確率と、（xⁿt+1, ..., xⁿT）が得られた場合のxⁿtの確率の積として、q(zⁿt-1j, zⁿtk)を算出する。具体的には、前向き後ろ向き確率合算処理部１０４４は、以下の式（６）の計算によってq(zⁿt-1j, zⁿtk)を算出する（式（６）左辺のq~の定義は式（７）を参照）。

この手順を隠れマルコフモデルを例に説明すると、前向き確率及び後ろ向き確率は、以下の式（７）の計算によって算出される。

ただし、f^tnk（式（７）の第１式）が前向き確率を表し、b^tnk（式（７）の第２式）が後ろ向き確率を表す。より具体的には、式（７）において、前向き確率と後ろ向き確率の両者とも、漸化式として記述されている。そして、前向き確率はt=1から順に算出することが可能であり、後ろ向き確率はt=Tから順に算出することが可能である。なお、正規化定数はζ^tnで算出される。後ろ向き確率計算処理部１０４３は、前向き確率計算処理部１０４１が前向き確率を算出する際に計算した正規化定数を利用して後ろ向き確率を算出すればよい。

また、式（５）の第３式において、δに関する乗算が含まれているが、これは、最適化基準Aで算出されるモデル複雑度を考慮していることを意味する。

モデル最適化処理部１０５は、式（４）に対してモデル（パラメータθ及びその種類S） を最適化する。具体的には、q及びq~を隠れ変数変分確率計算処理部１０４で算出された隠れ変数の変分分布（q⁽ⁱ⁾）に固定し、式（４）におけるGを最大化するモデルを算出する。この処理において重要な点は、式（４）によって定義されたGは、コンポーネントごとに最適化の関数を分解する事が可能なため、コンポーネント種類の組合せ(S1からSKのどの種類を指定するか)を考慮することなく、S1からSK及びパラメータφ₁からφ_Kを別々に最適化する事が可能な点である。これによって、コンポーネントの種類を最適化する際に、組み合わせ爆発を回避して最適化を実行する事が可能となる。

最適性判定処理部１０６は、式（４）で計算される最適化基準Aの収束を判定する。最適化基準Aが収束したと判定していない場合には、隠れ変数変分確率計算処理部１０４から最適性判定処理部１０６の処理を繰り返す。なお、式（４）で計算される最適化基準Aの算出は、隠れ状態が独立ではないためΣzn q(zn) log q(zn)の計算に指数時間の計算量を必要とするが、正規化定数記憶部１０４２に記憶されている正規化定数を利用して効率的に計算をする事が可能である。例えば隠れマルコフモデルの場合には、以下の式（８）のように計算される。

式（８）に示すＢは、式（４）において下線を付した部分である。

隠れ変数変分確率計算処理部１０４から最適性判定処理部１０６の処理を繰り返し、変分分布とモデルを更新する事で、適切なモデルを選択する事が可能となる。なお、この繰り返しによって最適化基準Aが単調に増加する事が保証される。

最適化基準Aが収束している場合、隠れ状態数設定部１０２で設定された隠れ状態数Kに対して、隠れ変数変分確率計算処理部１０４から最適性判定処理部１０６のループ処理で算出される最適化基準Aと、その１つ前のループ処理で算出された最適化基準Aのうち、大きい方の最適化基準Aに対応するモデルを最適なモデルとして設定する。全ての候補値についてモデルの最適化が完了した場合には、処理がモデル推定結果出力装置１０８へ移り、まだ最適化の済んでいない候補が存在する場合には、隠れ状態数設定部１０２へ処理が移る。

モデル推定結果出力装置１０８は、最適な隠れ状態数、観測確率の種類、パラメータ、変分分布などをモデル推定結果出力結果１１２として出力する。

隠れ状態数設定部１０２、初期化処理部１０３、隠れ変数変分確率計算処理部１０４（前向き確率計算処理部１０４１、正規化定数記憶部１０４２、後ろ向き確率計算処理部１０４３、前向き後ろ向き確率合算処理部１０４４）、モデル最適化処理部１０５、最適性判定処理部１０６、最適モデル選択処理部１０７およびモデル推定結果出力装置１０８は、例えば、隠れ変数モデル推定プログラムに従って動作するコンピュータのＣＰＵによって実現される。ＣＰＵが、隠れ変数モデル推定プログラムを記録したコンピュータ読み取り可能な記録媒体から隠れ変数モデル推定プログラムを読み込み、上記の各要素として動作すればよい。

また、隠れ状態数設定部１０２、初期化処理部１０３、隠れ変数変分確率計算処理部１０４、モデル最適化処理部１０５、最適性判定処理部１０６、最適モデル選択処理部１０７およびモデル推定結果出力装置１０８が別々のハードウェアで実現されていてもよい。また、隠れ変数変分確率計算処理部１０４においても、前向き確率計算処理部１０４１、正規化定数記憶部１０４２、後ろ向き確率計算処理部１０４３、前向き後ろ向き確率合算処理部１０４４が別々のハードウェアで実現されていてもよい。

図３は、本発明の第１の実施形態の処理経過の例を示すフローチャートである。データ入力装置１０１を介して入力データ１１１が入力される（ステップＳ１００）。

次に、隠れ状態数設定部１０２は、入力された隠れ状態数の候補値のうち、まだ最適化の行なわれていない隠れ状態数を選択し設定する（ステップＳ１０１）。

次に、初期化処理部１０３は、設定された隠れ状態数に対して、推定のため、パラメータや隠れ変数変分確率の初期化処理を実施する（ステップＳ１０２）。

次に、隠れ変数変分確率計算処理部１０４は、隠れ変数の変分確率を計算する（ステップＳ１０３）。

次に、モデル最適化処理部１０５は、各隠れ状態に対して観測確率の種類とパラメータの推定を実施する（ステップＳ１０４）。この処理は、各隠れ状態のモデルの最適化であると言うことができる。

次に、最適性判定処理部１０６は、最適化基準A が収束したかを判定する。（例えば、Ｓ１０５）。最適性判定処理部１０６は、直近のステップＳ１０３〜Ｓ１０５のループ処理において得られた最適化基準Aと、その１つ前のステップＳ１０３〜Ｓ１０５のループ処理で得られた最適化基準Aの差を計算し、その差の絶対値が予め定められた閾値以下になっていれば、最適化基準A が収束したと判定してよい。また、その差の絶対値が閾値より大きければ、最適性判定処理部１０６は、最適化基準A が収束していないと判定してよい。なお、絶対値による最適化基準Aの差の算出は一例であり、例えば相対的な差によって収束を判定する等の方法を採用してもよい。

ステップＳ１０５において、最適化基準A が収束していないと判定された場合には、ステップＳ１０３〜Ｓ１０５の処理を繰り返す。

ステップＳ１０５において、最適化基準A が収束したと判定された場合、最適モデル選択処理部１０７は、直近のステップＳ１０３〜Ｓ１０５のループ処理において最適化されたモデル（隠れ状態数、観測確率の種類、パラメータ）の最適化基準Aと、その１つ前のステップＳ１０３〜Ｓ１０５のループ処理で最適化されたモデルの最適化基準Aとを比較し、値の大きい方の最適化基準Aに対応するモデルを、最適なモデルとして設定する（ステップＳ１０６）。

次に、隠れ状態数設定部１０２は、推定されていない隠れ状態数の候補が残っているか否かを判定する（ステップＳ１０７）。隠れ状態数の候補が残っている場合には、ステップＳ１０２〜Ｓ１０７の処理を繰り返す。一方、隠れ状態数の候補が残っていない場合には、モデル推定結果出力装置１０８がモデル推定結果を出力し（ステップＳ１０８）、処理を完了する。

図４は、本実施形態における隠れ変数変分確率計算処理部１０４の動作（換言すれば、ステップＳ１０３の処理経過）を示すフローチャートである。

前向き確率計算処理部１０４１は、n番目のデータのt番目の時刻に対する前向き確率f^t(i)nkを算出する（ステップＳ１１１）。このとき、前向き確率計算処理部１０４１は、正規化定数も算出し、正規化定数記憶部１０４２に記憶させる（ステップＳ１１２）。

続いて、前向き確率計算処理部１０４１は、すべての時刻tに対して前向き確率の算出が完了したかを確認し（ステップＳ１１３）、未完了の場合にはステップＳ１１１，Ｓ１１２の処理を繰り返す。完了した場合には、ステップＳ１１４の処理に移る。

後ろ向き確率計算処理部１０４３は、n番目のデータのt番目の時刻に対する後ろ向き確率b^t(i)nkを算出する（ステップＳ１１４）。そして、後ろ向き確率計算処理部１０４３は、すべての時刻tに対して後ろ向き確率の算出が完了したかを確認し（ステップＳ１１５）、未完了の場合にはステップＳ１１４の処理を繰り返す。完了した場合には、ステップＳ１１６の処理に移る。

前向き後ろ向き確率合算処理部１０４４は、n番目のデータのすべての時刻に対して前向き確率、後ろ向き確率の合算処理を行い、変分分布の計算を行う（ステップＳ１１６）。

続いて、前向き後ろ向き確率合算処理部１０４４は、nに関してすべてのデータに対して変分分布の算出処理が完了しているかを確認する（ステップＳ１１７）。未完了の場合には、ステップＳ１１１以降の処理を繰り返し、完了した場合には、処理を終了する。

隠れ状態数および観測確率の種類と共にモデルの候補数が指数的に増加する場合であっても、上記のような本発明の動作（特に、隠れ変数変分確率計算処理部１０４の動作）によって、高速にモデル選択を実現することができる。

また、非特許文献４に記載の技術と本発明を比較すると、既に説明したように、非特許文献４に記載の技術では、隠れ変数間の独立性が必要となり、系列依存性のある隠れ変数モデルに適用できなかった。これに対し本発明は、系列依存性を持つ多変量データの隠れ変数モデルを推定することができる。

実施形態２．
図５は、本発明の隠れ変数モデル推定装置の第２の実施形態を示すブロック図である。第２の実施形態の隠れ変数モデル推定装置２００は、第１の実施形態の隠れ変数モデル推定装置１００（図１参照）と比較して、最適モデル選択処理部１０７を備えず、隠れ状態数選択処理部２０１を備える。

データ入力装置１０１、隠れ状態数設定部１０２、初期化処理部１０３、隠れ変数変分確率計算処理部１０４、モデル最適化処理部１０５、最適性判定処理部１０６およびモデル推定結果出力装置１０８に関しては、第１の実施形態と同様である。

また、第１の実施形態の隠れ変数モデル推定装置１００は、隠れ状態数の候補に対してモデル最適化を行い、最適化基準Aを最大化するモデルを選択する。これに対して、第２の実施形態の隠れ変数モデル推定装置２００では、隠れ変数変分確率計算処理部１０４の処理の後で、隠れ状態数選択処理部２０１が、小さくなった隠れ状態をモデルから除去する。

具体的には、隠れ状態数選択処理部２０１は、隠れ変数変分確率計算処理部１０４で算出されたq(zⁿtk)に対して、以下の式（９）の状態を満たす隠れ状態を除去する。

式（９）の右辺に示すεは、入力データ１１１と同時に入力される閾値である。すなわち、隠れ状態数選択処理部２０１は、閾値ε以下である隠れ状態を除去する。

式（９）による隠れ状態の除去は、以下の理由によって正しく動く事が説明される。まず、式（７）の前向き確率を観察すると、小さい隠れ状態（すなわち小さいδ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を基礎とする優先権を主張し、その開示の全てをここに取り込む。

以上、実施形態を参照して本願発明を説明したが、本願発明は上記の実施形態に限定されるものではない。本願発明の構成や詳細には、本願発明のスコープ内で当業者が理解し得る様々な変更をすることができる。

産業上の利用の可能性

本発明は、系列依存性を持つ多変量データの隠れ変数モデル推定装置に好適に適用される。

１０１ データ入力装置
１０２ 隠れ状態数設定部
１０３ 初期化処理部
１０４ 隠れ変数変分確率計算処理部
１０５ モデル最適化処理部
１０６ 最適性判定処理部
１０７ 最適モデル選択処理部
１０８ モデル推定結果出力装置
２０１ 隠れ状態数選択処理部

Claims (12)

1. 周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値を最大化することによって変分確率を計算する変分確率計算部と、
   各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定するモデル推定部と、
   変分確率計算部が変分確率を計算する際に用いた基準値が収束したか否かを判定する収束判定部とを備える
   ことを特徴とする隠れ変数モデル推定装置。
2. 変分確率計算部が変分確率を計算し、モデル推定部が最適な隠れ変数モデルを推定し、収束判定部が、基準値が収束したか否かを判定するループ処理を繰り返し、
   基準値が収束したときに、当該基準値と、その１つ前のループ処理における基準値とを比較し、大きい方の基準値に対応する隠れ変数モデルを最適な隠れ変数モデルとして選択する最適モデル選択部を備える
   請求項１に記載の隠れ変数モデル推定装置。
3. 変分確率計算部が変分確率を計算し、モデル推定部が最適な隠れ変数モデルを推定し、収束判定部が、基準値が収束したか否かを判定するループ処理を繰り返し、
   変分確率計算部の計算結果に応じて、所定の条件を満たす隠れ状態を除去する隠れ状態除去部を備える
   請求項１に記載の隠れ変数モデル推定装置。
4. モデル推定部は、隠れ変数モデルとして隠れマルコフモデルを推定する請求項１から請求項３のうちのいずれか１項に記載の隠れ変数モデル推定装置。
5. 周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値を最大化することによって変分確率を計算し、
   各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定し、
   変分確率を計算する際に用いた基準値が収束したか否かを判定する
   ことを特徴とする隠れ変数モデル推定方法。
6. 変分確率を計算し、最適な隠れ変数モデルを推定し、基準値が収束したか否かを判定するループ処理を繰り返し、
   基準値が収束したときに、当該基準値と、その１つ前のループ処理における基準値とを比較し、大きい方の基準値に対応する隠れ変数モデルを最適な隠れ変数モデルとして選択する
   請求項５に記載の隠れ変数モデル推定方法。
7. 変分確率を計算し、最適な隠れ変数モデルを推定し、基準値が収束したか否かを判定するループ処理を繰り返し、
   変分確率の計算結果に応じて、所定の条件を満たす隠れ状態を除去する
   請求項５に記載の隠れ変数モデル推定方法。
8. 隠れ変数モデルとして隠れマルコフモデルを推定する請求項５から請求項７のうちのいずれか１項に記載の隠れ変数モデル推定方法。
9. コンピュータに、
   周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値を最大化することによって変分確率を計算する変分確率計算処理、
   各隠れ状態に対して観測確率の種類とパラメータを推定することで最適な隠れ変数モデルを推定するモデル推定処理、および、
   変分確率計算処理で変分確率を計算する際に用いた基準値が収束したか否かを判定する収束判定処理
   を実行させるための隠れ変数モデル推定プログラム。
10. コンピュータに、
    変分確率計算処理、モデル推定処理、収束判定処理のループ処理を繰り返し実行させ、
    基準値が収束したときに、当該基準値と、その１つ前のループ処理における基準値とを比較し、大きい方の基準値に対応する隠れ変数モデルを最適な隠れ変数モデルとして選択する最適モデル選択処理
    を実行させる請求項９に記載の隠れ変数モデル推定プログラム。
11. コンピュータに、
    変分確率計算処理、モデル推定処理、収束判定処理のループ処理を繰り返し実行させ、
    変分確率計算処理の計算結果に応じて、所定の条件を満たす隠れ状態を除去する隠れ状態除去処理
    を実行させる請求項９に記載の隠れ変数モデル推定プログラム。
12. コンピュータに、
    モデル推定処理で、隠れ変数モデルとして隠れマルコフモデルを推定させる
    請求項９から請求項１１のうちのいずれか１項に記載の隠れ変数モデル推定プログラム。

Images