(57)【特許請求の範囲】
【請求項１】
因子化隠れマルコフモデルの各layer1隠れ変数によって決まるパラメータの線形結合として表現された観測モデルのパラメータに関するヘッセ行列の行列式の近似値を計算する計算部と、
前記行列式の近似値を利用して隠れ変数の変分確率を計算する変分確率計算部と、
変分分布に基づいて隠れ状態を削除する隠れ状態削除部と、
周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値に対してパラメータを最適化するとともに、基準値を計算するパラメータ最適化部と、
前記基準値が収束したか否かを判定する収束判定部とを
備えることを特徴とするモデル推定装置。
【請求項２】
計算部がヘッセ行列の行列式の近似値を計算し、変分確率計算部が隠れ変数の変分確率を計算し、隠れ状態削除部が隠れ状態を削除し、パラメータ最適化部がパラメータを最適化し、計算部がヘッセ行列の行列式の近似値を計算し、パラメータ最適化部が基準値を計算し、収束判定部が、前記基準値が収束したか否かを判定するループ処理を、収束判定部が、前記基準値が収束したと判定するまで繰り返す
請求項１に記載のモデル推定装置。
【請求項３】
モデル推定装置は、因子化隠れマルコフモデル推定装置である請求項１または請求項２に記載のモデル推定装置。
【請求項４】
コンピュータが、
因子化隠れマルコフモデルの各layer1隠れ変数によって決まるパラメータの線形結合として表現された観測モデルのパラメータに関するヘッセ行列の行列式の近似値を計算し、
前記行列式の近似値を利用して隠れ変数の変分確率を計算し、
変分分布に基づいて隠れ状態を削除し、
周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値に対してパラメータを最適化し、
ヘッセ行列の行列式の近似値を計算し、
前記基準値を計算し、
前記基準値が収束したか否かを判定する
ことを特徴とするモデル推定方法。
【請求項５】
コンピュータが、
ヘッセ行列の行列式の近似値を計算し、隠れ変数の変分確率を計算し、隠れ状態を削除し、パラメータを最適化し、ヘッセ行列の行列式の近似値を計算し、前記基準値を計算し、前記基準値が収束したか否かを判定するループ処理を、前記基準値が収束するまで繰り返す
請求項４に記載のモデル推定方法。
【請求項６】
モデル推定方法は、因子化隠れマルコフモデル推定方法である請求項４または請求項５に記載のモデル推定方法。
【請求項７】
コンピュータに、
因子化隠れマルコフモデルの各layer1隠れ変数によって決まるパラメータの線形結合として表現された観測モデルのパラメータに関するヘッセ行列の行列式の近似値を計算する計算処理、
前記行列式の近似値を利用して隠れ変数の変分確率を計算する変分確率計算処理、
変分分布に基づいて隠れ状態を削除する隠れ状態削除処理、
周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値に対してパラメータを最適化するパラメータ最適化処理、
前記基準値を計算する基準値計算処理、および、
前記基準値が収束したか否かを判定する収束判定処理
を実行させるためのモデル推定プログラム。
【請求項８】
コンピュータに、
計算処理、変分確率計算処理、隠れ状態削除処理、パラメータ最適化処理、計算処理、基準値計算処理、および収束判定処理のループ処理を、基準値が収束したと判定されるまで繰り返し実行させる
請求項７に記載のモデル推定プログラム。
【請求項９】
モデル推定プログラムは、因子化隠れマルコフモデル推定プログラムである請求項７または請求項８に記載のモデル推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、モデル推定装置、モデル推定方法およびモデル推定プログラムに関し、特に、モデル事後確率を近似し、その下界を最大化する事によって因子化隠れマルコフモデルを推定するモデル推定装置、モデル推定方法およびモデル推定プログラムに関する。
【背景技術】
【０００２】
自動車から取得されるセンサデータ、健康診断の検査値履歴、電力需要履歴などに代表されるデータは、すべて「時間的な依存性」を持つ多変量データである。このようなデータの分析は、産業上重要な数多くの分野に適用される。例えば、自動車から取得されるセンサデータの分析によって、自動車の故障原因を解析して素早い修理を実現することが考えられる。また、健康診断の検査値履歴の分析によって、疾患のリスクの推定および疾患の予防を実現することが考えられる。また、電力需要履歴の分析によって、電力の需要を予測して過不足に備えられるようにすることが考えられる。
【０００３】
このようなデータに対しては、時間依存性を持つ隠れ変数モデル（例えば、隠れマルコフモデル）を用いてモデル化する事が一般的である。例えば、隠れマルコフモデルを利用するためには、隠れ状態数、観測確率分布の種類、及び分布パラメータを決定する必要がある。隠れ状態数と観測確率分布の種類がわかっている場合には、Expectation Maximization法（例えば、非特許文献１参照）を利用してパラメータを推定する事が可能である。
【０００４】
隠れマルコフモデルでは、１つの隠れ状態に対して観測が独立に得られる（換言すれば、隠れ状態は時間的に依存している）という仮定をおく。しかし、実応用のさまざまな場面において、この仮定が成立していないケースが存在する。すなわち、複数の隠れ状態の重ね合わせによって観測が得られるケースが存在する。例えば、音声認識において、複数の人が話している状況の音声データに対しては、各人に対する隠れ状態が存在する。そして、観測される音声は、各隠れ状態からの観測の重ね合わせである。すなわち、観測される音声には、複数の人の声が混じっている。また、例えば、複数の赤外線センサが設置された屋内で観測されるセンサ反応から、人間の行動を推定する例を考える。このセンサ反応は、人間の位置情報を表す。複数の人が屋内で行動している場合には、人間の行動（隠れ状態）は時間的に遷移している。そして、赤外線センサの反応（観測）は複数の人の行動の重ねあわせによって決まる。
【０００５】
このような複雑な系列データを扱うために、因子化隠れマルコフモデルが提案されている（例えば、非特許文献２参照）。因子化隠れマルコフモデルでは、複数の隠れ状態の時間遷移が考慮され、観測モデルのパラメータが各隠れ状態に依存して決まる。
【０００６】
因子化隠れマルコフモデルでは、独立に時間遷移する上位（layer1）の隠れ状態と、各layer1の隠れ変数に対する状態を表す下位（layer2）の隠れ状態とがある。
【０００７】
隠れ状態数を決定する問題は、一般的に「モデル選択問題」や「システム同定問題」と呼ばれ、信頼性のあるモデルを構築するために極めて重要な問題である。そのための技術が種々提案されている。
【０００８】
例えば、非特許文献２では、隠れ状態数を決定する方法として、変分ベイズ法によって、変分自由エネルギーを最大化する方法が提案されている。以下、この方法を第１の公知技術と記す。
【０００９】
また、例えば、非特許文献３では、隠れ状態数を決定する方法として、階層Dirichlet過程事前分布を用いたノンパラメトリックベイズ法が提案されている。以下、この方法を第２の公知技術と記す。
【００１０】
隠れマルコフモデルでは、隠れ変数が時間依存性を持ち、隠れ変数に対してパラメータが独立である。そして、非特許文献４では、隠れマルコフモデルに適用される技術として、因子化漸近ベイズ推論（Factorized Asymptotic Bayesian Inference）という技術が提案されている。この技術は、変分ベイズ法やノンパラメトリックベイズ法に対して、速度・精度の面で優位性を持つ。
【００１１】
また、混合分布モデルでは、隠れ変数が時間依存性を持たない。そして、非特許文献５では、混合分布モデルに対する因子化漸近ベイズ推論が提案されている。
【００１２】
また、完全周辺尤度関数を近似し、その下界を最大化することが、非特許文献４、非特許文献５に記載されている。
【先行技術文献】
【非特許文献】
【００１３】
【非特許文献１】C. Bishop, "Pattern Recognition and Machine Learning", Springer, 2007.
【非特許文献２】Ghahramani, Zoubin; Jordan, Michael I. (1997). "Factorial Hidden Markov Models".  Machine Learning  29  (2/3): 245-273
【非特許文献３】J. V. Gael, Y. W. Teh, and Z. Ghahramani. "The infinite factorial hidden Markov model". In Neural Information Processing Systems (NIPS), 2008.
【非特許文献４】Ryohei Fujimaki, Kohei Hayashi: "Factorized Asymptotic Bayesian Hidden Markov Model". Proceedings of the 25th international conference on machine learning (ICML), 2012
【非特許文献５】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
【発明の概要】
【発明が解決しようとする課題】
【００１４】
第１の公知技術は、周辺化尤度関数の下界を最大化する際に、変分分布上における隠れ状態と分布パラメータの独立性を仮定する。そのため、第１の公知技術では、周辺化尤度の近似精度が悪くなるという問題がある。
【００１５】
第２の公知技術には、モデルが複雑になるため計算量が非常に多くなるという問題や、layer1とlayer2の隠れ状態数を同時には推定できないという問題がある。また、入力パラメータによって結果が大きく変わるという問題もある。
【００１６】
また、非特許文献４、非特許文献５等に記載された技術では、本質的に隠れ変数に関するパラメータの独立性が重要である。そのため、因子化隠れマルコフモデルのように隠れ変数に対してパラメータが依存関係をもつモデルに直接、因子化漸近ベイズ推論を適用することができなかった。
【００１７】
本発明は、因子化漸近ベイズ推論に基づいて、因子化隠れマルコフモデルに関するモデル選択問題を解決することができる因子化隠れマルコフモデル推定装置、因子化隠れマルコフモデル推定方法および因子化隠れマルコフモデル推定プログラムを提供することを目的とする。
【課題を解決するための手段】
【００１８】
本発明によるモデル推定装置は、因子化隠れマルコフモデルの各layer1隠れ変数によって決まるパラメータの線形結合として表現された観測モデルのパラメータに関するヘッセ行列の行列式の近似値を計算する計算部と、その行列式の近似値を利用して隠れ変数の変分確率を計算する変分確率計算部と、変分分布に基づいて隠れ状態を削除する隠れ状態削除部と、周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値に対してパラメータを最適化するとともに、基準値を計算するパラメータ最適化部と、基準値が収束したか否かを判定する収束判定部とを備えることを特徴とする。
【００１９】
また、本発明によるモデル推定方法は、コンピュータが、因子化隠れマルコフモデルの各layer1隠れ変数によって決まるパラメータの線形結合として表現された観測モデルのパラメータに関するヘッセ行列の行列式の近似値を計算し、行列式の近似値を利用して隠れ変数の変分確率を計算し、変分分布に基づいて隠れ状態を削除し、周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値に対してパラメータを最適化し、ヘッセ行列の行列式の近似値を計算し、基準値を計算し、基準値が収束したか否かを判定することを特徴とする。
【００２０】
また、本発明によるモデル推定プログラムを記録したコンピュータ読み取り可能な記録媒体は、コンピュータに、因子化隠れマルコフモデルの各layer1隠れ変数によって決まるパラメータの線形結合として表現された観測モデルのパラメータに関するヘッセ行列の行列式の近似値を計算する計算処理、その行列式の近似値を利用して隠れ変数の変分確率を計算する変分確率計算処理、変分分布に基づいて隠れ状態を削除する隠れ状態削除処理、周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値に対してパラメータを最適化するパラメータ最適化処理、基準値を計算する基準値計算処理、および、基準値が収束したか否かを判定する収束判定処理を実行させるためのモデル推定プログラムを記憶する。
【発明の効果】
【００２１】
本発明では、因子化漸近ベイズ推論に基づいて、因子化隠れマルコフモデルに関するモデル選択問題を解決することができる。
【図面の簡単な説明】
【００２２】
【図１】図１は、本発明の因子化隠れマルコフモデル推定装置の構成例を示すブロック図である。
【図２】図２は、本発明の処理経過の例を示すフローチャートである。
【図３】図３は、本発明の概要を示すブロック図である。
【発明を実施するための形態】
【００２３】
本発明の貢献を明確化するために、まず隠れマルコフモデルと因子化隠れマルコフモデルについて違いと、因子化漸近ベイズ推論がなぜ直接、因子化隠れマルコフモデルに適用できないかの問題点を説明する。
【００２４】
以下の説明では時間依存するデータ列xn(n=1,...,N)が入力されると仮定する。ここで、各xnは、長さTnの多変量データ列（xn=xn1, ..., xnTn，t=1,...,N)であるとする。また、各xntはD次元の観測ベクトルであり、xnt=(xnt1, ..., xntD)である。次に、観測変数xntに対するlayer1の隠れ変数znt=(znt1, ..., zntK)を定義する。ここで、Kはlayer1の隠れ状態数である。
【００２５】
隠れマルコフモデルでは、layer1の隠れ変数に対して、znt1がバイナリ変数であり、
【数１】
0006398990.tif 000002
である。すなわち、zntは、１つの要素だけが1となる。隠れマルコフモデルでは、xnとznの同時分布を、
【数２】
0006398990.tif 000003
と表す。ただし、j=(r，s，h)は、それぞれ隠れ状態の初期確率、隠れ状態の遷移確率、観測確率のパラメータである。また、観測確率は、以下の式１に示すように分解される。
【００２６】
【数３】
0006398990.tif 000004
【００２７】
ただし、h=(h1, ..., hk)である。ここで重要な点は、hkはkに関して依存関係を持たない点である。この性質によって、後述するように、対数同時尤度のヘッセ行列がブロック対角となり、因子化漸近ベイズ推論の理論的に優れた性質を利用することが可能となる。
【００２８】
次に、因子化隠れマルコフモデルを説明する。因子化隠れマルコフモデルでは、k番目のlayer1隠れ変数は、layer2の隠れ変数のベクトルzntk=(zntk1, ..., zntkMk)として表される。Mkは、k番目のlayer1隠れ状態数である。zntkmは、バイナリ変数であり、
【数４】
0006398990.tif 000005
である。すなわち、zntkは、１つの要素だけが1となる。因子化隠れマルコフモデルが、隠れマルコフモデルと最も異なる点は、観測確率のパラメータがlayer2の隠れ変数に依存する点である。これを、正規分布を例に説明する。なお、指数分布族など、より広いクラスの確率分布においても、以下に示す説明と同様の議論が成立する。観測モデルのパラメータを、各layer1隠れ変数によって決まるパラメータの線形結合として表現すると、以下の式２に示すようになる。
【００２９】
【数５】
0006398990.tif 000006
【００３０】
ただし、k番目のlayer1隠れ変数のm番目のlayer２隠れ変数に対応するパラメータを、hkmと表記する。また、N(x, a, b)は、xに関する平均a, 共分散bの正規分布であるとする。ここで重要な点は、各サンプルのd次元目の観測分布が隠れ状態に依存している点である。換言すれば、重要な点は、正規分布の平均が隠れ変数によって決まる点である。このモデルでは、隠れマルコフモデルとは異なり、
【数６】
0006398990.tif 000007
のヘッセ行列がブロック対角とならず、因子化漸近ベイズ推論の理論的に優れた性質（不要な隠れ状態の削除など）が失われてしまう。
【００３１】
上記の説明を、より具体的に説明する。なお、表記の簡単化のために、添え字nについては省略するが、本発明と非特許文献４で本質的に異なる点の説明には影響しない。まず、本発明では、ベイズ推論にしたがって、周辺化対数尤度を最大化することで、モデル及びパラメータを最適化する。ただし、周辺化対数尤度を直接最適化することは難しいため、まず周辺化対数尤度を式３のように変形する。
【００３２】
【数７】
0006398990.tif 000008
【００３３】
ただし、Mはモデルであり、q(z)はzに関する変分分布である。また、max_qは、qに関する最大値であることを表す。ここで、同時周辺化尤度p(x, z|M)は、パラメータに対する積分形式で、以下に示す式４のように変形可能である。
【００３４】
【数８】
0006398990.tif 000009
【００３５】
以下では、隠れマルコフモデルと、因子化隠れマルコフモデルが本質的に異なる観測モデルに限定して、説明を行う。なお、Nは、隠れマルコフモデルと因子化隠れマルコフモデルの違いに影響を与えないため、以下では、Nおよびnに関するノーテーションは全て省略する。
【００３６】
まず、隠れマルコフモデルでは、式１より
【数９】
0006398990.tif 000010
となる。そして、最尤推定量h’の周りで各項をテイラー展開し、３次以上の項を無視することで、log p(x| z，h)は、式５に示すように近似される。
【００３７】
【数１０】
0006398990.tif 000011
【００３８】
上記の式５は、非特許文献４に記載された（３）式の観測に関する項に対応する。ここで、Fkは、p(xt|hk)のヘッセ行列を
【数１１】
0006398990.tif 000012
で割った行列であり、p(x, z|j)のヘッセ行列のブロック対角項である。式５の近似に関して式４の積分を計算すると、関連する項は、以下に示す式６となる。なお、両辺にexpをかけて、左辺のlogを除いたものが、非特許文献４に記載された（５）式の観測に関する項に対応する。
【００３９】
【数１２】
0006398990.tif 000013
【００４０】
ただし、detは引数の行列式を表す。また、Dkは、hkの次元を表すものとする。ここで、Tの極限を考えると、
【数１３】
0006398990.tif 000014
log det(Fk)は相対的に小さくなるため、無視することが可能であり、最適化基準として非特許文献４に記載された（６）式及び（７）式が得られる。
【００４１】
なお、Tは、データの長さ（シーケンスデータのシーケンス長）である。
【００４２】
次に、因子化隠れ変数モデルでは、式２より、
【数１４】
0006398990.tif 000015
となる。そして、最尤推定量h’の周りで各項をテイラー展開し、３次以上の項を無視することで、log p(xd| z，h)は、式７に示すように近似される。
【００４３】
【数１５】
0006398990.tif 000016
【００４４】
ただし、Fdは、
【数１６】
0006398990.tif 000017
のヘッセ行列をTで割った行列である。式７の近似に関して式４の積分を計算すると、関連する項は、以下に示す式８となる。
【００４５】
【数１７】
0006398990.tif 000018
【００４６】
隠れマルコフモデルのケース（式６）と因子化隠れマルコフモデルのケース（式８）で本質的に異なるのは、式６の
【数１８】
0006398990.tif 000019
という項が、式８ではDd log Nになり、モデル複雑性が隠れ変数に依存しなくなっている点である。より具体的に説明する。非特許文献４で提案されている隠れマルコフモデルに対する因子化漸近ベイズ推論では、モデルの複雑性が隠れ変数に依存していることによって、不要な隠れ状態の削除など、理論的に優れた性質を持つ。この点は、例えば、非特許文献４の"Section 4.4 Automatic Hidden State Selection"で説明されている。一方、非特許文献４の手順を因子化隠れマルコフモデルに適用して得られる式８では、そのような性質が失われてしまう。
【００４７】
本発明が、非特許文献４に記載された技術と本質的に異なる点は、log det(Fd))  の近似処理にある。非特許文献４に記載の手順では、式８のlog det(Fd)を、漸近的に小さいとして、以下に示すように近似する。
【００４８】
【数１９】
0006398990.tif 000020
【００４９】
これに対し、本発明では、log det(Fd)を、以下の式９に示すように近似する。具体的には、本発明の因子化隠れマルコフモデル推定装置に設けられる後述の情報量基準近似処理部１０５が、式９によりlog det(Fd)を近似する。
【００５０】
【数２０】
0006398990.tif 000021
【００５１】
この処理で重要な点は、非特許文献４に記載の処理では無視されているlog det(Fd)が、隠れ変数に依存する点である。なお、（式９）の近似は、Tが無限大の極限に等しくなる。以下の式のような非特許文献４における近似処理では、モデルの複雑性を過小評価している。
【００５２】
【数２１】
0006398990.tif 000022
【００５３】
本発明では、式９に示す近似を行うことで、不要な隠れ状態の削除など、理論的に優れた性質を持ったモデル選択を、因子化隠れ変数モデルに対して実現することができる。
【００５４】
以下、本発明の実施形態を図面を参照して説明する。
【００５５】
図１は、本発明の因子化隠れマルコフモデル推定装置の構成例を示すブロック図である。本発明の因子化隠れマルコフモデル推定装置１００は、データ入力装置１０１と、隠れ状態数設定部１０２と、初期化処理部１０３と、隠れ変数変分確率計算処理部１０４と、情報量基準近似処理部１０５と、隠れ状態選択処理部１０６と、パラメータ最適化処理部１０７と、最適性判定処理部１０８と、モデル推定結果出力装置１０９とを備えている。因子化隠れマルコフモデル推定装置１００には、入力データ１１１が入力される。因子化隠れマルコフモデル推定装置１００は、入力データ１１１に対して因子化隠れマルコフモデルを最適化し、モデル推定結果１１２として出力する。
【００５６】
データ入力装置１０１は、入力データ１１１を入力するための装置である。観測確率の種類や、隠れ状態数の候補値など、モデルの推定に必要なパラメータが入力データ１１１として同時にデータ入力装置１０１に入力される。
【００５７】
隠れ状態数設定部１０２は、モデルのlayer1隠れ状態数Kを、入力データ１１１として入力された最大値Kmaxに設定する。また、隠れ状態数設定部１０２は、layer2隠れ状態数Mkを、入力データ１１１として入力された最大値Mmaxに設定する。すなわち、隠れ状態数設定部１０２によって、K=Kmax, Mk=Mmaxと設定される。
【００５８】
初期化処理部１０３は、推定のための初期化処理を実施する。なお、初期化は任意の方法によって実施することが可能である。例としては、各観測確率のパラメータjをランダムに設定する方法や、隠れ変数の変分確率をランダムに設定する方法が挙げられる。
【００５９】
隠れ変数変分確率計算処理部１０４は、隠れ変数の変分確率を計算する。ここで、パラメータjは初期化処理部１０３あるいはパラメータ最適化処理部１０７で算出されているため、隠れ変数変分確率計算処理部１０４は、その値を利用する。隠れ変数変分確率計算処理部１０４は、以下のように定義する最適化基準Aを最大化することによって変分確率を計算する。最適化基準Aとは、周辺化対数尤度関数を、完全変数に対する推定量（例えば最尤推定量や最大事後確率推定量）に関してラプラス近似した近似量の下界として定義される。
【００６０】
情報量基準近似処理部１０５は、隠れ変数変分確率計算処理部１０４やパラメータ最適か処理部１０７で必要な、ヘッセ行列の行列式の近似処理を行う。具体的には、情報量基準近似処理部１０５は、前述の式９による近似計算を行う。
【００６１】
隠れ状態選択処理部１０６は、隠れ状態のうち、小さな状態をモデルから除去する。具体的には、隠れ状態選択処理部１０６は、k番目の隠れ状態について、
【数２２】
0006398990.tif 000023
が、入力データ１１１として設定された閾値を下回った場合には、モデルから、layer2のm番目の隠れ状態を削除する。また、隠れ状態選択処理部１０６は、layer2の隠れ状態数が1になった場合には、layer1の隠れ状態を削除する。
【００６２】
パラメータ最適化処理部１０７は、隠れ変数の変分確率を固定した上で、最適化基準Aに対してjを最適化する。なお、最適化基準Aのjに関する項は、隠れ状態の変分分布によって重み付けされた対数同時尤度関数であり、任意の最適化アルゴリズムによって最適化することが可能である。例えば、上述の例における正規分布では、パラメータ最適化処理部１０７は、平均場近似法によって、パラメータを最適化することができる。また、パラメータ最適化処理部１０７は、最適化されたパラメータに対して最適化基準Aを同時に計算する。この際に、パラメータ最適化処理部１０７は、上述の情報量基準近似処理部１０５の近似計算を利用する。すなわち、パラメータ最適化処理部１０７は、式９によるヘッセ行列の行列式の近似結果を利用する。
【００６３】
最適性判定処理部１０８は、最適化基準Aの収束を判定する。収束の判定は、最適化基準Aの絶対変化量や相対変化量などに閾値を設け、閾値を用いることで実現できる。
【００６４】
モデル推定結果出力装置１０９は、最適な隠れ状態数、観測確率のパラメータ、変分分布などをモデル推定結果出力結果１１２として出力する。
【００６５】
隠れ状態数設定部１０２、初期化処理部１０３、隠れ変数変分確率計算処理部１０４、情報量基準近似処理部１０５、隠れ状態選択処理部１０６、パラメータ最適化処理部１０７および最適性判定処理部１０８は、例えば、因子化隠れマルコフモデル推定プログラムに従って動作するコンピュータのＣＰＵによって実現される。この場合、ＣＰＵが因子化隠れマルコフモデル推定プログラムを読み込み、そのプログラムに従って、隠れ状態数設定部１０２、初期化処理部１０３、隠れ変数変分確率計算処理部１０４、情報量基準近似処理部１０５、隠れ状態選択処理部１０６、パラメータ最適化処理部１０７および最適性判定処理部１０８として動作すればよい。因子化隠れマルコフモデル推定プログラムは、コンピュータが読み取り可能な記録媒体に記憶されていてもよい。また、上記の各要素１０２〜１０８がそれぞれ別々のハードウェアで実現されていてもよい。
【００６６】
図２は、本発明の処理経過の例を示すフローチャートである。データ入力装置１０１を介して入力データ１１１が入力される（ステップＳ１００）。
【００６７】
次に、隠れ状態数設定部１０２は、隠れ状態数の初期値として、入力データ１１１として入力された隠れ状態数の最大値を設定する（ステップＳ１０１）。すなわち、隠れ状態数設定部１０２は、モデルのlayer1隠れ状態数Kを、入力された最大値Kmaxに設定する。また、隠れ状態数設定部１０２は、layer2隠れ状態数Mkを、入力された最大値Mmaxに設定する。
【００６８】
次に、初期化処理部１０３は、指定された隠れ状態数に対して、推定のためパラメータ（例えば、各観測確率のパラメータj）や隠れ変数変分確率の初期化処理を実施する（ステップＳ１０２）。
【００６９】
次に、情報量基準近似処理部１０５は、ヘッセ行列の行列式の近似処理を行う（ステップＳ１０３）。情報量基準近似処理部１０５は、式９の計算によりヘッセ行列の行列式の近似値を計算する。
【００７０】
隠れ変数変分確率計算処理部１０４は、計算されたヘッセ行列の行列式の近似値を利用して、隠れ変数の変分確率を計算する（ステップＳ１０４）。
【００７１】
次に、隠れ状態選択処理部１０６は、前述の閾値判定により、不要な隠れ状態をモデルから除去する（ステップＳ１０５）。すなわち、隠れ状態選択処理部１０６は、k番目の隠れ状態について、
【数２３】
0006398990.tif 000024
が、入力データ１１１として設定された閾値を下回った場合には、その状態のモデルから、layer2のm番目の隠れ状態を削除する。また、隠れ状態選択処理部１０６は、layer2の隠れ状態数が1になった場合には、layer1の隠れ状態を削除する。
【００７２】
次に、パラメータ最適化処理部１０７は、最適化基準Aを最適とするパラメータを算出する（ステップＳ１０６）。なお、パラメータ最適化処理部１０７が最初にステップＳ１０６を実行する場合に用いる最適化基準Aは、例えば、初期化処理部１０３によってランダムに定められていてもよい。あるいは、初期化処理部１０３が隠れ変数の変分確率をランダムに設定し、１回目のステップＳ１０３〜Ｓ１０９ａのループ処理（図２参照）では、ステップＳ１０６の処理を実行しなくてもよい。
【００７３】
次に、情報量基準近似処理部１０５は、ヘッセ行列の行列式の近似処理を行う（ステップＳ１０７）。情報量基準近似処理部１０５は、式９の計算によりヘッセ行列の行列式の近似値を計算する。
【００７４】
次に、パラメータ最適化処理部１０７は、ステップＳ１０６で最適化されたパラメータを利用し、最適化基準Aの値を算出する（ステップＳ１０８）。
【００７５】
次に、最適性判定処理部１０８は、最適化基準A が収束したか否かを判定する（ステップＳ１０９）。例えば、最適性判定処理部１０８は、直近のステップＳ１０３〜Ｓ１０９ａのループ処理で得られた最適化基準Aと、その１つ前のステップＳ１０３〜Ｓ１０９ａのループ処理で得られた最適化基準Aとの差を計算し、その差の絶対値が予め定められた閾値以下になっていれば、最適化基準Aが収束したと判定してよい。また、最適性判定処理部１０８は、その差の絶対値が閾値より大きければ、最適化基準Aが収束していないと判定してよい。
【００７６】
最適化基準Aが収束していないと判定した場合（ステップＳ１０９ａのＮｏ）、因子化隠れマルコフモデル推定装置１００は、ステップＳ１０３以降の処理を繰り返す。最適化基準Aが収束したと判定した場合（ステップＳ１０９ａのＹｅｓ）、モデル推定結果出力装置１０９は、モデル推定結果を出力し、処理を完了する（ステップＳ１１０）。ステップＳ１１０では、モデル推定結果出力装置１０９は、最適化基準Aが収束したと判定された時点における隠れ状態数、およびその時点で得られているパラメータ、変分分布を出力する。
【００７７】
以下では、本発明で提案する因子化隠れ変数モデル推定装置の応用例について、屋内に設置された位置センサから人間の行動を推定する場合を例に説明する。本例では、xとして、D次元のセンサ反応時系列を考える。このようなデータから、人間の行動（例えば、食事中、睡眠中、外出中など）を推定する場合に、屋内にいる人間が１名であれば、隠れマルコフモデルの隠れ状態として人間の行動を抽出可能である。一方、屋内にいる人間が複数存在するのであれば、各人の位置が同時に観測される。そのため、因子化隠れマルコフモデルによるモデル化が適切である。ここで、因子化隠れマルコフモデルのlayer1隠れ状態が人間に相当し、layer2隠れ状態が行動に相当する。そして、本発明によって因子化隠れマルコフモデルを推定することで、行動中の人間の数（layer1の隠れ状態数）及び、各人間の行動（layer2の隠れ状態）をデータから自動的に抽出することが可能である。
【００７８】
次に、本発明の概要について説明する。図３は、本発明の概要を示すブロック図である。本発明の因子化隠れマルコフモデル推定装置１００は、近似値計算部７１と、変分確率計算部７２と、隠れ状態削除部７３と、パラメータ最適化部７４と、収束判定部７５とを備える。
【００７９】
近似値計算部７１（例えば、情報量基準近似処理部１０５）は、因子化隠れマルコフモデルの各layer1隠れ変数によって決まるパラメータの線形結合として表現された観測モデルのパラメータに関するヘッセ行列の行列式の近似値を計算する（例えば、式９の近似計算を行う）。
【００８０】
変分確率計算部７２（例えば、隠れ変数変分確率計算処理部１０４）は、その行列式の近似値を利用して隠れ変数の変分確率を計算する。
【００８１】
隠れ状態削除部７３（例えば、隠れ状態選択処理部１０６）は、変分分布に基づいて隠れ状態を削除する。
【００８２】
パラメータ最適化部７４（例えば、パラメータ最適化処理部１０７）は、周辺化対数尤度関数を完全変数に対する推定量に関してラプラス近似した近似量の下界として定義される基準値（例えば、最適化基準A）に対してパラメータを最適化するとともに、基準値を計算する。
【００８３】
収束判定部７５（例えば、最適性判定処理部１０８）は、基準値が収束したか否かを判定する。
【００８４】
また、近似値計算部７１がヘッセ行列の行列式の近似値を計算し、変分確率計算部７２が隠れ変数の変分確率を計算し、隠れ状態削除部７３が隠れ状態を削除し、パラメータ最適化部７４がパラメータを最適化し、近似値計算部７１がヘッセ行列の行列式の近似値を計算し、パラメータ最適化部７４が基準値を計算し、収束判定部７５が、基準値が収束したか否かを判定するループ処理を、収束判定部７５が、基準値が収束したと判定するまで繰り返す構成であることが好ましい。
【００８５】
この出願は、２０１３年５月２０日に出願された米国特許出願１３／８９８０８７を基礎とする優先権を主張し、その開示の全てをここに取り込む。
【００８６】
以上、実施形態を参照して本願発明を説明したが、本願発明は上記の実施形態に限定されるものではない。本願発明の構成や詳細には、本願発明のスコープ内で当業者が理解し得る様々な変更をすることができる。
【符号の説明】
【００８７】
１０１  データ入力装置
１０２  隠れ状態数設定部
１０３  初期化処理部
１０４  隠れ変数変分確率計算処理部
１０５  情報量基準近似処理部
１０６  隠れ状態選択処理部
１０７  パラメータ最適化処理部
１０８  最適性判定処理部
１０９  モデル推定結果出力装置
【図１】
0006398990.tif 000025
【図２】
0006398990.tif 000026
【図３】
0006398990.tif 000027