DEIM Forum 2019 A4-1

打ち切りデータを入力とする例に基づく混合モデル
幸島 匡宏†

松林 達史†

戸田 浩之†† 日本電信電話株式会社 NTT サービスエボリューション研究所〒 239-0847 神奈川県横須賀市光の丘 1-1
E-mail: †{kohjima.masahiro,matsubayashi.tatsushi,toda.hiroyuki}@lab.ntt.co.jp
あらまし 本研究では, 多次元の打ち切りデータを分析する新たな機械学習手法を提案する. 提案手法は例に基づくモデルのアプローチを利用することで, 大域的最適解への収束が保証でき, かつ数値計算を繰り返し実行する必要のない,
混合モデルの推定が可能となる手法である. 提案手法の有効性を検証した実験結果について報告する.
キーワード

打ち切りデータ, 混合モデル, 例に基づくモデル, KL ダイバージェンス, L2 ダイバージェンス, 最小二乗保証されるものの収束先は初期値に依存して変わるため, 異な

1 はじめに打ち切りデータとは, 値がある閾値以上のデータについては値が観測されず, 閾値以上であるという情報しか得ることができないデータのことを指す [1] [2]. 打ち切りデータの代表例は生存分析と呼ばれる領域で分析対象となる機器の故障時間デー
タである (図 1(a)). 観測期間内に故障せず観測が打ち切られた機器の故障までに要する時間 (生存時間) は, その時点までの稼
働時間以上としかわからないため, 打ち切りデータとして表現される. 近年の生存分析では E コマースサイト利用者のサービ
ス利用期間データ [3] や学生の就職活動の継続期間データ [4] なども分析対象になっている.
上記で挙げた例では, 基本的に 1 次元の打ち切りデータから,

1 次元の確率分布を推定する問題となる. しかしながら, 多次
元の打ち切りデータから, 多次元確率分布を推定することが必
要な問題も存在する. たとえば, 家族の構成員 (夫, 妻, 子など)の生存期間を推定する問題を考える. 子が死亡した場合に妻の死亡確率が増すという報告 [5] も存在することから, 精度の良い推定のためには構成員それぞれの生存期間を独立に扱うのは不十分であり, 構成員数次元の時刻上の同時確率分布を推定する必要があると考えられる. 動画, 書籍, 音楽など複数の定額配信サービスにおけるユーザ利用期間を分析する問題 (図 1(b)) も同様である. たとえば動画サービスがユーザにとって非常に魅力的なサービスとなっている場合, 動画サービスの解約時にはその他のサービスの解約も同時にまとめて起こる可能性が高いと考えられる. よって精度の良い推定のためにはサービス数次
元の同時確率分布を推定する必要がある.
確率モデルとして混合モデルを用いて打ち切りデータを扱う
場合, パラメタの推定には Expectation Maximization (EM) ア
ルゴリズム [6] を拡張した EMCM (Expectation Maximization

for Censored Mixture models) [7] や Variational Bayes (VB)
アルゴリズム [8] [9] を拡張した VBCM (Variational Bayes for

Censored Mixture models) [10] が利用できる. しかしながら,これらのアプローチは次の 2 つの問題がある. 1 つ目は局所最適解の存在である. EMCM も VBCM も目的関数の単調減少はる初期値からの繰り返し実行が必要となる. 2 つ目はモデルで利用する確率分布の切断分布 (ある一定の範囲でしか値をとらないように変更した分布) の統計量の計算の必要性である. 1 次
元正規分布の切断分布である 1 次元切断正規分布などの例外的なものを除き, この統計量は解析的に計算することができず, 何
らかの数値計算を利用することが必要となる. EMCM, VBCM
ともにパラメタの更新と統計量の計算を何度も繰り返すアルゴ
リズムであるから, 各反復で数値計算を繰り返すことは避けることが望ましいと考えられる.
そこで本研究では, 多次元の打ち切りデータを入力とし, 大域的最適解への収束が保証できかつ数値計算の繰り返し実行を必要としない混合モデル, 例に基づく打ち切り混合モデルを提案する. 上記を可能とするアイディアは, 手法名に含まれる例に基づくモデル (Examplar based model) のアプローチ [11] [12]の利用にある. このアプローチでは, 混合モデルの各コンポー
ネントのパラメタは陽にパラメタとしては扱わず, データ点の存在する点に各コンポーネントを配置する. これによって, 混
合比 (各コンポーネントの重みパラメタ) のみの推定を行う定式化をすることで大域的最適解へ収束するアルゴリズムが構築できる. 我々の提案する手法は, このアプローチを打ち切りデー
タを入力として利用できるよう発展させたものとみなすことができる. さらに前述した EMCM と VBCM で必要となる切断分布の統計量はコンポーネントのパラメタの推定に必要なものであったため, 数値計算の繰り返し実行は必要なくなる. よって多次元の打ち切りデータを分析するという問題において, 例に基づくモデルのアプローチは, 2 つの困難さを同時に解決する, 一石二鳥のアプローチと言える.
提案手法のパラメタ推定のため, 我々は次に示す目的関数の異なる 2 つの推定手法を構築した. 1 つ目は確率分布間のKullback-Leibler (KL) ダイバージェンス最小化に基づく推定手法, 2 つ目は確率分布間の L2 ダイバージェンス最小化に基づく推定手法である. 1 つ目の推定手法では非常にシンプルな繰
り返し計算でパラメタの推定が可能であり, 2 つ目の推定手法
では, 制約付きの 2 次の最適化問題を解くことでパラメタの推定が可能である.

