JP7107246B2 - 推定装置、推定方法、及びプログラム - Google Patents

推定装置、推定方法、及びプログラム Download PDF

Info

Publication number
JP7107246B2
JP7107246B2 JP2019029769A JP2019029769A JP7107246B2 JP 7107246 B2 JP7107246 B2 JP 7107246B2 JP 2019029769 A JP2019029769 A JP 2019029769A JP 2019029769 A JP2019029769 A JP 2019029769A JP 7107246 B2 JP7107246 B2 JP 7107246B2
Authority
JP
Japan
Prior art keywords
observed
data
probability distribution
censored
censored data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2019029769A
Other languages
English (en)
Other versions
JP2020135554A (ja
Inventor
匡宏 幸島
達史 松林
浩之 戸田
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
NTT Inc
NTT Inc USA
Original Assignee
Nippon Telegraph and Telephone Corp
NTT Inc USA
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nippon Telegraph and Telephone Corp, NTT Inc USA filed Critical Nippon Telegraph and Telephone Corp
Priority to JP2019029769A priority Critical patent/JP7107246B2/ja
Priority to PCT/JP2020/004909 priority patent/WO2020170867A1/ja
Priority to US17/431,791 priority patent/US20220138375A1/en
Publication of JP2020135554A publication Critical patent/JP2020135554A/ja
Application granted granted Critical
Publication of JP7107246B2 publication Critical patent/JP7107246B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N99/00—Subject matter not provided for in other groups of this subclass
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Software Systems (AREA)
  • Operations Research (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Evolutionary Biology (AREA)
  • Probability & Statistics with Applications (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Computation (AREA)
  • Computing Systems (AREA)
  • Geometry (AREA)
  • Computer Hardware Design (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Medical Informatics (AREA)
  • Complex Calculations (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Description

本発明は、推定装置、推定方法、及びプログラムに係り、特に、打ち切りデータから混合モデルのパラメタを推定するための推定装置、推定方法、及びプログラムに関する。
打ち切りデータとは、観測値がある閾値以上（またはある閾値以下）であるサンプルについては、値が観測されず、閾値以上である、という情報しか得られないデータのことを指す。病気の発症や人の死亡などを記述する臨床データや、インターネット回線利用者の契約履歴データ、Ｅコマースサイトのサービス利用履歴データなど多くのデータが打ち切りデータとして表現される。打ち切りデータを用いる問題の代表例は、例えばある機器が故障するまでに要する時間の分布を推定する、生存時間分析である。生存時間分布は、初期故障や劣化故障などの存在によりしばしば多峰性をもつため、混合モデルを用いた分布推定が広く用いられる。
打ち切りデータから混合モデルのパラメタを推定するには、文献（非特許文献１）で提案されているＥＭＣＭ（Ｅｘｐｅｃｔａｔｉｏｎ Ｍａｘｉｍｉｚａｔｉｏｎ ｆｏｒ Ｃｅｎｓｏｒｅｄ Ｍｉｘｔｕｒｅｍｏｄｅｌｓ）アルゴリズムが適用できる。
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.
＜＜準備＞＞
＜＜＜打ち切りデータ＞＞＞
まず、打ち切りデータについて説明する。図１に、１次元の打ち切りデータの例を示す。図１の例に示すように、打ち切りデータの代表例である機器故障データを用いて説明する。
機器故障データでは、各機器の設置された時刻と故障が起こった際の故障の時刻が記録されている。機器１及び２に関しては観察期間中に故障が起きているため、その時刻が記録されている。一方、機器３及び４に関しては観察期間中には故障せずに、観察が打ち切られたため、故障の時刻が記録されてはいない。
しかしながら、機器３及び機器４の何れも、いつかは必ず故障するとして、故障の時刻は観察終了時刻以降にあるということは読み取れる。このように機器１と２のように観測値（故障の時刻）が分かるデータと機器３と４のように観測値（故障の時刻）がある値以上であることが分かるデータの組からなるデータを、打ち切りデータと呼ぶ。
図１の機器故障データは、１次元の打ち切りデータであるが、本手法では２次元以上の打ち切りデータも扱うことができるためここで説明する。
図２は、ある２つのサービスの利用者の利用期間を表す２次元の打ち切りデータである。図２のケースでは、利用者の（少なくとも１つの）サービスを利用開始した時刻と各サービスを追加契約または解約した時刻が記録されている。
利用者１は両サービスを同時に利用開始し、観察期間中に同時に解約しているため両方のサービスの解約時刻が記録されている。利用者２は両サービスを同時に利用開始し、観察期間中にサービス２のみを解約している。利用者３は両サービスを同時に利用開始し、観察期間中にサービス１のみを解約している。利用者４はサービス２をまず利用開始し、観察期間中にサービス１を追加契約している。よって、観測の打ち切りによって、利用者２のサービス１の解約時刻、利用者３のサービス２の解約時刻、利用者４のサービス１と２の解約時刻は記録されていない。
このように２次元の打ち切りデータにおいては、打ち切られている値の次元が異なる、３通りの打ち切りの種類が存在する。一般にｎ次元の打ち切りデータでは２ｎ－１種類の打ち切りが存在する。なお、以下、各次元に観測が打ち切られるかが決まる例を用いて説明するが、どれか一つの要素が打ち切られたときに全ての要素が観測されなくなるという状況でも同様のアプローチで対応することができる。
ここで、打ち切りデータの定義を与える。扱いを簡単にするため、観測データは図１のようなカレンダー時刻ではなく、生存時間（機器の設置から故障までの時間、サービスの契約から解約までの時間）を用いて表現する。図１の観測データを生存時間で表現したものを図３、図２の観測データを生存時間で表現したものを図４に示す。多次元の打ち切りデータとして定義するため、以下では図４を説明の例に用いる。
打ち切りデータを
Figure 0007107246000001

と書く。ここで、
Figure 0007107246000002

と、
Figure 0007107246000003

とは共にｄｘ次元のベクトルであり、
Figure 0007107246000004

、
Figure 0007107246000005

である。ｘｉｊが利用者ｉのサービスｊの利用時間、ｗｉｊが利用者ｉのサービスｊの解約時刻が記録されたか（ｗｉｊ＝１）、打ち切りにより記録されなかったか（ｗｉｊ＝０）を表す。同様に、ｉ番目の利用者の、観測終了時刻までの時間の長さを
Figure 0007107246000006

と書く。ｖｉｊが利用者ｉのサービスｊ利用開始時刻から観測終了時刻までの長さを表す。打ち切りにより観測値が観測されなかったとき（ｗｉｊ＝０）には、ｘｉｊ＝ｖｉｊと設定されているとする。
＜＜＜混合モデル＞＞＞
次に、本発明の実施の形態で用いるモデルについて説明する。混合モデルで表される、観測値の確率密度関数は一般に下記式（１）で定義される。
Figure 0007107246000007
ここで、Ｋは混合数、
Figure 0007107246000008

がｋ番目コンポーネントの確率分布を表す。コンポーネントの確率分布
Figure 0007107246000009

には、例えば下記式で表されるガウス分布が利用できる。
Figure 0007107246000010
ここで、
Figure 0007107246000011

とσとは、ガウス分布の平均と標準偏差を表す。ただし、ｅｂｍ（参考文献１）のアプローチに従い、コンポーネントの確率分布のパラメタ
Figure 0007107246000012

はＫ＝ｎとして、それぞれが観測データ点に対応するよう
Figure 0007107246000013

と設定されたもの等であるとする。この場合には、観測値の確率密度関数は、観測データの各々に対する、各コンポーネントの混合モデルで表され、各コンポーネントに含まれる、ガウス分布の平均を、対応する観測値とする。本手法は打ち切りデータを扱うため、Ｋ＝ｎとし、値が観測されていれば（ｗｉｊ＝１）、μｉｊ＝ｘｉｊ、そうでなければ（ｗｉｊ＝０）、μｉｊ＝ｘｉｊ＋εと設定してもよい。εは０以上の値を取る確率分布（例えば指数分布）からランダムに生成した値を表す。また、データ数が多い場合には、例えばランダムに選んだ１００個のデータのみを用いてもよいし、事前知識に基づいて設定したコンポーネントを用いてもよい。標準偏差σは交差検証法等により決定することが可能である。
上記モデルを用いたときの打ち切りデータ
Figure 0007107246000014

の生成過程は次のように記述できる。まず、各サンプルについて観測終了までの時間の長さ
Figure 0007107246000015

が既知のもと、打ち切りが起こるか否かを表す変数
Figure 0007107246000016

が、下記式（２）の確率分布に従い生成される。
Figure 0007107246000017
ただし、
Figure 0007107246000018

と、
Figure 0007107246000019

とのうち、ｗｉｊ＝１である、観測値が観測された要素の集合を
Figure 0007107246000020

とし、ｗｉｊ＝０である、観測値が観測されなかった集合を
Figure 0007107246000021

とする。また、観測値が観測された要素と、観測値が観測されなかった要素とを区別しない場合に、観測データと呼ぶこととする。コンポーネントの確率分布
Figure 0007107246000022

に、上記式（２）のガウス分布のように標準的な確率分布を利用すると、
Figure 0007107246000023

を累積密度関数を用いて解析的に計算できる。
Figure 0007107246000024

であり、少なくとも一つ観測された要素がある場合、
Figure 0007107246000025

は下記式（４）で表される分布にしたがって生成される。
Figure 0007107246000026
ここで、
Figure 0007107246000027

はデルタ関数であり、ｆｔｒはｆを観測されなかった要素に関して周辺化した分布の切断分布であり、下記式で表される。
Figure 0007107246000028
ただし、ｆｔｒは下記式（５）～（６）に従う。
Figure 0007107246000029
コンポーネントが上記式（２）のガウス分布である場合、
Figure 0007107246000030

は
Figure 0007107246000031

となる。ただし、
Figure 0007107246000032

はそれぞれ
Figure 0007107246000033

から
Figure 0007107246000034

と
Figure 0007107246000035

に対応する次元の要素を抜き出したベクトルである。
Figure 0007107246000036

であり、一つも観測された要素がない場合、下記式（７）で示すように、デルタ関数のみの表現となる。
Figure 0007107246000037
よって、上記をまとめると、各打ち切りデータの生成確率は、下記式（８）で与えられる。
Figure 0007107246000038
＜＜ＫＬダイバージェンスとＬ２ダイバージェンス＞＞
次に、提案手法の目的関数を定義する際に利用するダイバージェンスについて記す。良く知られるように、確率分布ｐ（ｘ）とｑ（ｘ）とに対するカルバックライブラー（ＫＬ）ダイバージェンスは下記式（９）で定義される。
Figure 0007107246000039
これに加え、本発明の実施の形態では、下記で定義されるＬ２ダイバージェンス（参考文献２）も利用する場合についても説明する（下記式（１０））。
Figure 0007107246000040
Ｌ２ダイバージェンスは、２つの確率密度関数の２乗誤差として定義されている。何れのダイバージェンスを用いるべきかについては、問題に応じて異なる。よって、本発明の実施の形態では、何れのダイバージェンスを利用する場合にも適用できるように、２種類の手法を構築することとした。
＜＜ＫＬダイバージェンスの最適化による推定＞＞
まず、目的関数としてＫＬダイバージェンスを利用する場合の提案手法を示す。打ち切りデータの確率分布を表すモデル
Figure 0007107246000041

と、打ち切りデータから得られる真の確率分布
Figure 0007107246000042

とのＫＬダイバージェンスは、上記式（９）の定義に従い、下記式（１１）で与えられる。
Figure 0007107246000043
これは、下記式（１２）のように式変形することができる。
Figure 0007107246000044
定数項を除去し、未知である真の分布
Figure 0007107246000045

に関する期待値を標本平均で置き換えれば、下記の式（１３）が導かれる。
Figure 0007107246000046
ここで、
Figure 0007107246000047

は、データから計算できる量であり、これを目的関数とすることでアルゴリズムを導く。具体的には、下記式（１４）の最適化問題を解けばよい。
Figure 0007107246000048
ここで、制約条件（パラメタの要素が０以上、かつ、和が１）は、混合モデルｆが確率分布となるためのものである。ラグランジュの未定乗数法を用いると、上記の最適化問題の解は、下記式（１５）を満たすことが分かる。
Figure 0007107246000049
よって、上記式（１５）に基づいて
Figure 0007107246000050

の更新を繰り返すことにより、最適化が可能となる。なお、参考文献１の手法と同様に、計算量の削減と収束を早めるため、パラメタ更新の際にθｋがある閾値（例えば１０－３／ｎ）より小さい場合には、θｋ＝０と設定した後に、全体を和が１になるように調整する、再正規化操作を行ってもよい。
＜＜Ｌ２ダイバージェンスの最適化による推定＞＞
次に、Ｌ２ダイバージェンスを利用する場合の提案手法を示す。ＫＬダイバージェンスとは異なり、Ｌ２ダイバージェンスの定義から直接目的関数を定義することはせず、ＫＬダイバージェンスを用いた際の目的関数に注目することにより新たな目的関数を定義する。
ＫＬダイバージェンスを用いた際の目的関数（式（１２））に注目すると、混合モデルｆの周辺分布
Figure 0007107246000051

と、変数が与えられた下での観測値の真の分布
Figure 0007107246000052

のＫＬダイバージェンスに対応する項と、観測値が観測されなかったことを表す変数のモデルの分布
Figure 0007107246000053

と、真の分布
Figure 0007107246000054

を用いた対数尤度比に対応する項との２つの項をそれぞれ
Figure 0007107246000055

で重み付き和を取ることで構成されていることが分かる。
この洞察に基づき、これら２つの各項でＫＬダイバージェンス／対数尤度比を用いている箇所をＬ２ダイバージェンスに置き換えることで、目的関数を下記式のように設計することができる。
Figure 0007107246000056
これを式変形すれば、
Figure 0007107246000057

となる。更に、定数項（Ｃｏｎｓｔ）を除去し、真の分布に関する平均を標本平均で置き換えた下記式の最適化を考える。
Figure 0007107246000058
ただし、
Figure 0007107246000059

とおいた。これを目的関数とすることで、アルゴリズムを導出する。
Figure 0007107246000060

は、行列・ベクトル形式で、下記式（１６）、（１７）のように表すことができる。
Figure 0007107246000061
ただし、下記式（１８）～（２０）に従う。
Figure 0007107246000062
また、
Figure 0007107246000063

は、コンポーネントの分布にガウス分布を用いる場合、解析的に計算できる値であり、下記式（２１）のように表すことができる。
Figure 0007107246000064
上記式（２１）中で、
Figure 0007107246000065

に関して総和をとる際、
Figure 0007107246000066

や
Figure 0007107246000067

は
Figure 0007107246000068

の値によって異なる値であることを明記しておく。これにより、目的関数は
Figure 0007107246000069

に関する２次の形式で表現されることが分かる。よって、パラメタの推定値は下記で示す制約付きの２次最適化問題を解くことで得ることができる。
Figure 0007107246000070
数値ソルバーを用いて直接上記式（２２）の最適化問題を解くことができる。その際、正則化項を加えた下記の最適化問題を解くことにしてもよい。
Figure 0007107246000071
ただし、βがハイパーパラメタを表す。また、例えば次に示す近似的な方法を用いてもよい。上記式（２２）の最適化問題から制約条件を外した問題の最適解は下記式（２３）、（２４）のように求められる。
Figure 0007107246000072
ただし、
Figure 0007107246000073

は正則化項、βはハイパーパラメタであり、パラメタの発散を防ぐ効果がある。
Figure 0007107246000074

は、値が０以上かつ和が１であるから、下記式（２５）の処理で条件を満たす
Figure 0007107246000075

が得られる。
Figure 0007107246000076
ただし、
Figure 0007107246000077

は、
Figure 0007107246000078

ノルムを表す
本発明の実施の形態では、以上説明した２つの手法の何れかにより、モデルのパラメタを推定することにより、計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができる。
＜本発明の実施の形態に係る推定装置の構成＞
次に、図５及び図６を参照して、本発明の実施の形態に係る推定装置１の構成について説明する。図５は、本発明の実施の形態に係る推定装置１として機能するコンピュータの概略構成を示すブロック図である。図６は、本発明の実施の形態に係る推定装置１の構成を示すブロック図である。
図５に示すように、推定装置１は、ＣＰＵ１１０と、ＲＡＭ等のメモリ１２０と、通信インターフェース（ＩＦ）部１３０と、キーボード等の入力部１４０と、ディスプレイ等の表示部１５０と、後述する推定処理ルーチンを実行するためのプログラム１７０を記憶したＲＯＭ等の記憶部１６０とを備えたコンピュータで構成されている。また、ＣＰＵ１１０、メモリ１２０、通信ＩＦ部１３０、入力部１４０、表示部１５０、及び記憶部１６０は、バス１００を介して接続されている。また、通信ＩＦ部１３０は、ＬＡＮケーブル等の通信回線により外部装置２と接続されている。なお、ネットワーク（図示しない）を介して外部装置２と接続される構成としてもよい。
図６に示すように、本実施形態に係る推定装置１は、データ処理部１０と、パラメタ推定部２０と、パラメタ出力部３０と、記憶部４０と、入力部５０と、出力部６０とを備えて構成される。
データ処理部１０は、入力部５０が受け付けた観測値が観測されたサンプルの観測データ
Figure 0007107246000079

と、観測値が観測されなかったサンプルの観測データ
Figure 0007107246000080

と、各サンプルについて観測値が観測されたか否かを表す変数
Figure 0007107246000081

とを含む打ち切りデータ
Figure 0007107246000082

をデータ記憶部４１に格納する。サンプルは、例えば、上述の図１、３の例では各機器のこと、図２、４の例では各利用者のことである。
パラメタ推定部２０は、入力部５０が受け付けた打ち切りデータ
Figure 0007107246000083

の各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、観測データの確率密度関数を用いて表される、打ち切りデータ
Figure 0007107246000084

の確率分布を表すモデル
Figure 0007107246000085

と、入力部５０が受け付けた打ち切りデータ
Figure 0007107246000086

から得られる、当該打ち切りデータ
Figure 0007107246000087

の真の確率分布
Figure 0007107246000088

とのダイバージェンスである目的関数を最適化することにより、モデルのパラメタ
Figure 0007107246000089

を推定する。
具体的には、パラメタ推定部２０は、まず、当該打ち切りデータ
Figure 0007107246000090

の真の確率分布
Figure 0007107246000091

を求める。
次に、パラメタ推定部２０は、目的関数を、当該打ち切りデータ
Figure 0007107246000092

の確率分布を表すモデル
Figure 0007107246000093

と、当該打ち切りデータ
Figure 0007107246000094

の真の確率分布
Figure 0007107246000095

とのＫＬダイバージェンス、又はＬ２ダイバージェンスとして、パラメタを推定する。
ここで、打ち切りデータの確率分布を表すモデル
Figure 0007107246000096

は、上記式（８）で示すように、観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、観測データの確率分布
Figure 0007107246000097

と、観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、変数が与えられた下での、変数の確率分布
Figure 0007107246000098

と、を用いて表される。
パラメタ推定部２０は、ＫＬダイバージェンスを用いる場合、上記式（１５）のパラメタ更新を繰り返すことにより、パラメタを推定する。また、パラメタ推定部２０は、Ｌ２ダイバージェンスを用いる場合、上記式（２４）、（２５）を用いてパラメタを推定する。
そして、パラメタ推定部２０は、推定したパラメタ
Figure 0007107246000099

をパラメタ記憶部４２に格納する。
パラメタ出力部３０は、パラメタ記憶部４２のパラメタ
Figure 0007107246000100

を取得し、パラメタ出力部３０に当該パラメタ
Figure 0007107246000101

を渡す。
記憶部４０は、データ記憶部４１と、パラメタ記憶部４２とを含む。データ記憶部４１には、打ち切りデータ
Figure 0007107246000102

が格納される。また、パラメタ記憶部４２には、モデルのパラメタ
Figure 0007107246000103

が格納される。
入力部５０は、外部装置２から、打ち切りデータ
Figure 0007107246000104

の入力を受け付ける。そして、入力部５０は、受け付けた打ち切りデータ
Figure 0007107246000105

を、データ処理部１０に渡す。
出力部６０は、パラメタ出力部３０から受け取ったモデルのパラメタ
Figure 0007107246000106

を、外部装置２へ出力する。
＜本発明の実施の形態に係る推定装置の作用＞
図７は、本発明の実施の形態に係る推定処理ルーチンを示すフローチャートである。
入力部５０に打ち切りデータ
Figure 0007107246000107

が入力されると、推定装置１において、図７に示す推定処理ルーチンが実行される。
まず、ステップＳ１００において、入力部５０は、観測値が観測されたサンプルの観測データ
Figure 0007107246000108

と、観測値が観測されなかったサンプルの観測データ
Figure 0007107246000109

と、各サンプルについて観測値が観測されたか否かを表す変数
Figure 0007107246000110

とを含む打ち切りデータ
Figure 0007107246000111

の入力を受け付ける。また、入力部５０は、各サンプルｉについての観測終了時刻までの時間の長さの入力を受け付ける。
ステップＳ１１０において、パラメタ推定部２０は、目的関数を、当該打ち切りデータ
Figure 0007107246000112

の確率分布を表すモデル
Figure 0007107246000113

と、当該打ち切りデータ
Figure 0007107246000114

の確率分布
Figure 0007107246000115

とのＫＬダイバージェンス、又はＬ２ダイバージェンスとして、上記式（１５）のパラメタ更新を繰り返すことにより、パラメタを推定するか、あるいは、上記式（２４）、（２５）を用いてパラメタを推定する。
ステップＳ１２０において、出力部６０は、上記ステップＳ１１０により推定されたパラメタ
Figure 0007107246000116

を出力する。
以上説明したように、本発明の実施形態に係る推定装置によれば、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、観測データの確率密度関数を用いて表される、当該打ち切りデータの確率分布を表すモデルと、当該打ち切りデータから得られる、当該打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、モデルのパラメタを推定するため、計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができる。
なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
例えば、上述の実施の形態では、ダイバージェンスとしてＫＬダイバージェンス又はＬ２ダイバージェンスを用いる場合について説明したが、これに限定されるものではなく、他のダイバージェンスを用いることも可能である。
また、上述の実施の形態では、時系列データである打ち切りデータを前提に記載しているが、これに限定されるものではなく、時系列データではない任意の打ち切りデータに対しても本発明の適用が可能である。
また、上述の実施の形態に係る推定装置１は、各部の処理をプログラムとして構築し、推定装置として利用されるコンピュータにインストールして実行させる構成として説明したが、ネットワークを介して流通させる構成としてもよい。
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
１ 推定装置
２ 外部装置
１０ データ処理部
２０ パラメタ推定部
３０ パラメタ出力部
４０ 記憶部
４１ データ記憶部
４２ パラメタ記憶部
５０ 入力部
６０ 出力部
１００ バス
１１０ ＣＰＵ
１２０ メモリ
１３０ 通信ＩＦ部
１４０ 入力部
１５０ 表示部
１６０ 記憶部
１７０ プログラム

Claims (5)

  1. 観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する推定装置であって、
    前記打ち切りデータの入力を受け付ける入力部と、
    前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定するパラメタ推定部と、
    を含む推定装置。
  2. 前記打ち切りデータの確率分布を表すモデルは、
    前記観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、前記変数の確率分布と、
    前記観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、前記変数が与えられた下での、前記観測データの確率分布と、
    を用いて表される請求項１記載の推定装置。
  3. 前記目的関数は、前記打ち切りデータの確率分布を表すモデルと、前記打ち切りデータの確率分布とのカルバックライブラーダイバージェンス、又はＬ２ダイバージェンスである
    請求項１又は２記載の推定装置。
  4. 観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する推定方法であって、
    入力部が、前記打ち切りデータの入力を受け付け、
    パラメタ推定部が、前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定する
    推定方法。
  5. 観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する処理をコンピュータに実行させるプログラムであって、
    入力部が、前記打ち切りデータの入力を受け付け、
    パラメタ推定部が、前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定する
    ことを含む処理をコンピュータに実行させるプログラム。
JP2019029769A 2019-02-21 2019-02-21 推定装置、推定方法、及びプログラム Active JP7107246B2 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2019029769A JP7107246B2 (ja) 2019-02-21 2019-02-21 推定装置、推定方法、及びプログラム
PCT/JP2020/004909 WO2020170867A1 (ja) 2019-02-21 2020-02-07 推定装置、推定方法、及びプログラム
US17/431,791 US20220138375A1 (en) 2019-02-21 2020-02-07 Estimation device, estimation method, and program

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2019029769A JP7107246B2 (ja) 2019-02-21 2019-02-21 推定装置、推定方法、及びプログラム

Publications (2)

Publication Number Publication Date
JP2020135554A JP2020135554A (ja) 2020-08-31
JP7107246B2 true JP7107246B2 (ja) 2022-07-27

Family

ID=72144974

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2019029769A Active JP7107246B2 (ja) 2019-02-21 2019-02-21 推定装置、推定方法、及びプログラム

Country Status (3)

Country Link
US (1) US20220138375A1 (ja)
JP (1) JP7107246B2 (ja)
WO (1) WO2020170867A1 (ja)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP7107246B2 (ja) * 2019-02-21 2022-07-27 日本電信電話株式会社 推定装置、推定方法、及びプログラム
US11921488B2 (en) * 2020-12-15 2024-03-05 Xerox Corporation System and method for machine-learning-enabled micro-object density distribution control with the aid of a digital computer
US20230359882A1 (en) * 2022-05-06 2023-11-09 International Business Machines Corporation Training a neural network to achieve average calibration
CN119128368B (zh) * 2024-11-14 2025-02-21 中国建设银行股份有限公司 回归模型的构建方法及装置、程序产品、存储介质

Family Cites Families (21)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040117051A1 (en) * 2000-06-06 2004-06-17 Ford Dean M. Method of determining a cumulative distribution function confidence bound
US8438126B2 (en) * 2007-07-13 2013-05-07 The Regents Of The University Of California Targeted maximum likelihood estimation
US8204714B2 (en) * 2010-02-26 2012-06-19 International Business Machines Corporation Method and computer program product for finding statistical bounds, corresponding parameter corners, and a probability density function of a performance target for a circuit
US20140149205A1 (en) * 2012-11-29 2014-05-29 Adobe Systems Incorporated Method and Apparatus for an Online Advertising Predictive Model with Censored Data
JP6194450B2 (ja) * 2013-04-15 2017-09-13 株式会社メガチップス 状態推定装置、プログラムおよび集積回路
JP6366999B2 (ja) * 2014-05-22 2018-08-01 株式会社メガチップス 状態推定装置、プログラムおよび集積回路
US10474770B2 (en) * 2014-08-27 2019-11-12 Nec Corporation Simulation device, simulation method, and memory medium
US20190095556A1 (en) * 2016-03-31 2019-03-28 Nec Corporation Information processing device, simulation method, and non-transitory recording medium storing simulation program
US20190057065A1 (en) * 2016-04-01 2019-02-21 Hitachi, Ltd. Modeling device
WO2018138880A1 (ja) * 2017-01-27 2018-08-02 三菱日立パワーシステムズ株式会社 モデルパラメータ値推定装置及び推定方法、プログラム、プログラムを記録した記録媒体、モデルパラメータ値推定システム
US10528644B1 (en) * 2017-06-30 2020-01-07 Cadence Design Systems, Inc. Estimation and visualization of a full probability distribution for circuit performance obtained with Monte Carlo simulations over scaled sigma sampling
JP6935765B2 (ja) * 2018-02-13 2021-09-15 日本電信電話株式会社 動的分布推定装置、方法、及びプログラム
JP7064356B2 (ja) * 2018-03-14 2022-05-10 株式会社日立製作所 将来状態推定装置および将来状態推定方法
WO2019217876A1 (en) * 2018-05-10 2019-11-14 Equifax Inc. Training or using sets of explainable machine-learning modeling algorithms for predicting timing of events
JP7172616B2 (ja) * 2019-01-11 2022-11-16 日本電信電話株式会社 データ解析装置、方法、及びプログラム
WO2020148940A1 (ja) * 2019-01-18 2020-07-23 株式会社ヒデ・ハウジング 日射量出現確率分布解析法、日射量出現確率分布解析システム、日射量出現確率分布解析プログラム、日射量正規化統計解析システム、日射量正規化統計解析法および日射量正規化統計解析プログラム
JP7107246B2 (ja) * 2019-02-21 2022-07-27 日本電信電話株式会社 推定装置、推定方法、及びプログラム
JP7215579B2 (ja) * 2019-06-26 2023-01-31 日本電信電話株式会社 パラメタ推定装置、パラメタ推定方法、及びパラメタ推定プログラム
US11010222B2 (en) * 2019-08-29 2021-05-18 Sap Se Failure mode specific analytics using parametric models
JP7268752B2 (ja) * 2019-10-02 2023-05-08 日本電信電話株式会社 パラメタ推定装置、パラメタ推定方法、及びパラメタ推定プログラム
JP7374868B2 (ja) * 2020-08-28 2023-11-07 株式会社東芝 情報処理装置、情報処理方法およびプログラム

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
BASU, Srabashi et al.，"Robust and efficient parametric estimation for censored survival data"，Annals of the Institute of Statistical Mathematics [online]，The Institute of Statistical Mathmatics，2006年，Vol. 58, No. 2，p.341-355，[2020年03月18日検索]，インターネット＜ＵＲＬ: https://www.ism.ac.jp/editsec/aism/pdf/058_2_0341.pdf＞，ISSN 0020-3157
幸島匡宏，ほか３名，"打ち切りデータに対する混合モデルのオンラインEM法の導出と大規模集客イベントにおける到着時間分布推定"，第10回データ工学と情報マネジメントに関するフォーラム（第16回日本データベース学会年次大会） ［Ｏｎｌｉｎｅ］，電子情報通信学会データ工学研究専門委員会 日本データベース学会 情報処理学会データベースシステム研究会，2018年03月，［2018年04月17日検索］，インターネット＜ＵＲＬ：http://db-event.jpn.org/deim2018/data/papers/272.pdf＞

Also Published As

Publication number Publication date
JP2020135554A (ja) 2020-08-31
US20220138375A1 (en) 2022-05-05
WO2020170867A1 (ja) 2020-08-27

Similar Documents

Publication Publication Date Title
JP7107246B2 (ja) 推定装置、推定方法、及びプログラム
Taylor et al. Performance estimation toolbox (PESTO): Automated worst-case analysis of first-order optimization methods
Henseler On the convergence of the partial least squares path modeling algorithm
Song et al. A Bayesian modeling approach for generalized semiparametric structural equation models
Nakayama Confidence intervals for quantiles using sectioning when applying variance-reduction techniques
Jenkinson et al. Numerical integration of the master equation in some models of stochastic epidemiology
Fried et al. Retrospective Bayesian outlier detection in INGARCH series
CN114492823A (zh) 消除量子噪声的方法及装置、电子设备和介质
US8813009B1 (en) Computing device mismatch variation contributions
Mohammadi et al. An Introduction to the BDgraph for Bayesian Graphical Models
Gåsemyr On an adaptive version of the Metropolis–Hastings algorithm with independent proposal distribution
Antoniano‐Villalobos et al. A nonparametric model for stationary time series
Song et al. Wasserstein generative regression
de la Fuente et al. An efficient nonlinear programming strategy for PCA models with incomplete data sets
Liu et al. On Random Batch Methods (RBM) for interacting particle systems driven by Lévy processes
Mudholkar et al. Transformation of the bathtub failure rate data in reliability for using Weibull-model analysis
Karimi et al. An approximate expectation maximization algorithm for estimating parameters, noise variances, and stochastic disturbance intensities in nonlinear dynamic models
EP3712784A2 (en) System and method for signal pre-processing based on data driven models and data dependent model transformation
Wang et al. Constrained spline regression in the presence of AR (p) errors
Yurko et al. Demonstration of Emulator‐Based Bayesian Calibration of Safety Analysis Codes: Theory and Formulation
Zhou Convergence estimates of nonrestarted and restarted block‐Lanczos methods
Chan et al. Model predictive control of Hammerstein systems with multivariable nonlinearities
Allal et al. Parameter estimation for first-order random coefficient autoregressive (RCA) models based on Kalman filter
Xiao et al. Unified uncertainty analysis using the maximum entropy approach and simulation
Wang et al. Gaussian variational approximation for Bayesian Lasso quantile regression model with zero-or-one inflated proportional data

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20210527

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20220614

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20220627

R150 Certificate of patent or registration of utility model

Ref document number: 7107246

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

S533 Written request for registration of change of name

Free format text: JAPANESE INTERMEDIATE CODE: R313533

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350