(57)【要約】【課題】計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができるようにする。【解決手段】パラメタ推定部２０が、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、観測データの確率密度関数を用いて表される、当該打ち切りデータの確率分布を表すモデルと、当該打ち切りデータから得られる、当該打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、モデルのパラメタを推定する。【選択図】図６
2020135554.tif 000002
【特許請求の範囲】
【請求項１】
観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する推定装置であって、
前記打ち切りデータの入力を受け付ける入力部と、
前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定するパラメタ推定部と、
を含む推定装置。
【請求項２】
前記打ち切りデータの確率分布を表すモデルは、
前記観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、前記変数の確率分布と、
前記観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、前記変数が与えられた下での、前記観測データの確率分布と、
を用いて表される請求項１記載の推定装置。
【請求項３】
前記目的関数は、前記打ち切りデータの確率分布を表すモデルと、前記打ち切りデータの確率分布とのカルバックライブラーダイバージェンス、又はＬ２ダイバージェンスである
請求項１又は２記載の推定装置。
【請求項４】
観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する推定方法であって、
入力部が、前記打ち切りデータの入力を受け付け、
パラメタ推定部が、前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定する
推定方法。
【請求項５】
観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する処理をコンピュータに実行させるプログラムであって、
入力部が、前記打ち切りデータの入力を受け付け、
パラメタ推定部が、前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定する
ことを含む処理をコンピュータに実行させるプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、推定装置、推定方法、及びプログラムに係り、特に、打ち切りデータから混合モデルのパラメタを推定するための推定装置、推定方法、及びプログラムに関する。
【背景技術】
【０００２】
打ち切りデータとは、観測値がある閾値以上（またはある閾値以下）であるサンプルについては、値が観測されず、閾値以上である、という情報しか得られないデータのことを指す。病気の発症や人の死亡などを記述する臨床データや、インターネット回線利用者の契約履歴データ、Ｅコマースサイトのサービス利用履歴データなど多くのデータが打ち切りデータとして表現される。打ち切りデータを用いる問題の代表例は、例えばある機器が故障するまでに要する時間の分布を推定する、生存時間分析である。生存時間分布は、初期故障や劣化故障などの存在によりしばしば多峰性をもつため、混合モデルを用いた分布推定が広く用いられる。
【０００３】
打ち切りデータから混合モデルのパラメタを推定するには、文献（非特許文献１）で提案されているＥＭＣＭ（Ｅｘｐｅｃｔａｔｉｏｎ  Ｍａｘｉｍｉｚａｔｉｏｎ  ｆｏｒ  Ｃｅｎｓｏｒｅｄ  Ｍｉｘｔｕｒｅｍｏｄｅｌｓ）アルゴリズムが適用できる。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】Didier Chauveau. A stochastic em algorithm for mixtures with censored data. Journal of statistical planning and inference, Vol.46, No.1, 1995, pp.1-25.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかし、このアプローチには以下の２つの問題があった。
【０００６】
１つ目は、局所最適解の存在である。すなわち、ＥＭＣＭは目的関数の単調減少は保証されるものの収束先は初期値に依存して変わるため、異なる初期値からの繰り返し実行が必要となってしまう、という問題があった。
【０００７】
２つ目は、モデルで利用する確率分布の切断分布（ある一定の範囲でしか値をとらないように変更した分布）の統計量の計算の必要性である。すなわち、１次元正規分布の切断分布である１次元切断正規分布などの例外的なものを除き、この統計量は解析的に計算することができず、モンテカルロ法などの数値計算を利用することが必要となる、という問題があった。ＥＭＣＭは、パラメタの更新と統計量の計算を何度も繰り返すアルゴリズムであるから、各反復で数値計算を繰り返すことは避けることが望ましいと考えられる。
【０００８】
本発明は上記の点に鑑みてなされたものであり、計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができる推定装置、推定方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００９】
本発明に係る推定装置は、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する推定装置であって、前記打ち切りデータの入力を受け付ける入力部と、前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定するパラメタ推定部と、を備えて構成される。
【００１０】
また、本発明に係る推定方法は、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する推定方法であって、入力部が、前記打ち切りデータの入力を受け付け、パラメタ推定部が、前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定する。
【００１１】
また、本発明に係るプログラムは、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する処理をコンピュータに実行させるプログラムであって、入力部が、前記打ち切りデータの入力を受け付け、パラメタ推定部が、前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定することを含む処理をコンピュータに実行させるプログラムである。
【００１２】
本発明に係る推定装置、推定方法及びプログラムによれば、入力部が、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの入力を受け付ける。
【００１３】
そして、パラメタ推定部が、入力部が受け付けた打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、観測データの確率密度関数を用いて表される、打ち切りデータの確率分布を表すモデルと、入力部が受け付けた打ち切りデータから得られる、打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、モデルのパラメタを推定する。
【００１４】
このように、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、観測データの確率密度関数を用いて表される、当該打ち切りデータの確率分布を表すモデルと、当該打ち切りデータから得られる、当該打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、モデルのパラメタを推定することにより、計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができる。
【００１５】
また、本発明に係る推定装置の前記打ち切りデータの確率分布を表すモデルは、前記観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、前記変数の確率分布と、前記観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、前記変数が与えられた下での、前記観測データの確率分布と、を用いて表されることができる。
【００１６】
また、本発明に係る推定装置の前記目的関数は、前記打ち切りデータの確率分布を表すモデルと、前記打ち切りデータの確率分布とのカルバックライブラーダイバージェンス、又はＬ２ダイバージェンスであることができる。
【発明の効果】
【００１７】
本発明の推定装置、推定方法、及びプログラムによれば、計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができる。
【図面の簡単な説明】
【００１８】
【図１】１次元の打ち切りデータの例を表すイメージ図である。
【図２】２次元の打ち切りデータの例を表すイメージ図である。
【図３】１次元の打ち切りデータの生存時間表現の例を表すイメージ図である。
【図４】２次元の打ち切りデータの生存時間表現の例を表すイメージ図である。
【図５】本発明の実施の形態に係る推定装置として機能するコンピュータの概略構成を示すブロック図である。
【図６】本発明の実施の形態に係る推定装置の構成を示すブロック図である。
【図７】本発明の実施の形態に係る推定装置の推定処理ルーチンを示すフローチャートである。
【発明を実施するための形態】
【００１９】
以下、本発明の実施の形態について図面を用いて説明する。
【００２０】
＜本発明の実施の形態に係る推定装置の原理＞
まず、本発明の実施形態の原理について説明する。
【００２１】
本発明の実施の形態では、次に示す目的関数の異なる２つの手法を構築した。１つ目は確率分布間のＫｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ（ＫＬ：カルバック−ライブラー）ダイバージェンス最小化に基づく推定手法、２つ目は確率分布間のＬ２ダイバージェンス最小化に基づく推定手法である。１つ目の推定手法では非常にシンプルな繰り返し計算でパラメタの推定が可能であり、２つ目の推定手法では、繰り返しの必要すらなく解析的にパラメタの推定が可能である。
【００２２】
本発明の実施の形態に係る手法を構築する上でポイントとなるのは、ｅｘａｍｐｌａｒ  ｂａｓｅｄ  ｍｏｄｅｌ（以下、ｅｍｂ）と呼ばれるアプローチ（参考文献１）の利用である。
［参考文献１］Danial Lashkari and Polina Golland, "Convex clustering with exemplar-based models", In Advances in neural information processing systems, 2008, pp.825-832.
【００２３】
ｅｍｂアプローチでは、混合モデルの各コンポーネントのパラメタを陽にパラメタとしては扱わず、データ点の存在する点に各コンポーネントを配置する。これによって、混合比（各コンポーネントの重みパラメタ）のみの推定を行う定式化をすることで大域的最適解へ収束するアルゴリズムが構築できる。上記の手法は、このアプローチに対して打ち切りデータを入力として利用できるよう発展させたものとみなすことができる。
【００２４】
更に、上述したＥＭＣＭで必要となる切断分布の統計量はコンポーネントのパラメタの推定に必要なものであったため、このアプローチの採用により数値計算の繰り返し実行の必要のない手法を構築することができる。
【００２５】
なお、文献（参考文献１）の方法では、通常のデータを用いてＫＬダイバージェンス最小化によるｅｍｂを行う手法を提案しているが、打ち切りデータを扱うことができない。
【００２６】
確率分布の誤差の関数にＬ２ダイバージェンスを用いるアプローチは文献（参考文献２）でも考えられているが、打ち切りデータではなく通常のデータを入力する状況を考えている。
［参考文献２］DavidW Scott, "Parametric statistical modeling by minimum integrated square error", Technometrics, Vol.43, No.3, 2001, pp.274-285.
【００２７】
また、打ち切りデータを利用し、かつＬ２ダイバージェンスを用いる研究には文献（参考文献３）が存在するが、指数分布、ワイブル分布のように単峰の単純な分布をモデルとして用いる場合のみが考えられている。
［参考文献３］Srabashi Basu, Ayanendranath Basu, and MCJones, "Robust and efficient parametric estimation for censored survival data", Annals of the Institute of Statistical Mathematics, Vol.58, No.2, 2006, pp.341-355.
【００２８】
＜＜準備＞＞
＜＜＜打ち切りデータ＞＞＞
まず、打ち切りデータについて説明する。図１に、１次元の打ち切りデータの例を示す。図１の例に示すように、打ち切りデータの代表例である機器故障データを用いて説明する。
【００２９】
機器故障データでは、各機器の設置された時刻と故障が起こった際の故障の時刻が記録されている。機器１及び２に関しては観察期間中に故障が起きているため、その時刻が記録されている。一方、機器３及び４に関しては観察期間中には故障せずに、観察が打ち切られたため、故障の時刻が記録されてはいない。
【００３０】
しかしながら、機器３及び機器４の何れも、いつかは必ず故障するとして、故障の時刻は観察終了時刻以降にあるということは読み取れる。このように機器１と２のように観測値（故障の時刻）が分かるデータと機器３と４のように観測値（故障の時刻）がある値以上であることが分かるデータの組からなるデータを、打ち切りデータと呼ぶ。
【００３１】
図１の機器故障データは、１次元の打ち切りデータであるが、本手法では２次元以上の打ち切りデータも扱うことができるためここで説明する。
【００３２】
図２は、ある２つのサービスの利用者の利用期間を表す２次元の打ち切りデータである。図２のケースでは、利用者の（少なくとも１つの）サービスを利用開始した時刻と各サービスを追加契約または解約した時刻が記録されている。
【００３３】
利用者１は両サービスを同時に利用開始し、観察期間中に同時に解約しているため両方のサービスの解約時刻が記録されている。利用者２は両サービスを同時に利用開始し、観察期間中にサービス２のみを解約している。利用者３は両サービスを同時に利用開始し、観察期間中にサービス１のみを解約している。利用者４はサービス２をまず利用開始し、観察期間中にサービス１を追加契約している。よって、観測の打ち切りによって、利用者２のサービス１の解約時刻、利用者３のサービス２の解約時刻、利用者４のサービス１と２の解約時刻は記録されていない。
【００３４】
このように２次元の打ち切りデータにおいては、打ち切られている値の次元が異なる、３通りの打ち切りの種類が存在する。一般にｎ次元の打ち切りデータでは２ｎ−１種類の打ち切りが存在する。なお、以下、各次元に観測が打ち切られるかが決まる例を用いて説明するが、どれか一つの要素が打ち切られたときに全ての要素が観測されなくなるという状況でも同様のアプローチで対応することができる。
【００３５】
ここで、打ち切りデータの定義を与える。扱いを簡単にするため、観測データは図１のようなカレンダー時刻ではなく、生存時間（機器の設置から故障までの時間、サービスの契約から解約までの時間）を用いて表現する。図１の観測データを生存時間で表現したものを図３、図２の観測データを生存時間で表現したものを図４に示す。多次元の打ち切りデータとして定義するため、以下では図４を説明の例に用いる。
【００３６】
打ち切りデータを
2020135554.tif 000003
と書く。ここで、
2020135554.tif 000004
と、
2020135554.tif 000005
とは共にｄｘ次元のベクトルであり、
2020135554.tif 000006
、
2020135554.tif 000007
である。ｘｉｊが利用者ｉのサービスｊの利用時間、ｗｉｊが利用者ｉのサービスｊの解約時刻が記録されたか（ｗｉｊ＝１）、打ち切りにより記録されなかったか（ｗｉｊ＝０）を表す。同様に、ｉ番目の利用者の、観測終了時刻までの時間の長さを
2020135554.tif 000008
と書く。ｖｉｊが利用者ｉのサービスｊ利用開始時刻から観測終了時刻までの長さを表す。打ち切りにより観測値が観測されなかったとき（ｗｉｊ＝０）には、ｘｉｊ＝ｖｉｊと設定されているとする。
【００３７】
＜＜＜混合モデル＞＞＞
次に、本発明の実施の形態で用いるモデルについて説明する。混合モデルで表される、観測値の確率密度関数は一般に下記式（１）で定義される。
【００３８】
【数１】
2020135554.tif 000009
【００３９】
ここで、Ｋは混合数、
2020135554.tif 000010
がｋ番目コンポーネントの確率分布を表す。コンポーネントの確率分布
2020135554.tif 000011
には、例えば下記式で表されるガウス分布が利用できる。
【００４０】
【数２】
2020135554.tif 000012
【００４１】
ここで、
2020135554.tif 000013
とσとは、ガウス分布の平均と標準偏差を表す。ただし、ｅｂｍ（参考文献１）のアプローチに従い、コンポーネントの確率分布のパラメタ
2020135554.tif 000014
はＫ＝ｎとして、それぞれが観測データ点に対応するよう
2020135554.tif 000015
と設定されたもの等であるとする。この場合には、観測値の確率密度関数は、観測データの各々に対する、各コンポーネントの混合モデルで表され、各コンポーネントに含まれる、ガウス分布の平均を、対応する観測値とする。本手法は打ち切りデータを扱うため、Ｋ＝ｎとし、値が観測されていれば（ｗｉｊ＝１）、μｉｊ＝ｘｉｊ、そうでなければ（ｗｉｊ＝０）、μｉｊ＝ｘｉｊ＋εと設定してもよい。εは０以上の値を取る確率分布（例えば指数分布）からランダムに生成した値を表す。また、データ数が多い場合には、例えばランダムに選んだ１００個のデータのみを用いてもよいし、事前知識に基づいて設定したコンポーネントを用いてもよい。標準偏差σは交差検証法等により決定することが可能である。
【００４２】
上記モデルを用いたときの打ち切りデータ
2020135554.tif 000016
の生成過程は次のように記述できる。まず、各サンプルについて観測終了までの時間の長さ
2020135554.tif 000017
が既知のもと、打ち切りが起こるか否かを表す変数
2020135554.tif 000018
が、下記式（２）の確率分布に従い生成される。
【００４３】
【数３】
2020135554.tif 000019
【００４４】
ただし、
2020135554.tif 000020
と、
2020135554.tif 000021
とのうち、ｗｉｊ＝１である、観測値が観測された要素の集合を
2020135554.tif 000022
とし、ｗｉｊ＝０である、観測値が観測されなかった集合を
2020135554.tif 000023
とする。また、観測値が観測された要素と、観測値が観測されなかった要素とを区別しない場合に、観測データと呼ぶこととする。コンポーネントの確率分布
2020135554.tif 000024
に、上記式（２）のガウス分布のように標準的な確率分布を利用すると、
2020135554.tif 000025
を累積密度関数を用いて解析的に計算できる。
【００４５】
2020135554.tif 000026
であり、少なくとも一つ観測された要素がある場合、
2020135554.tif 000027
は下記式（４）で表される分布にしたがって生成される。
【００４６】
【数４】
2020135554.tif 000028
【００４７】
ここで、
2020135554.tif 000029
はデルタ関数であり、ｆｔｒはｆを観測されなかった要素に関して周辺化した分布の切断分布であり、下記式で表される。
【００４８】
【数５】
2020135554.tif 000030
【００４９】
ただし、ｆｔｒは下記式（５）〜（６）に従う。
【００５０】
【数６】
2020135554.tif 000031
【００５１】
コンポーネントが上記式（２）のガウス分布である場合、
2020135554.tif 000032
は
2020135554.tif 000033
となる。ただし、
2020135554.tif 000034
はそれぞれ
2020135554.tif 000035
から
2020135554.tif 000036
と
2020135554.tif 000037
に対応する次元の要素を抜き出したベクトルである。
【００５２】
2020135554.tif 000038
であり、一つも観測された要素がない場合、下記式（７）で示すように、デルタ関数のみの表現となる。
【００５３】
【数７】
2020135554.tif 000039
【００５４】
よって、上記をまとめると、各打ち切りデータの生成確率は、下記式（８）で与えられる。
【００５５】
【数８】
2020135554.tif 000040
【００５６】
＜＜ＫＬダイバージェンスとＬ２ダイバージェンス＞＞
次に、提案手法の目的関数を定義する際に利用するダイバージェンスについて記す。良く知られるように、確率分布ｐ（ｘ）とｑ（ｘ）とに対するカルバックライブラー（ＫＬ）ダイバージェンスは下記式（９）で定義される。
【００５７】
【数９】
2020135554.tif 000041
【００５８】
これに加え、本発明の実施の形態では、下記で定義されるＬ２ダイバージェンス（参考文献２）も利用する場合についても説明する（下記式（１０））。
【００５９】
【数１０】
2020135554.tif 000042
【００６０】
Ｌ２ダイバージェンスは、２つの確率密度関数の２乗誤差として定義されている。何れのダイバージェンスを用いるべきかについては、問題に応じて異なる。よって、本発明の実施の形態では、何れのダイバージェンスを利用する場合にも適用できるように、２種類の手法を構築することとした。
【００６１】
＜＜ＫＬダイバージェンスの最適化による推定＞＞
まず、目的関数としてＫＬダイバージェンスを利用する場合の提案手法を示す。打ち切りデータの確率分布を表すモデル
2020135554.tif 000043
と、打ち切りデータから得られる真の確率分布
2020135554.tif 000044
とのＫＬダイバージェンスは、上記式（９）の定義に従い、下記式（１１）で与えられる。
【００６２】
【数１１】
2020135554.tif 000045
【００６３】
これは、下記式（１２）のように式変形することができる。
【００６４】
【数１２】
2020135554.tif 000046
【００６５】
定数項を除去し、未知である真の分布
2020135554.tif 000047
に関する期待値を標本平均で置き換えれば、下記の式（１３）が導かれる。
【００６６】
【数１３】
2020135554.tif 000048
【００６７】
ここで、
2020135554.tif 000049
は、データから計算できる量であり、これを目的関数とすることでアルゴリズムを導く。具体的には、下記式（１４）の最適化問題を解けばよい。
【００６８】
【数１４】
2020135554.tif 000050
【００６９】
ここで、制約条件（パラメタの要素が０以上、かつ、和が１）は、混合モデルｆが確率分布となるためのものである。ラグランジュの未定乗数法を用いると、上記の最適化問題の解は、下記式（１５）を満たすことが分かる。
【００７０】
【数１５】
2020135554.tif 000051
【００７１】
よって、上記式（１５）に基づいて
2020135554.tif 000052
の更新を繰り返すことにより、最適化が可能となる。なお、参考文献１の手法と同様に、計算量の削減と収束を早めるため、パラメタ更新の際にθｋがある閾値（例えば１０−３／ｎ）より小さい場合には、θｋ＝０と設定した後に、全体を和が１になるように調整する、再正規化操作を行ってもよい。
【００７２】
＜＜Ｌ２ダイバージェンスの最適化による推定＞＞
次に、Ｌ２ダイバージェンスを利用する場合の提案手法を示す。ＫＬダイバージェンスとは異なり、Ｌ２ダイバージェンスの定義から直接目的関数を定義することはせず、ＫＬダイバージェンスを用いた際の目的関数に注目することにより新たな目的関数を定義する。
【００７３】
ＫＬダイバージェンスを用いた際の目的関数（式（１２））に注目すると、混合モデルｆの周辺分布
2020135554.tif 000053
と、変数が与えられた下での観測値の真の分布
2020135554.tif 000054
のＫＬダイバージェンスに対応する項と、観測値が観測されなかったことを表す変数のモデルの分布
2020135554.tif 000055
と、真の分布
2020135554.tif 000056
を用いた対数尤度比に対応する項との２つの項をそれぞれ
2020135554.tif 000057
で重み付き和を取ることで構成されていることが分かる。
【００７４】
この洞察に基づき、これら２つの各項でＫＬダイバージェンス／対数尤度比を用いている箇所をＬ２ダイバージェンスに置き換えることで、目的関数を下記式のように設計することができる。
【００７５】
【数１６】
2020135554.tif 000058
【００７６】
これを式変形すれば、
2020135554.tif 000059
となる。更に、定数項（Ｃｏｎｓｔ）を除去し、真の分布に関する平均を標本平均で置き換えた下記式の最適化を考える。
【００７７】
【数１６】
2020135554.tif 000060
【００７８】
ただし、
2020135554.tif 000061
とおいた。これを目的関数とすることで、アルゴリズムを導出する。
2020135554.tif 000062
は、行列・ベクトル形式で、下記式（１６）、（１７）のように表すことができる。
【００７９】
【数１７】
2020135554.tif 000063
【００８０】
ただし、下記式（１８）〜（２０）に従う。
【００８１】
【数１８】
2020135554.tif 000064
【００８２】
また、
2020135554.tif 000065
は、コンポーネントの分布にガウス分布を用いる場合、解析的に計算できる値であり、下記式（２１）のように表すことができる。
【００８３】
【数１９】
2020135554.tif 000066
【００８４】
上記式（２１）中で、
2020135554.tif 000067
に関して総和をとる際、
2020135554.tif 000068
や
2020135554.tif 000069
は
2020135554.tif 000070
の値によって異なる値であることを明記しておく。これにより、目的関数は
2020135554.tif 000071
に関する２次の形式で表現されることが分かる。よって、パラメタの推定値は下記で示す制約付きの２次最適化問題を解くことで得ることができる。
【００８５】
【数２０】
2020135554.tif 000072
【００８６】
数値ソルバーを用いて直接上記式（２２）の最適化問題を解くことができる。その際、正則化項を加えた下記の最適化問題を解くことにしてもよい。
【００８７】
【数２１】
2020135554.tif 000073
【００８８】
ただし、βがハイパーパラメタを表す。また、例えば次に示す近似的な方法を用いてもよい。上記式（２２）の最適化問題から制約条件を外した問題の最適解は下記式（２３）、（２４）のように求められる。
【００８９】
【数２２】
2020135554.tif 000074
【００９０】
ただし、
2020135554.tif 000075
は正則化項、βはハイパーパラメタであり、パラメタの発散を防ぐ効果がある。
2020135554.tif 000076
は、値が０以上かつ和が１であるから、下記式（２５）の処理で条件を満たす
2020135554.tif 000077
が得られる。
【００９１】
【数２３】
2020135554.tif 000078
【００９２】
ただし、
2020135554.tif 000079
は、
2020135554.tif 000080
ノルムを表す
【００９３】
本発明の実施の形態では、以上説明した２つの手法の何れかにより、モデルのパラメタを推定することにより、計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができる。
【００９４】
＜本発明の実施の形態に係る推定装置の構成＞
次に、図５及び図６を参照して、本発明の実施の形態に係る推定装置１の構成について説明する。図５は、本発明の実施の形態に係る推定装置１として機能するコンピュータの概略構成を示すブロック図である。図６は、本発明の実施の形態に係る推定装置１の構成を示すブロック図である。
【００９５】
図５に示すように、推定装置１は、ＣＰＵ１１０と、ＲＡＭ等のメモリ１２０と、通信インターフェース（ＩＦ）部１３０と、キーボード等の入力部１４０と、ディスプレイ等の表示部１５０と、後述する推定処理ルーチンを実行するためのプログラム１７０を記憶したＲＯＭ等の記憶部１６０とを備えたコンピュータで構成されている。また、ＣＰＵ１１０、メモリ１２０、通信ＩＦ部１３０、入力部１４０、表示部１５０、及び記憶部１６０は、バス１００を介して接続されている。また、通信ＩＦ部１３０は、ＬＡＮケーブル等の通信回線により外部装置２と接続されている。なお、ネットワーク（図示しない）を介して外部装置２と接続される構成としてもよい。
【００９６】
図６に示すように、本実施形態に係る推定装置１は、データ処理部１０と、パラメタ推定部２０と、パラメタ出力部３０と、記憶部４０と、入力部５０と、出力部６０とを備えて構成される。
【００９７】
データ処理部１０は、入力部５０が受け付けた観測値が観測されたサンプルの観測データ
2020135554.tif 000081
と、観測値が観測されなかったサンプルの観測データ
2020135554.tif 000082
と、各サンプルについて観測値が観測されたか否かを表す変数
2020135554.tif 000083
とを含む打ち切りデータ
2020135554.tif 000084
をデータ記憶部４１に格納する。サンプルは、例えば、上述の図１、３の例では各機器のこと、図２、４の例では各利用者のことである。
【００９８】
パラメタ推定部２０は、入力部５０が受け付けた打ち切りデータ
2020135554.tif 000085
の各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、観測データの確率密度関数を用いて表される、打ち切りデータ
2020135554.tif 000086
の確率分布を表すモデル
2020135554.tif 000087
と、入力部５０が受け付けた打ち切りデータ
2020135554.tif 000088
から得られる、当該打ち切りデータ
2020135554.tif 000089
の真の確率分布
2020135554.tif 000090
とのダイバージェンスである目的関数を最適化することにより、モデルのパラメタ
2020135554.tif 000091
を推定する。
【００９９】
具体的には、パラメタ推定部２０は、まず、当該打ち切りデータ
2020135554.tif 000092
の真の確率分布
2020135554.tif 000093
を求める。
【０１００】
次に、パラメタ推定部２０は、目的関数を、当該打ち切りデータ
2020135554.tif 000094
の確率分布を表すモデル
2020135554.tif 000095
と、当該打ち切りデータ
2020135554.tif 000096
の真の確率分布
2020135554.tif 000097
とのＫＬダイバージェンス、又はＬ２ダイバージェンスとして、パラメタを推定する。
【０１０１】
ここで、打ち切りデータの確率分布を表すモデル
2020135554.tif 000098
は、上記式（８）で示すように、観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、観測データの確率分布
2020135554.tif 000099
と、観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、変数が与えられた下での、変数の確率分布
2020135554.tif 000100
と、を用いて表される。
【０１０２】
パラメタ推定部２０は、ＫＬダイバージェンスを用いる場合、上記式（１５）のパラメタ更新を繰り返すことにより、パラメタを推定する。また、パラメタ推定部２０は、Ｌ２ダイバージェンスを用いる場合、上記式（２４）、（２５）を用いてパラメタを推定する。
【０１０３】
そして、パラメタ推定部２０は、推定したパラメタ
2020135554.tif 000101
をパラメタ記憶部４２に格納する。
【０１０４】
パラメタ出力部３０は、パラメタ記憶部４２のパラメタ
2020135554.tif 000102
を取得し、パラメタ出力部３０に当該パラメタ
2020135554.tif 000103
を渡す。
【０１０５】
記憶部４０は、データ記憶部４１と、パラメタ記憶部４２とを含む。データ記憶部４１には、打ち切りデータ
2020135554.tif 000104
が格納される。また、パラメタ記憶部４２には、モデルのパラメタ
2020135554.tif 000105
が格納される。
【０１０６】
入力部５０は、外部装置２から、打ち切りデータ
2020135554.tif 000106
の入力を受け付ける。そして、入力部５０は、受け付けた打ち切りデータ
2020135554.tif 000107
を、データ処理部１０に渡す。
【０１０７】
出力部６０は、パラメタ出力部３０から受け取ったモデルのパラメタ
2020135554.tif 000108
を、外部装置２へ出力する。
【０１０８】
＜本発明の実施の形態に係る推定装置の作用＞
図７は、本発明の実施の形態に係る推定処理ルーチンを示すフローチャートである。
【０１０９】
入力部５０に打ち切りデータ
2020135554.tif 000109
が入力されると、推定装置１において、図７に示す推定処理ルーチンが実行される。
【０１１０】
まず、ステップＳ１００において、入力部５０は、観測値が観測されたサンプルの観測データ
2020135554.tif 000110
と、観測値が観測されなかったサンプルの観測データ
2020135554.tif 000111
と、各サンプルについて観測値が観測されたか否かを表す変数
2020135554.tif 000112
とを含む打ち切りデータ
2020135554.tif 000113
の入力を受け付ける。また、入力部５０は、各サンプルｉについての観測終了時刻までの時間の長さの入力を受け付ける。
【０１１１】
ステップＳ１１０において、パラメタ推定部２０は、目的関数を、当該打ち切りデータ
2020135554.tif 000114
の確率分布を表すモデル
2020135554.tif 000115
と、当該打ち切りデータ
2020135554.tif 000116
の確率分布
2020135554.tif 000117
とのＫＬダイバージェンス、又はＬ２ダイバージェンスとして、上記式（１５）のパラメタ更新を繰り返すことにより、パラメタを推定するか、あるいは、上記式（２４）、（２５）を用いてパラメタを推定する。
【０１１２】
ステップＳ１２０において、出力部６０は、上記ステップＳ１１０により推定されたパラメタ
2020135554.tif 000118
を出力する。
【０１１３】
以上説明したように、本発明の実施形態に係る推定装置によれば、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、観測データの確率密度関数を用いて表される、当該打ち切りデータの確率分布を表すモデルと、当該打ち切りデータから得られる、当該打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、モデルのパラメタを推定するため、計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができる。
【０１１４】
なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１１５】
例えば、上述の実施の形態では、ダイバージェンスとしてＫＬダイバージェンス又はＬ２ダイバージェンスを用いる場合について説明したが、これに限定されるものではなく、他のダイバージェンスを用いることも可能である。
【０１１６】
また、上述の実施の形態では、時系列データである打ち切りデータを前提に記載しているが、これに限定されるものではなく、時系列データではない任意の打ち切りデータに対しても本発明の適用が可能である。
【０１１７】
また、上述の実施の形態に係る推定装置１は、各部の処理をプログラムとして構築し、推定装置として利用されるコンピュータにインストールして実行させる構成として説明したが、ネットワークを介して流通させる構成としてもよい。
【０１１８】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【０１１９】
１     推定装置
２     外部装置
１０   データ処理部
２０   パラメタ推定部
３０   パラメタ出力部
４０   記憶部
４１   データ記憶部
４２   パラメタ記憶部
５０   入力部
６０   出力部
１００ バス
１１０ ＣＰＵ
１２０ メモリ
１３０ 通信ＩＦ部
１４０ 入力部
１５０ 表示部
１６０ 記憶部
１７０ プログラム
【図１】
2020135554.tif 000119
【図２】
2020135554.tif 000120
【図３】
2020135554.tif 000121
【図４】
2020135554.tif 000122
【図５】
2020135554.tif 000123
【図６】
2020135554.tif 000124
【図７】
2020135554.tif 000125