(a) 1 次元打ち切りデータ

図1(b) 2 次元打ち切りデータ

打ち切りデータの例. (a) 機器の故障時刻を表す 1 次元の打ち切りデータを表す. 黒色の点が故障時刻が観測されたデータ, 赤色の点が観測打ち切りにより故障時刻が観測されなかったデータを表す. (b) 利用者のサービス利用時間を表す 2 次元の打ち切りデータを表す. 黒色の点が全サービスの解約時刻が観測されたデータ, 緑色/水色/赤色の点が, 観測打
ち切りによりそれぞれサービス 1/サービス 2/サービス 1 と 2 両方の解約時刻が観測されなかったデータを表す.

本稿の構成は次の通りである. § 2 で関連研究, § 3 で提案手

計算の繰り返しを避ける手法を構築した.

法構築の準備として打ち切りデータと例に基づくモデル, ダイ

例に基づくモデルに関する研究は文献 [11] [12] [21] などでこバージェンスについて記す. § 4 で KL ダイバージェンスを利れまで進められており, 局所最適解のない, 大域的最適解が求ま用する提案手法, § 5 で L2 ダイバージェンスを利用する提案手

るアルゴリズムが構築されている. 文献 [11] では KL ダイバー法を示す. § 6 でその有効性を実験的に検証し, § 7 でまとめる.

ジェンスを最小化することで推定する手法, 文献 [12] [21] では
確率分布間の二乗誤差を最小化することで推定する手法を提案

2 関連研究近年, 打ち切りデータを分析する機械学習のモデルやアルゴ
リズムが多数提案されている [10] [13] [14] [15]. この研究も上記と同様に例に基づく混合モデルのアプローチを用いた打ち切り
データの分析手法を提案する研究として位置付けられる.
打ち切りデータを扱うモデルでは, 利用する確率分布の切断
分布 (ある一定の範囲でしか値をとらないように変更した分布)の統計量や元となる分布の累積確率の計算が必要となることが多い. 古くから研究は進められてはいるものの [16] [17] [18] [19],これらの計算をする環境は十分整えられているとはいえない 1 .
また存在しても 2 , しばしば数値計算を内部で必要とする. たRでとえば MathWorks 社の数値計算ソフトウェア MATLAB⃝

している

4

. しかしながら上記の手法のいずれも打ち切りデー

タを入力とすることはできない. 我々の研究はこれら 2 つのアプローチを入力として打ち切りデータが利用できるように拡張したものとみなすことができる.
打ち切りデータを利用し, かつ 2 乗誤差を用いる研究には文献 [22] が存在する. しかしこの文献では, 指数分布, ワイブル分
布のように単峰の表現力の低い分布をモデルとして用いる場合が考えられている. 本研究は, 混合モデルのように表現力の高
い分布を用いた場合を考えている点が異なる.

3 準

備

3. 1 打ち切りデータ

は, 多次元正規分布の累積確率の計算に 3 次元以下の場合は適

まず打ち切りデータについて説明する. はじめに図 1(a) の打応求積法 [18], 4 次元以上の場合は準モンテカルロの積分アルゴ

ち切りデータの代表例である機器故障データを用いて説明する.

リズム [17] を利用している 3 . したがってこれらの計算を繰りこのデータでは, 各機器の設置された時刻と故障が起こった際返すことは避けることが望ましいと考えられる. 前述したようの故障の時刻が記録されている. 機器 1 と 2 に関しては観察期に我々は例に基づくモデルのアプローチを利用することで数値

間中に故障が起きているため, その時刻が記録されているが, 機

1：た と え ば python の 数 値 計 算 ラ イ ブ ラ リ scipy の 統 計 関 数 モ ジュー ル
scipy.stats で実装されている切断分布は 1 次元の切断正規分布と切断指数器 3 と 4 に関しては観察期間中には故障せずに, 観察が打ち切られたため, 故障の時刻が記録されてはいない. しかしながら,分布のみである.

2：たとえば R の多次元切断正規分布パッケージ “tmvtnorm” [20]3：https://jp.mathworks.com/help/stats/mvncdf.html

4：厳密にはこの手法は条件付確率分布を推定する手法として提案されている

(a) 1 次元打ち切りデータの生存時間表現

図2(b) 2 次元打ち切りデータの生存時間表現

打ち切りデータ (図 1) の生存時間表現. (a) 横軸を機器故障までの生存時間, (b) 横軸と縦軸が各サービスのユーザの利用時間を表す.

機器 3 と 4 も (いつかは必ず故障するとして) 故障の時刻は観ルであり, xi ∈ Rdx , wi ∈ {0, 1}dx である. xij が利用者 i の察終了時刻以降にあるということは読み取れる. このように機サービス j の利用時間, wij が利用者 i のサービス j の解約時

器 1 と 2 のように値 (故障の時刻) が分かるデータと機器 3 と刻が記録されたか (wij = 0), 打ち切りにより記録されなかった

4 のように値 (故障の時刻) がある値以上であること分かるデー

か (wij = 1) を表す. 同様に i 番目データの打ち切りの閾値を5タの組からなるデータを打ち切りデータと呼ぶ .

vi ∈ Rdx と書く. vij が利用者 i の サービス j 利用開始時刻か

上記の機器故障データは 1 次元の打ち切りデータであるが,ら観測終了時刻までの長さを表す. 打ち切りにより値が観測さ本研究では 2 次元以上の打ち切りデータも扱うためここで説明れなかったとき (wij = 1) には, xij = vij と設定されているとする. 図 1(b) はある 2 つのサービスのユーザの利用期間を表する.

す 2 次元の打ち切りデータである. 利用者の (少なくとも 1 つの) サービスを利用開始した時刻と各サービスを追加契約または解約した時刻が記録されている. 利用者 1 は両サービスを同時に利用開始し, 観察期間中に同時に解約しているため両方の

3. 2 例に基づく混合モデル
次に提案手法で用いるモデルについて説明する. 混合モデルの確率密度関数は一般に次の式で定義される.

サービスの解約時刻が記録されている. 利用者 2 は両サービス

f (x|θ) =を同時に利用開始し, 観察期間中にサービス 2 のみを解約して∑K
k=1

θk ψk (x) = θ T ψ(x).

(1)

いる. 利用者 3 は両サービスを同時に利用開始し, 観察期間中

ただし, K は混合数, ψk (x) が k 番目コンポーネントの確率分にサービス 1 のみを解約している. 利用者 4 はサービス 2 を布を表す. コンポーネントの確率分布には例えばガウス分布

まず利用開始し, 観察期間中にサービス 1 を追加契約している.
よって観測の打ち切りによって, 利用者 2 のサービス 1 の解約( ∥x − µ ∥2 )1kψk (x) = N (x|µk , σ 2 ) = √
exp −
2σ 2
( 2πσ 2 )dx

時刻, 利用者 3 のサービス 2 の解約時刻, 利用者 4 のサービス(2)

1 と 2 の解約時刻は記録されていない. このように 2 次元の打
ち切りデータにおいては, 打ち切られている値の次元が異なる,が利用できる. µk と σ はガウス分布の平均と標準偏差を表

3 通りの打ち切りの種類が存在する. 一般に n 次元の打ち切り

す. Examplar based model [11] のアプローチでは, コンポー

データでは 2 − 1 種類の打ち切りが存在する.

ネントの確率分布のパラメタ (ガウス分布の場合, µ1 , · · · , µK )n6

ここで打ち切りデータの定義を与える. 扱いを簡単にするたは K = n として, それぞれが観測データ点に対応するようめ, データは図 1(a) のようなカレンダー時刻ではなく, 生存時(µ1 , · · · , µn ) = (x1 , · · · , xn ) と設定されたものなどが用いら

間 (機器の設置から故障までの時間, サービスの契約から解約れる. 7 本研究では, このアプローチに従いかつ打ち切りデー

までの時間) を用いて表現する. 図 1 のデータを生存時間で表タを扱うため, K = n とし, 値が観測されていれば (wij = 0),

現したものを図 2 に示す. 多次元の打ち切りデータとして定義µij = xij , そうでなければ (wij = 1), µij = xij + ϵ と設定しするため以下では図 2(b) を説明の例に用いる. 打ち切りデー

た. ϵ は 0 以上の値をとる確率分布からランダムに生成した値

タを D={xi , wi }ni=1 と書く. xi , wi はともに dx 次元のベクトであり, 後の分布では指数分布を利用した. なお, σ は交差検証法などで決定することが可能である.

5：この例の場合, より正確には, 右側打ち切りされたデータと呼ばれる
6：本研究では各次元ごとに観測が打ち切られるかが決まる状況を考えるが, ど

7：事前知識に基づいて設定されたコンポーネントも利用されてもよいし, デー

れか一つの要素が打ち切られたときに全ての要素が観測されなくなるという状況

タ数が多い場合には, たとえばランダムに選んだ 100 個のデータのみを用いてもでも同様のアプローチで対応することができる.

よい.

上記モデルを用いたときの打ち切りデータ D = {xi , wi } の生成過程は次のように記述できる. まず打ち切りの閾値 vi が既知のもと, 打ち切りがおこるか否かを表す変数 wi が, 下記の確率分布に従い生成される.

P (wi |θ) = F (wi |θ; vi ) ==∫ ∞ {∫ vio
vic[∫ ∞ {∫ vic

∑K
k=1

θk
vic

−∞

}
f (xi |θ)dxoi dxci}]
ψk (x)dxo dxc .

(3)

−∞

ただし, xi と vi のうち値が観測された要素の集合を xoi ={xij |wij = 0},vio = {vij |wij = 0}, それ以外の値が観測されなかった集合を xci = {xij |wij = 1}, vic = {vij |wij = 1} と書いた. また, コンポーネントの分布 ψk (x) に式 (2) のガウス分布のように標準的な確率分布を利用すれば, 式中の積分は累積密

図3

正規分布 N (x|µp , σ 2 ) と N (x|µq , σ 2 ) の KL ダイバージェンスと L2 ダイバージェンス. 横軸が µp , 縦軸が µq = 0, σ = 1.0と固定した際の各ダイバージェンスの値を表す.

度関数を用いて計算できる.

wi =| 1 であり, 少なくとも一つ観測された要素がある場合,
xi は以下の分布にしたがって生成される.

3. 3 KL ダイバージェンスと L2 ダイバージェンス
ここで次の章以降に提案手法の目的関数を定義する際に利P (xi |wi =| 1, θ) = δ(xci − vic )ftr (xoi |wi , θ; vi ).

(4)用するダイバージェンスについて記す. 良く知られるように確率分布 p(x) と q(x) に対するカルバックライブラー (KL) ダイ

δ(·) はデルタ関数であり, ftr は f を観測されなかった要素に関して周辺化した分布の切断分布である.
 o
 g(xi |wi ,θ;vi )F (wi |θ;vi )ftr (xoi |wi , θ; vi ) =0(if xoi <= vi )(otherwise)

バージェンスは以下の式で定義される.
[
q(x) ]KL(p, q) = Eq log
(10)
p(x)
∫
∫
= q(x) log q(x)dx − q(x) log p(x)dx.
これに加えて, 本研究では下記で定義される L2 ダイバージェン

ただし,g(xoi |wi , θ; vi ) =
=

∫ ∞
vic

f (xi |θ)dxci
(5)

∑K

T

Ψi (xoi , wi , vi ),θk Ψik = θ∫ ∞Ψik = Ψk (xoi , wi , vi ) =ψk (xi )dxci .
k=1

ス [23] も利用する.
∫
L2 (p, q) = {p(x) − q(x)}2 dx
(11)
∫
∫
∫
= {p(x)}2 dx + {q(x)}2 dx − 2 p(x)q(x)dx.

(6)

vic

L2 ダイバージェンスは 2 つの確率密度関数の 2 乗誤差としてコンポーネントが式 (2) のガウス分布の場合, Ψik は下記で与定義されている.
上記の違いを明確にするため p(x) と q(x) がそれぞれ平均

えられる.
xo

Ψik = N (xoi |µki , σ 2 )
xo

∫ ∞
vic

µp , µq , 分散 σ 2 の正規分布であるときの KL ダイバージェンスxc

N (xci |µki , σ 2 )dxci .

と L2 ダイバージェンスを以下に記す. KL ダイバージェンスは

xc

ただし, µki , µki はそれぞれ µk から xoi と xci に対応する次元の要素を抜き出したベクトルである. また, 式 (3)(5) の定義より, wi = 1 のときは,F (wi = 1|θ; vi ) = g(xoi |wi = 1, θ; vi )(7)が成立することを後ほど利用する.

wi = 1 であり, 一つも観測された要素がない場合はデルタ関数のみの表現となる.

P (xi |wi = 1, θ) = δ(xci − vi )で与えられる.

さらに L2 ダイバージェンスは()
L2 N (x|µp , σ 2 ), N (x|µq , σ 2 ){( ∥µ − µ ∥2 )}
2
p
q
= √
1 − exp −4σ 2(2 πσ 2 )dx
で与えられる.

µq = 0 と固定したときの各ダイバージェンスの値を図 3 に(8)

示す. KL ダイバージェンスは平均パラメタ µp が離れるに従い2 乗で大きくなるのに対し, L2 ダイバージェンスは平均パラメ

よって上記をまとめると, 各打ち切りデータの生成確率はP (xi , wi |θ) = P (xi |wi , θ)P (wi |θ)()∥µp − µq ∥2
KL N (x|µp , σ 2 ), N (x|µq , σ 2 ) =
,
2σ 2

(9)

タが十分離れると一定値となる.
どちらのダイバージェンスを用いるべきかは問題に応じて異なる. よって本研究では, いずれのダイバージェンスを利用する
場合にも適用できるように 2 種類の手法を構築することとした.

Algorithm 1 KL ダイバージェンスを利用した打ち切りデー

4 提案 1: KL ダイバージェンス最小化による推定
この章では, 目的関数として KL ダイバージェンスを利用する場合の提案手法を示す.

タを用いた例に基づく混合モデルの学習
Input: D: 入力データ, α: 閾値Output: θ̂: KL ダイバージェンス式 (14) を最小化するパラメタ
1: repeat2:

4. 1 目 的 関 数§ 3 で示したモデル P (x, w|θ) と打ち切りデータを生成する
真の確率分布 P ∗ (x, w) との KL ダイバージェンスは, 式 (10)3:4:for k = 1 to K do
θk を式 (16) に従い更新.
end for

5: until 収束条件が満たされるまで繰り返しの定義に従い,[P ∗ (xi , wi ) ]DKL (θ) =EP ∗ (xi ,wi ) log
P (xi , wi |θ)(12)

す. なお文献 [11] と同様に, 計算量の削減と収束を早めるため,で与えられる. これは以下のように式変形することができる.

DKL (θ)[∑ ∗P ∗ (xi |wi )P ∗ (wi ) ]=P (wi )EP ∗ (xi |wi ) log
P (xi |wi , θ)P (wi |θ)
wi =1|

P ∗ (wi = 1)
+ P (wi = 1) log
P (wi = 1|θ)[∑ ∗P ∗ (xoi |wi ) ]=P (wi )EP ∗ (xoi |wi ) log
g(xoi |wi , θ; vi )が可能となる. アルゴリズムの疑似コードを Algorithm 1 に示
パラメタ更新の際に θk がある閾値 α (たとえば 10−3 /n) より小さい場合には θk = 0 と設定した後に全体を和が 1 になるよ
うに調整する, 再正規化操作を行ってもよい.

5 提案 2: L2 ダイバージェンス最小化による推定

∗

wi =1
|

P ∗ (wi = 1)+C
+ P (wi = 1) log
P (wi = 1|θ)
[
]
∑ ∗
=
P (wi )EP ∗ (xo |wi ) − log g(xoi |wi , θ; vi )この章では, 目的関数として L2 ダイバージェンスを利用する場合の提案手法を示す.

(13)

∗

5. 1 目 的 関 数
ここでは, § 4 とは異なり式 (11) の定義から直接目的関数を定義することはせず, § 4 の KL ダイバージェンスを用いた際の目的関数に注目することで新たな目的関数を定義する. 8

wi =1
|

KL ダイバージェンスを用いた際の目的関数 (式 (13)) に注目− P ∗ (wi = 1) log F (wi = 1|θ; vi ) + C ′ .

すると, モデルの周辺分布 g と真の分布 P ∗ の KL ダイバージェ

定数項を除去し, 未知である真の分布 P ∗ に関する期待値を標ンスに対応する項と, モデル P (wi = 1) と真の分布 P ∗ (wi = 1)

本平均で置き換えれば, 下記の式が導かれる.

の対数尤度比に対応する項, これら 2 つの項の P ∗ (wi ) を用い

1 {∑
log g(xoi |wi , θ; vi ){i|wi =1}|n}∑+log F (wi = 1|θ; vi )(14)
{i|wi =1}
{
}
∑
(
)
n
1
log θ T Ψi (xoi , wi , vi ) .
=−
i=1
n

D̂KL (θ) = −最終行の式変形には式 (7) の関係を用いた. D̂KL (θ) はデータから計算できる量であり, これを目的関数とすることでアルゴ

た重み付き和で構成されていることが分かる.
この洞察に基づき, これら 2 つの各項で KL ダイバージェン
ス/対数尤度比を用いている箇所を L2 ダイバージェンスに置き換えることで, 目的関数を次のように設計することができる.
∫ (
)2
∑ ∗
DL2 (θ) =P (wi )
g(xoi |wi , θ; vi ) − P ∗ (xoi |wi ) dxoi
wi =1
|

(
)2
+ P ∗ (wi = 1) F (wi = 1|θ; vi ) − P ∗ (wi = 1) .

リズムを導く.
これを式変形すれば,

4. 2 アルゴリズム
前節で導出した目的関数を利用し, パラメタを推定する. 具∑P ∗ (wi )∫{}2g(xoi |wi , θ; vi ) dxoi

wi =1
|

体的には, 下記の最適化問題を解けばよい.

θ̂ = arg min D̂KL (θ) s.t. θk >= 0, |θ| = 1.0.

DL2 (θ) =

−2
(15)

θ

制約条件 (パラメタの要素が 0 以上でかつ和が 1) は, モデル f∑ ∫P ∗ (xoi , wi )g(xoi |wi , θ; vi )dxoi

wi =1
|

{
}2
+ P ∗ (wi = 1) F (wi = 1|θ; vi )が確率分布となるためのものである.
ラグランジュの未定乗数法を用いると, 上記の最適化問題の8：式 (11) から直接目的関数を下記のように定義することには 2 つ問題がある.
}2∑∫ { ∗
P (xi , wi ) − P (xi , wi |θ) dxi
D L2 (θ) =解は以下の式を満たすことが分かる.

1 {∑ θk Ψk (xoi , wi , vi ) }n i=1 θ T Ψi (xoi , wi , vi )
n

θ̂k =

wi

(16)

よって, 式 (16) に基づいて θ の更新を繰り返すことで最適化

1 つ目は, 値が定義できないデルタ関数の 2 乗の項が現れる. 2 つ目は, w の値
に関する総和をとる項があるため, 次元 dx のとき 2dx 通りの和をとることが必要となり, 10 次元程度でも計算が非常に困難になる.

{}2− 2 P ∗ (wi = 1) F (wi = 1|θ; vi ) + Const.

Algorithm 2 L2 ダイバージェンスを利用した打ち切りデータ

定数項を除去し, 真の分布 P ∗ に関する平均を標本平均で置き換えた, 下記の式の最適化を考える.
∫
{
}2
1{ ∑
D̂L2 (θ) =g(xoi |wi , θ; vi ) dxoinを用いた例に基づく混合モデルの学習
Input: D: 入力データ, β: ハイパーパラメタ
Output: θ̂: パラメタ
1: 式 (19) の最適化問題を解き θ̂ を得る{i|wi =1}|∑

−2

人工データ: 真の確率分布として用いる混合正規分布から

g(xoi |wi , θ; vi )生成することで作成した. 真の確率分布は, 混合数 K ∗ = 4,

{i|wi =1}
|

∑

+
−2
ただし, n1 =µ4 = (2, 2), 標準偏差は σ ∗ = 1.0 とした. 打ち切りの閾値は全}F (wi |θ; vi ) .

∑

n1
n

平均パラメタ µ1 = (−2, −2), µ2 = (2, −2), µ3 = (−2, 2),

F (wi |θ; vi )2{i|wi =1}

データ共通で vi = (1.5, 1.5) とした. 真の分布と生成した学習
データの例を図 4(a), 4(b) に示す. 評価用のテストデータは推{i|wi =1}∑ni=1 I(wi = 1) とおいた.

これを目的関数とす

定結果と真の確率分布との近さを測るため, 打ち切りはなく上記の閾値を超えた値も値が観測されるデータとなるよう作成しることでアルゴリズムを導出する.

た. 学習データとテストデータの数はそれぞれ 100, 1000 とし

5. 2 アルゴリズム

た. 上記の学習データとテストデータの生成を繰り返し, 計 5

D̂L2 (θ) は行列・ベクトル形式で次のように表すことができる.

実データ: Old Faithful Geyser Data 9 を利用した. この

D̂L2 (θ) = θ Ĥθ − 2ĥ θTつのデータセットを作成した.

T

(17)

データは米国イエローストーン国立公園にある間欠泉の噴出までの待ち時間と持続時間を表したものである. これは打ち切り

ただし,のないデータであるが, 待ち時間と持続時間がそれぞれの 3/4

n ∫1∑
Ψi (xoi , wi , vi )Ψi (xoi , wi , vi )T dxoi
n i=1
n1 ∑
1 ∑
Ψi + 2
ĥ =
Ψi .
n
n

分位点を閾値として, これ以上の値の場合には値が観測されなĤ ={i|wi =1}|いものとした. また x 軸と y 軸のスケールを合わせるために

y 軸の値は 0.1 倍した. 元データと作成したデータを図 5(a),{i|wi =1}

5(b) に示す. データを 5 分割し, 学習データとテストデータかコンポーネントの分布にガウス分布を用いる場合, Ĥ 中の積分
項は (累積密度関数を用いて) 計算できる値である.
(∫

∫ ∞
vic

xc

N (xci |µki , σ 2 )dxci

比較手法: 比較手法には, 提案手法の基となった, 打ち切り
データは本来扱えない “通常データ”を用いる Examplar Based)Ψi (xoi , wi , vi )Ψi (xoi , wi , vi )T dxo
i
kk′
∫
=
Ψk (xoi , wi , vi )Ψk′ (xoi , wi , vi )T dxo
i)( xo xo= N µk′i |µki , 2σ 2らなる計 5 つのデータセットを作成した.

Model の手法 [11] を利用した. この手法を利用するために打ち切りデータを下記の２つの方針で (見かけ上) 打ち切りのない∫ ∞vic

データに変換をした. (i) Baseline 1: 少なくとも 1 つ打ち切りにxc

N (xci |µk′i , σ 2 )dxciが観測されているデータのみを用いて推定する. (ii) Baseline 2:

これより, 目的関数は θ に関する 2 次の形式で表現されることが分かる. よってパラメタの推定値は下記で示す制約付きの 2
次最適化問題を解くことで得ることができる.

θ̂ = arg min D̂L2 (θ) s.t. θk >= 0, |θ| = 1.0.

より値が観測されていない要素があるデータは扱わず, 全ての値打ち切りにより値が観測されない要素は, 打ち切りの閾値を観測値として扱うことで打ち切りのないデータを作成し, これを用いて推定する. なお, 事前実験により人工データ実験, 実デー

(18)

θ

タ実験のそれぞれで各手法の標準偏差は共通に σ = 1.0, 0.25 と設定した. 提案手法 Proposed (SE) のハイパーパラメタは両実

これは標準的な数値ソルバーを用いて解くことができる. 以後の実験では上記の目的関数に正則化項を加えた下記の問題を解

験で β = 5.0 とした.
評価指標: 評価指標には下記の負の対数尤度 G を利用した.

くことで推定を行った.

[

]

G=

θ̂ = arg min D̂L2 (θ) + βθ θ s.t. θk >= 0, |θ| = 1.0. (19)Tθ

ただし, β はハイパーパラメタである.

6 実

験

1

n∑
test

ntest i=1

− log P (xtest
).
i

test
ただし, テストデータを {xtest}nii=1 と書いた. ntest はテスト

データ数を表す. 値が小さいほどモデルが正しく真の確率分布を表現できていることを表す.

6. 1 実 験 結 果
この章では提案手法の有効性を実験的に検証する. 値が観測されなかったデータを利用することで推定性能が向上することが確かめる実験を行った.

定量評価: 人工データと実データを用いた実験の定量評価結9：https://www.stat.cmu.edu/~larry/all-of-statistics/

表 1 実験結果: 負の対数尤度の平均と標準偏差 (小さいほど良い)
人工データ

実データ

Proposed (KL) 4.22 ± 0.05 2.13 ± 0.08
Proposed (L2)

4.27 ± 0.02 2.26 ± 0.15

Baseline 1

5.51 ± 0.44 3.08 ± 0.23

Baseline 2

4.33 ± 0.01 2.32 ± 0.16

[7]

[8]

[9]

果を表 1 に示す. 2 つの提案手法 Proposed (KL, L2) は比較手法 Baseline1, 2 よりも優れた性能を示していることが分かる.
これは提案手法は打ち切りデータを考慮した目的関数の設計が[10]されているためであると考えられる.
定性評価: 図 4 に人工データを用いた実験における推定結果[11]の確率密度関数を示す. 提案手法の推定結果 図 4(c), 4(d) は真の確率分布 図 4(a) に近いものとなっており, 真の確率分布[12]のピークのある位置を正確に捉えることができている. それに対し, Baseline 1 の推定結果 図 4(e) は, 全ての値が観測されたデータしか用いていないため図中の左下のピークの位置しか捉[13]

えられていない. また, Baseline 2 図 4(f) は, 4 つのピークがある構造は捉えられているものの, 値が観測されなかったデータ[14]については閾値を観測値としてしまっているため, 図中の左上,
右下, 右上のピークの位置が真の位置よりも原点に寄った結果となってしまっている. 提案手法は打ち切りデータを適切に扱[15]うことで, Baseline 1, 2 よりも真の確率分布の構造を捉えることができ, それが性能向上に貢献していると考えられる. 人工[16]

データを用いた実験結果を表す図 5 においても同様の結果が得られていることが分かる.

[17]

7 ま と め本研究では, 多次元の打ち切りデータを分析する, 例に基づく[18]モデルのアプローチを用いた手法を提案した. 提案手法は, 大域的最適解への収束が保証でき, かつ数値計算を繰り返し実行する必要なく混合モデルの推定が可能となる手法である. 人工

[19]
[20]

データを実データを用いた実験により提案手法の有効性を確かめた. 今後の発展の方向性としては, オンラインアルゴリズムの構築 [24] [25] などが挙げられる.

[21]

[22]

文

献

[1] David G Kleinbaum and Mitchel Klein. Survival analysis,Vol. 3. Springer, 2010.
[2] サバイバルデータの解析: 生存時間とイベントヒストリデータ.
バイオ統計シリーズ. 近代科学社, 2010.
[3] Shouichi Nagano, Yusuke Ichikawa, Noriko Takaya, Tadasu
Uchiyama, and Makoto Abe. Nonparametric hierarchal
bayesian modeling in non-contractual heterogeneous survival data. In KDD, pp. 668–676, 2013.
[4] 早川真央, 三川健太, 荻原大陸, 後藤正幸. 層別木と混合ワイブル
分布に基づく就職活動終了時期の分析モデルの構築. 情報処理学
会論文誌, Vol. 58, No. 5, pp. 1189–1206, may 2017.
[5] Javier Espinosa and William N Evans. Maternal bereavement: The heightened mortality of mothers after the deathof a child. Economics & Human Biology, Vol. 11, No. 3, pp.
371–381, 2013.
[6] Arthur P Dempster, Nan M Laird, and Donald B Rubin.
Maximum likelihood from incomplete data via the em algo-

[23]

[24]

[25]

rithm. Journal of the royal statistical society. Series B, pp.
1–38, 1977.
Didier Chauveau. A stochastic em algorithm for mixtureswith censored data. Journal of statistical planning and inference, Vol. 46, No. 1, pp. 1–25, 1995.
Hagai Attias. Inferring parameters and structure of latentvariable models by variational bayes. In UAI, pp. 21–30,
1999.
Michael I Jordan, Zoubin Ghahramani, Tommi S Jaakkola,and Lawrence K Saul. An introduction to variational methods for graphical models. Machine learning, Vol. 37, No. 2,
pp. 183–233, 1999.
Masahiro Kohjima, Tatsushi Matsubayashi, and Hiroyuki
Toda. Variational bayes for mixture models with censoreddata. In ECMLPKDD, pp. 605–620, 2018.
Danial Lashkari and Polina Golland. Convex clustering with
exemplar-based models. In NIPS, pp. 825–832, 2008.
Masashi Sugiyama, Ichiro Takeuchi, Taiji Suzuki, Takafumi
Kanamori, Hirotaka Hachiya, and Daisuke Okanohara. Conditional density estimation via least-squares density ratio
estimation. In AISTATS, pp. 781–788, 2010.
Sebastian Pölsterl, Nassir Navab, and Amin Katouzian.
Fast training of support vector machines for survival analysis. In ECMLPKDD, pp. 243–259, 2015.
Farkhondeh Kiaee, Hamid Sheikhzadeh, and Samaneh Eftekhari
Mahabadi. Relevance vector machine for survival analysis.
IEEE transactions on neural networks and learning systems, Vol. 27, No. 3, pp. 648–660, 2016.
Tamara Fernández, Nicolás Rivera, and Yee Whye Teh.
Gaussian processes for survival analysis. In NIPS, pp. 5021–
5029, 2016.
Georges M Tallis. The moment generating function of thetruncated multi-normal distribution. Journal of the Royal
Statistical Society. Series B, pp. 223–229, 1961.
Alan Genz and Frank Bretz. Numerical computation ofmultivariate t-probabilities with application to power calculation of multiple contrasts. Journal of Statistical Computation and Simulation, Vol. 63, No. 4, pp. 103–117, 1999.
Alan Genz. Numerical computation of rectangular bivariate and trivariate normal and t probabilities. Statistics andComputing, Vol. 14, No. 3, pp. 251–260, 2004.
Manjunath BG and Stefan Wilhelm. Moments calculation
for the double truncated multivariate normal density. 2009.
Stefan Wilhelm and BG Manjunath. tmvtnorm: A package
for the truncated multivariate normal distribution. sigma,
Vol. 2, No. 2, 2010.
Masashi Sugiyama, Taiji Suzuki, and Takafumi Kanamori.
Density ratio estimation in machine learning. Cambridge
University Press, 2012.
Srabashi Basu, Ayanendranath Basu, and MC Jones. Robust and eﬃcient parametric estimation for censored survival data. Annals of the Institute of Statistical Mathematics, Vol. 58, No. 2, pp. 341–355, 2006.
David W Scott. Parametric statistical modeling by minimum integrated square error. Technometrics, Vol. 43, No. 3,
pp. 274–285, 2001.
Olivier Cappé and Eric Moulines. On-line expectation–
maximization algorithm for latent data models. Journalof the Royal Statistical Society: Series B, Vol. 71, No. 3,
pp. 593–613, 2009.
幸島匡宏, 清武寛, 松林達史, 塩原寿子, 戸田浩之. 打ち切りデー
タを用いた混合モデルのオンライン em アルゴリズム. 電子情報通信学会論文誌, Vol. J102-D, No. 4, 2019. to appear.

(a) 真の確率分布(b) データ(c) 提案手法 (KL)(d) 提案手法 (L2)(e) Baseline 1(f) Baseline 2

図 4

人工データを用いた実験における (a) 真の確率分布と (b) 学習に用いたデータ. 白抜きの点が 2 次元全ての値が観測されたデータ, 青/緑の点が x 軸/y 軸の値のみが観測されたデータ, 赤点が両方の値が打ち切りにより観測されなかったデータを表す. (c) KL ダイ
バージェンスを用いた提案手法, (d) L2 ダイバージェンスを用いた提案手法, (e)(f) 各ベー
スライン手法の推定結果.

(a) 元データ(b) データ(c) 提案手法 (KL)(d) 提案手法 (L2)(e) Baseline 1(f) Baseline 2

図5

実データを用いた実験における (a) 元データと (b) 学習に用いるデータ. 白抜きの点が 2次元全ての値が観測されたデータ, 青/緑の点が x 軸/y 軸の値のみが観測されたデータ, 赤
点が両方の値が打ち切りにより観測されなかったデータを表す. (c) KL ダイバージェンスを用いた提案手法, (d) L2 ダイバージェンスを用いた提案手法, (e)(f) 各ベースライン手法の推定結果.