全文 - terrapub

日本統計学会誌
第 44 巻, 第 1 号, 2014 年 9 月
97 頁 ∼ 111 頁
特集
密度比によるダイバージェンス推定とその応用
金森敬文∗
Divergence Estimation Using Density-Ratio and its Applications
Takafumi Kanamori∗
密度比は 2 つの確率密度の比として定義され，さまざまな統計的推論と関連し，応用されてい
る．本稿では，まず密度比の推定法について紹介する．また密度比を用いて，2 つの確率密度の
非類似度を測る尺度であるダイバージェンスを推定する方法について解説する．さらにダイバー
ジェンス推定を 2 標本検定に応用した結果について述べる．理論的考察や数値実験を通して，提
案法の優位性を示す．
A density ratio is deﬁned by the ratio of two probability densities, and it is widely applied
to real world statistical problems. In this paper, we introduce moment matching estimators
of density ratios and f -divergence estimators based on density-ratios. Then, we exploit the
f -divergence estimators to the two-sample homogeneity test. We prove that the proposed test
dominates the existing empirical likelihood score test. Through numerical studies, we illustrate
the adequacy of the asymptotic theory for ﬁnite-sample inference.
キーワード: 密度比，ダイバージェンス，2 標本検定，検出力
はじめに
1.
2 つの確率密度関数の比を密度比という．密度比は数学的にはラドン・ニコディム微分
と等価であり，また統計学においても，症例対照研究などで古くからその推定法が議論さ
れている (Qin (1998))．また近年，おもに機械学習やデータマイニングと呼ばれる分野に
おいて，密度比がさまざまな統計的推論の問題と関連することが明らかにされ活発に研究
が進められている (Sugiyama et al. (2012))．
本稿では，まず密度比について簡単に紹介し，さらにその推定方法を示す．密度比の応
用として，2 つの確率密度の非類似度を測る尺度であるダイバージェンスの推定問題を考
える．ダイバージェンスは統計科学や情報理論など数理科学全般において重要な概念であ
る (Kullback and Leibler (1951), Cover and Thomas (2006))．統計学においては外れ値検
∗
名古屋大学情報科学研究科：〒 464-8601 愛知県名古屋市千種区不老町 (E-mail: [email protected]).
98
日本統計学会誌第44巻第1号 2014
出やデータの次元削減などに広く応用され，また情報理論では符号長の冗長度として解釈
されている．このため，ダイバージェンスの推定法を確立することは重要な課題となって
いる．本稿では密度比を用いたダイバージェンスの推定方法を提案する．また 2 標本検定
の問題にダイバージェンス推定を応用し，理論的な結果を示す．数学的条件の詳細につい
ては Kanamori et al. (2012) に記されている．
密度比
2.
確率密度関数 pd (x), pn (x) に対して，密度比 r(x) は
r(x) =
pn (x)
pd (x)
(2.1)
と定義される．密度比の分子 (numerator) を pn ，分母 (denominator) を pd としている．密
度関数の定義域上で pd > 0 を仮定し，r(x) の値が常に定義されるとする．以下のように，
それぞれの確率分布から標本が独立に得られている状況を考える．
(n)
x1 , . . . , x(n)
mn ∼i .i .d . pn ,
(d)
x1 , . . . , x(d)
md ∼i .i .d . pd .
(2.2)
標本から密度比を推定する方法と応用例について解説する．
本稿で用いる記号を以下に示す．実数の集合を R, 非負実数の集合を R+ とする．確率
密度 pd (pn ) の下での期待値を Ed [ · ] (En [ · ])，分散を Vd [ · ] (Vn [ · ]) とする．また標本 (2.2)
の下での期待値を E[ · ]，分散を V[ · ]，共分散を Cov[ ·, · ] とする．
2.1
密度比の推定法
密度比に対する統計モデル
r(x; θ),
θ ∈ Θ ⊂ Rd
を設定する．ここでパラメータ空間 Θ はユークリッド空間 Rd の開集合とする．密度比モ
デル r(x; θ) は密度比に対してはパラメトリックモデルであるが，確率分布に対してはセミ
パラメトリックモデルである．すなわち，密度比 r(x) に対して r(x) = pn (x)/pd (x) を満た
す密度関数の組 pd , pn は一意には定まらず，一般に無限次元の自由度がある．このため密
度比モデルは，確率分布に対して柔軟なモデリングとなっていると解釈できる．ここでは
簡単のため，密度比 r(x) はモデルに含まれており，r(x) = r(x; θ∗ ), θ∗ ∈ Θ と表せると仮
定する．この仮定が正確には成り立たない場合の統計的推論についても，4 節で考察する．
標本から密度比 r(x) を推定するためにモーメント法を用いることができる (Qin (1998))．
関数 η(x) に対して，密度比の定義から次の恒等式
∫
∫
η(x)r(x)pd (x)dx − η(x)pn (x)dx = 0
密度比によるダイバージェンス推定とその応用
99
が成り立つ．関数 η(x) がベクトル値関数であっても同様である．上式を標本平均で近似す
ることで，密度比推定のためのモーメント法を得ることができる．一般にパラメータ θ に
依存する関数 η(x; θ) に対して，推定関数 Qη (θ) を
md
mn
1 ∑
1 ∑
(d)
(d)
(n)
η(xi ; θ)r(xi ; θ) −
Qη (θ) =
η(xj ; θ)
md i=1
mn j=1
とおく．大数の法則より Qη (θ) は md , mn → ∞ のとき Ed [η(x; θ)r(x; θ)] − En [η(x; θ)] に
確率収束する．この式は θ = θ∗ で 0 になる．したがって近似的に Qη (θ∗ ) = 0 が成り立つ
と考えれる．以上より，パラメータ θ∗ に対する推定方程式として
Qη (θ) = 0,
θ∈Θ
(2.3)
を考える．ここでモーメント関数 η(x; θ) は Rd に値を取るベクトル値関数としておくと，
d 次元パラメータ θ に対して d 本の方程式が得られる．正則条件の下で解が定まり，方程
式 (2.3) から θ∗ を推定することができる．実際，この推定法はフィッシャー一致性を持ち，
さらに適当な正則条件の下で漸近一致性を持つことが分かる (Qin (1998), Kanamori et al.
(2012))．本稿では，標本数 md , mn に対して md , mn → ∞ のとき md /mn −→ ρ ∈ (0, ∞)
が成り立つことを仮定する．
推定方程式 (2.3) による推定量の精度については Qin (Qin (1998)) によって詳細に調べ
られている．以下，その結果を簡単に紹介する．密度比モデル r(x; θ) は正の定数倍に関し
て閉じており，
r(x; θ) = eθ1 r1 (x; θ¯r ), θ = (θ1 , θ2 , . . . , θd ) ∈ Rd , θ¯r = (θ2 , . . . , θd ) ∈ ×Rd−1
と表せるとする．推定方程式 (2.3) をパラメータ θ∗ のまわりで漸近展開する．密度比モデ
ルのパラメータ θ に関するベクトル微分演算を ∇ で表し，行列 Uη (θ) を
T
Uη (θ) = En [η(x; θ)∇ log r(x; θ) ] ∈ Rd×d
(2.4)
とする．また m を md と mn の調和平均 1/(1/md + 1/mn ) とする．行列 Uη (θ) が可逆な
ら，方程式 (2.3) の解 θb に対して
√
(
)
√
d
∗
−1
−1 ρVd [rη] + Vn [η]
T −1
b
m(θ − θ ) = − mUη Qη + op (1) −→ Nd 0, Uη
(Uη )
ρ+1
が成り立つ．ここで関数値は θ = θ∗ での評価値とする．漸近分散の大きさは，モーメント
法に用いる関数 η(x; θ) に依存する．漸近分散を一様に最小にする関数 η(x; θ) は
ηopt (x; θ) =
1
∇ log r(x; θ)
1 + ρr(x; θ)
(2.5)
100
日本統計学会誌第44巻第1号 2014
となることが，Qin (Qin (1998)) により示されている．実際，簡単な計算から任意のモー
メント関数 η に対して，
m · Cov[(Uη−1 Qη − Uη−1
Qηopt ), Uη−1
Qηopt ] = o(1)
opt
opt
となることが分かり，漸近的に
0 ≤ V [Uη Qη − Uη−1
Qηopt ]
opt
= V [Uη−1 Qη ] − V [Uη−1
Qηopt ] + Cov[(Uη−1 Qη − Uη−1
Qηopt ), Uη−1
Qηopt ]
opt
opt
opt
= V [Uη−1 Qη ] − V [Uη−1
Qηopt ]
opt
が成り立つ．なお d ≥ 2 のときには，上の不等式は行列に対する非負定値性を意味する．
関数 ηopt (x; θ) はロジスティック回帰の対数尤度関数に一致する．
2.2
密度比の応用
密度比を用いる学習の例として，共変量シフトの下での回帰分析を挙げる．詳細は
Shimodaira (2000), Sugiyama and Kawanabe (2012) に詳しい．教師つき学習データ
(x1 , y1 ), . . . , (xn , yn ) から回帰関数を推定する問題を考える．学習データの同時確率密度関数
を p(y|x)pd (x) とする．条件付き確率密度関数 p(y|x) の条件付き期待値 f (x) = E[Y |X = x]
が二乗損失ものとで最適な回帰関数を与える．回帰関数に対する統計モデルとして，パラ
メータ α をもつ f (x; α) を仮定する．通常の最小 2 乗推定量 f (x; α
b) は次の最適化問題を解
くことで得られる．
1∑
(yi − f (xi ; α))2 .
n i=1
n
min
α
テストデータが学習データと同じ分布 p(y|x)pd (x) にしたがうとき，最小 2 乗法で推定し
た関数 f (x; α
b) を用いることで，テスト入力 x に対する出力 y を精度よく予測できると考
えられる．
共変量シフトでは，テストデータが学習データの分布 p(y|x)pd (x) とは異なる分布
p(y|x)pn (x) にしたがう問題設定を考える．ただし分布が全く異なる訳ではなく，条件付き
確率 p(y|x) は学習データとテストデータで共通であり，入力 x の分布が異なると仮定して
いる (図 1 左図)．現実の問題設定において，学習時とテスト時の入力分布が異なることは
珍しいことではない．教師ありデータ (x, y) を収集するコストが入力 x に依存する場合，
ラベル y を付けやすい (コストが低い) 入力データのみが教師ありデータとして観測されや
すくなる．このような場合には，テスト時の x の分布が学習時の分布とは異なると考えら
れる．
密度比によるダイバージェンス推定とその応用
図1
101
左図：pd (x) と pn (x) の確率密度関数のプロット．右図：最小 2 乗法による推定結果 (破線) と，密度比
による重み付き最小 2 乗法による推定結果 (実線)．点 (＋) は p(y|x)pd (x) から生成された学習データ．点 (○)
は p(y|x)pn (x) から生成された学習データ．
学習データとテストデータの分布が同じという仮定を置いて通常の推定を行うと，共変
量シフトの下では推定量に大きなバイアスが生じる．その理由を以下に示す．通常の最小
2 乗推定量の損失関数は，学習データ数が十分大きいとき
∫
(y − f (x; α))2 p(y|x)pd (x)dydx
(2.6)
に確率収束する．一方テストデータの分布 p(y|x)pn (x) の下では，回帰関数 f (x; α) による
予測の期待損失は
∫
(y − f (x; α))2 p(y|x)pn (x)dydx
(2.7)
で与えられる．損失 (2.6) の下で誤差が小さい回帰関数であっても，損失 (2.7) の下では損
失が大きくなることがある．したがって通常の最小 2 乗法では，共変量シフトの下で適切
な予測や推定を行うことができないと考えられる．
このバイアスを補正するために，pd (x) と pn (x) の密度比を用いた重み付き最小 2 乗法
1 ∑ pn (xi )
(yi − f (xi ; α))2
n i=1 pd (xi )
n
min
α
(2.8)
が有用である．学習データが分布 p(y|x)pd (x) にしたがうとき，(2.8) の重み付き 2 乗損失
は p(y|x)pn (x) の下での期待 2 乗誤差 (2.7) に確率収束する．したがって (2.8) による推定
102
日本統計学会誌第44巻第1号 2014
法により，p(y|x)pd (x) にしたがう学習データを用いて，p(y|x)pn (x) にしたがうテストデー
タに対してバイアスのない回帰関数の推定を行うことができる (図 1 右図)．確率密度関数
pd (x), pn (x) が未知であっても標本 (2.2) が利用可能なら，密度比 pn (x)/pd (x) の推定量を
用いて回帰関数の重み付き最小 2 乗推定量を得ることができる．
共変量シフト以外にも，密度比を用いた統計的推論の例として，外れ値検出や次元削減
などがある．詳細は Sugiyama et al. (2012) に詳しい．
密度比によるダイバージェンスの推定
3.
前節で紹介した密度比を，確率分布間のダイバージェンスの推定に用いることができる．
ダイバージェンスとは距離を一般化した概念であり，主に数理統計や情報理論などの分野
で，確率分布間の非類似度を測る尺度として広く用いられている．ダイバージェンスの推
定は 2 標本検定やデータの次元削減などさまざまな応用を持つため，重要な課題となって
いる．ダイバージェンスの重要なクラスのひとつである f -ダイバージェンスについて以下
で紹介する．
3.1
f -ダイバージェンス
密度比と関連が深いダイバージェンスである f -ダイバージェンス (Ali and Silvey (1966),
Csisz´ar (1967)) を紹介する．関数 f : R+ → R を f (1) = 0 を満たす強凸関数とする．確率
密度関数 pd (x), pn (x) の間の f -ダイバージェンスを
(
)
∫
pn (x)
Df (pd , pn ) = pd (x) f
dx
pd (x)
と定義する．関数 f の凸性と Jensen の不等式より
(∫
)
pn (x)
Df (pd , pn ) ≥ f
pd (x)
dx = f (1) = 0
pd (x)
となり，非負性が保証される．また f の強凸性を用いると Df (pd , pn ) = 0 から pd = pn が
導かれる．したがって Df (pd , pn ) の値は pd = pn のとき 0 であり，pd 6= pn のとき正の値
をとるため，Df (pd , pn ) は 2 つの確率分布の違いを定量的に表していると解釈できる．関
数 f (z) が z = 1 で微分可能なら，一般性を失うことなく f (1) = f 0 (1) = 0 を仮定すること
ができる．
f -ダイバージェンスの例を以下に挙げる．
例 3.1 (カルバック・ライブラーダイバージェンス) 関数 f (r) を f (r) = r − log r − 1
とすると
∫
Df (pd , pn ) =
pd (x) log
pd (x)
dx
pn (x)
密度比によるダイバージェンス推定とその応用
103
となる．これはカルバック・ライブラーダイバージェンスと呼ばれ，統計的推論における
最尤推定量と関連している．また情報理論の分野では符号長の冗長度を表す．
例 3.2 (全変動距離) 関数 f (r) を f (r) = |1 − r| とすると
∫
|pd (x) − pn (x)|dx
Df (pd , pn ) =
となる．これは確率密度関数の間の L1 距離であり，全変動距離とも呼ばれる．Pinsker の
不等式 (Pinsker (1964)) を通して，全変動距離はカルバック・ライブラーダイバージェン
スの下限を与える．
例 3.3 (相互情報量) ρ > 0 に対して関数 f (r) を
f (r) =
1+ρ
ρr
r(1 + ρ)
1
log
+
log
1+ρ
1 + ρr 1 + ρ
1 + ρr
と定義する．確率分布 pd (x), pn (x) に対して同時確率 p(x, y), y ∈ {n, d} を
p(x, d) = pd (x)
1
,
1+ρ
p(x, n) = pn (x)
ρ
1+ρ
と定義する．関数 f に対応する f -ダイバージェンスは
Df (pd , pn ) =
∫ ∑
p(x, y) log
y=n,d
p(x, y)
dx
p(x)p(y)
と表せる．ここで p(x), p(y) は p(x, y) の周辺分布を表す．これは同時分布 p(x, y) の相互
情報量に他ならない．確率分布 p(x, y) において，x と y が独立であることと pd = pn が成
り立つことは等価である．したがって pd = pn のとき，またそのときのみ Df (pd , pn ) = 0
となる．
3.2
f -ダイバージェンスの推定
f -ダイバージェンスは，密度比 r(x) = pn (x)/pd (x) を用いて
∫
Df (pd , pn ) = pd (x)f (r(x))dx
と表せる．この関係を用いて，密度比から f -ダイバージェンスを推定することができる．
まず関数 f の凸共役関数を用いてダイバージェンスを推定する方法を紹介する．関数
f (r) は r ≥ 0 に対して定義されているが，以下ではさらに r < 0 に対して f (r) = +∞ と
定義しておく．凸関数 f の凸共役関数 f ∗ は
f ∗ (w) = sup{rw − f (r)}
r∈R
104
日本統計学会誌第44巻第1号 2014
で与えられる．凸解析の結果から，真閉凸関数 f に対して (f ∗ )∗ = f が成り立つ (Rockafellar
(1970))．真閉凸関数 f (r) が r > 0 で微分可能なら，
f (r) = sup {rw − f ∗ (w)} = rf 0 (r) − f ∗ (f 0 (r)),
r>0
(3.1)
w∈R
となる．したがって，密度比モデルを r(x; θ) として (3.1) の w を f 0 (r(x; θ)) で置き換えれ
ば，r(x) = r(x; θ∗ ) より
∫
pd (x) sup {r(x; θ)f 0 (r(x; θ)) − f ∗ (f 0 (r(x; θ)))} dx
θ
{∫
}
∫
0
∗ 0
= sup
pn (x)f (r(x; θ))dx − pd (x)f (f (r(x; θ)))dx
Df (pd , pn ) =
θ∈Θ
と表せる．これを標本で近似して，Df (pd , pn ) の推定量
b sup
D


md
mn

 1 ∑
∑
1
(n)
(d)
f 0 (r(xi ; θ)) −
f ∗ (f 0 (r(xj ; θ)))
= sup

md j=1
θ∈Θ  mn i=1
(3.2)
を得る．このような推定法は Keziou and Leoni-Aubin (2005), Nguyen et al. (2010) など
で提案されている．
b sup との関連について調べる．パラメータ θ について最大
密度比推定と (3.2) の推定量 D
化することで f -ダイバージェンスを推定しているので，θ について微分すると，
md
mn
1 ∑
1 ∑
(n)
(n)
(d)
(d)
(d)
f 00 (r(xi ; θ))∇r(xi ; θ) −
r(xj ; θ)f 00 (r(xj ; θ))∇r(xj ; θ) = 0
mn i=1
md j=1
b sup では，関数
となる．したがって上の D
η(x; θ) = f 00 (r(x; θ))∇r(x; θ)
(3.3)
によるモーメント法で密度比を推定していると解釈できる．上のモーメント関数は，密度
比推定の最適なモーメント関数である ηopt (x; θ) とは一般に一致しない．モーメント関数
(3.3) が最適なモーメント関数 ηopt (x; θ) に一致するとき，f (r) は例 3.3 に示した相互情報
量に対応する関数になることが分かる．
凸共役関数を用いた f -ダイバージェンス推定の方法を一般化する (Kanamori et al.
(2012))．この推定法では，密度比推定とダイバージェンス推定を分けて考える．まず，
関数 η(x; θ) を用いたモーメント法で密度比を推定する．次に，推定された密度比を用いて
f -ダイバージェンスを推定する．そのために関数 f を，関数 fn (r), fd (r) を用いて
f (r) = rfn (r) + fd (r)
密度比によるダイバージェンス推定とその応用
105
と分解する．具体的な分解の仕方については 3.3 節で詳しく述べる．このとき f -ダイバー
ジェンスは，密度比 r(x) = pn (x)/pd (x) を用いて
∫
Df (pd , pn ) =
∫
pn (x)fn (r(x))dx +
pd (x)fd (r(x))dx
と表せる．推定した密度比を上式の r(x) に代入し，期待値を標本平均に置きかえることで，
f -ダイバージェンスを推定することができる．具体的には，モーメント関数 η(x; θ) による
b として，
密度比の推定量を r(x; θ)
md
mn
∑
1 ∑
(d) b
(n) b
b= 1
fd (r(xj ; θ)).
D
fn (r(xi ; θ))
+
mn i=1
md j=1
により Df (pd , pn ) を推定する．凸共役関数を用いた f -ダイバージェンスの推定は，密度比
推定にモーメント関数 (3.3) を用い，関数 f (r) を fn (r) = f 0 (r), fd (r) = −f ∗ (f 0 (r)) と分
解して推定する方法に対応している．
3.3
推定の最適性
前節で述べた f -ダイバージェンスの推定においてモーメント関数 η(x; θ) と分解 f (r) =
rfn (r) + fd (r) を適切に選べば，精度の高い推定量が得られると考えられる．以下では f ダイバージェンスに対する最適な推定法について考察する．
モーメント関数 η¯(x; θ) と関数 f の分割 f (r) = rf¯n (r)+ f¯d (r) を用いて，D = Df (pd , pn ) の
¯ が得られたとする．また，他のモーメント関数 η(x; θ) と分割 f (r) = rfn (r)+fd (r)
推定量 D
b とおく．密度比推定の最適性の場合と同様に，もし任意
を用いたときの D の推定量を D
b に対して D
¯が
の推定量 D
b − D,
¯ D]
¯ = o(1)
m · Cov[D
(3.4)
を満たすなら，推定分散に関して
¯ ≤ m · V [D]
b
m · V [D]
¯ が最適であることが分かる．した
が漸近的に成り立ち，考えている推定方式のなかで D
がって (3.4) を満たすようなモーメント関数 η¯ と分割 f¯n , f¯d を見付けることが目標となる．
√
¯ − D) と √m(D
b − D) を漸近展開
密度比推定の最適性を導出したときと同様に， m(D
して (3.4) を計算する．ベクトル cfn を
cfn = En [{f 0 (r(x; θ∗ )) − fn (r(x; θ∗ ))} ∇ log r(x; θ∗ )] ∈ Rd
106
日本統計学会誌第44巻第1号 2014
と定義し，行列 Uη は (2.4) で与えられるとすると
b − D,
¯ D
¯]
m(1 + ρ−1 ) · Cov[ D
[{
}
]
= En f¯n (r) − fn (r) + cTf¯n Uη¯−1 η¯ − cTfn Uη−1 η {f (r) − (r + ρ−1 )(f¯n (r) + cTf¯n Uη¯−1 η¯)} + o(1)
(3.5)
となる．上式右辺の第 1 項が 0 となる η¯ と f¯n , f¯d を求める．すでに密度比については最適
なモーメント関数 ηopt (x; θ) が得られているので，これを η¯(x; θ) とする．天下り的ではあ
るが，f の分割として
ρf (r)
,
f¯n (r) =
1 + ρr
f (r)
f¯d (r) =
1 + ρr
(3.6)
を考えると，(3.5) の第 1 項が 0 となることが確認できる．実際，(3.5) において b =
T
−ρUη¯−1 cf¯n ∈ Rd とおくと f (r) − (r + ρ−1 )(f¯n (r) + cTf¯n Uη¯−1 η¯) = ∇ log r(x; θ∗ ) b となり，
また
{
}
En [ f¯n (r) − fn (r) + cTf¯n Uη¯−1 η¯ − cTfn Uη−1 η ∇log rT b]
= En [(f¯n (r) − fn (r))∇log rT + cfT¯n − cTfn ]b
=0
b − D,
¯ D]
¯ = o(1) を得る．以上の結果は一般化することができ，次
となるので，m · Cov[D
の 2 条件
η¯(x; θ) = ηopt (x; θ),
−1
f (r(x; θ)) − (r(x; θ) + ρ
{
)f¯n (r(x; θ)) ∈ span
∂ log r(x; θ)
∂ log r(x; θ)
,...,
∂θ1
∂θd
}
(3.7)
が任意の θ ∈ Θ に対して成り立つとき，(3.5) が漸近的に 0 となることが分かる．2 つ目
の条件は，x の関数として左辺の関数が右辺の線形空間に属していることを表す．最適な
f -ダイバージェンス推定に関して，他の十分条件についても考察されている (Kanamori et
al. (2012))．
以下に具体例を示す．
例 3.4
ベクトル値関数 φ(x) = (φ1 (x), . . . , φd (x)) に対して密度比モデルを r(x; θ) =
exp{θT φ(x)}, θ ∈ Rd とおく．ここで φ1 (x) = 1 とする．正定数 ρ > 0 に対して f (r)
を例 3.3 で示した相互情報量に対応する関数とする．この f -ダイバージェンスを推定す
るとき，ηopt (x; θ) をモーメント関数として，fn (r) = ρf (r)/(1 + ρr), fd (r) = f (r)/(1 +
ρr) とすれば，漸近分散の意味で最適である．一方，凸共役関数から定まる分割 fn (r) =
f 0 (r), fd (r) = −f ∗ (f 0 (r)) は，密度比モデル r(x; θ) = exp{θT φ(x)} の下で f (r(x; θ)) −
107
密度比によるダイバージェンス推定とその応用
表1
b f の平均 2 乗誤差 mE[(D
b f − Df )2 ] を
カルバック・ライブラーダイバージェンス Df に対する推定量 D
示す．pn , pd はそれぞれ N (0, 1), N (µ, 1) の確率密度としている．標本数を md = mn = 100 として，µ の値ご
とに結果を示す．
µ
0.1
0.3
0.5
0.7
0.9
1.1
最適な推定
0.022
0.111
0.294
0.614
0.958
1.540
凸共役
0.023
0.119
0.321
0.719
1.187
2.094
(r(x; θ) + ρ−1 )f 0 (r(x; θ)) = −(1 + ρ)−1 θT φ(x) となるので，一般化された条件 (3.7) を満
たす．したがって凸共役関数から定まる推定法も，相互情報量の推定に関して最適である．
凸共役関数を用いる方法は Keziou and Leoni-Aubin (2005) で提案されている．
例 3.5
ベクトル値関数 φ(x) = (φ1 (x), . . . , φd (x)) に対して密度比モデルを r(x; θ) =
exp{θT φ(x)}, θ ∈ Rd とおく．ここで φ1 (x) = 1 とする．関数 f (r) を例 3.1 で示した関
数とすると，対応する f -ダイバージェンスはカルバック・ライブラーダイバージェンスと
なる．この場合には，凸共役関数を用いる推定法について f -ダイバージェンス推定の最
適性が保証されない．実際 fn (r) = f 0 (r) とすると，一般化された条件 (3.7) は eθ
T
φ(x)
∈
span{φ1 (x), . . . , φd (x)} と等価になる．一般にこの条件は成り立たないため，最適性のため
の十分条件は満たされない．最適な推定法と凸共役関数を用いた推定法を数値実験で比較
した結果を表 1 に示す．凸共役関数を用いた場合，平均 2 乗誤差が大きくなる傾向がある．
4.
2 標本検定への応用
f -ダイバージェンスの推定量を 2 標本検定に応用する．標本 (2.2) が得られたとき，次
の仮説検定を考える．
H0 : pd = pn ,
H1 : pd 6= pn .
ここで確率分布に対して密度比モデル r(x; θ) を仮定し，pn (x)/pd (r) = r(x; θ∗ ) が成り立
つとする．帰無仮説が正しいとき，密度比は r(x; θ∗ ) = 1 となる．密度比モデルの仮定は，
確率密度に対してはセミパラメトリックモデルを考えていることになる．
検定統計量として，密度比を用いた f -ダイバージェンスの推定量を用いる．推定値が 0
より十分大きいとき，帰無仮説を棄却する．以下の定理は帰無仮説の下での検定統計量の
漸近分布を与える．これより棄却域を構成することができる．
定理 4.1
モーメント関数 ηopt (x; θ) と f の分割 (3.6) から定まる Df (pd , pn ) の推定量
¯ f とする．正則条件を仮定するとき，帰無仮説 H0 の下で
をD
d − 1 のカイ 2 乗分布にしたがう．
2m ¯
f 00 (1) Dr
は漸近的に自由度
108
日本統計学会誌第44巻第1号 2014
自由度 d のカイ 2 乗分布の 100α パーセント点を χ2d (α) とする．帰無仮説 H0 の棄却域
¯f ≥
をD
f 00 (1) 2
2m χd−1 (1
− α) とすれば，有意水準 α の検定を行うことができる．
密度比モデルを用いた検定法としては，ワルド型の検定統計量が提案されている (Fokianos
et al. (2001))．密度比モデルを r(x; θ) = exp{θ1 + θ¯rT φ(x)} とすると，帰無仮説の下では
θ¯r = 0 が成り立つ．経験尤度最大化で密度比のパラメータ θ を推定する．帰無仮説の下で，
√
d
¯ −1 ) となる．ここで ∇r
¯ はパラ
θ¯r = 0 の推定量 θbr の漸近分布は mθbr −→ Nd−1 (0, Vn [∇r]
メータ θ¯r に対する r(x; θ) の θ = θ∗ = 0 ∈ Rd での勾配であり，d − 1 次元ベクトルとな
¯ の一致推定量 Vbn を用いて，経験尤度スコア S を S = mθbrT Vbn θbr と定義す
る．分散 Vn [∇r]
ると，S の漸近分布は自由度 d − 1 のカイ 2 乗分布であることが分かる．この結果を用い
て棄却域を構成することができる．S を用いた 2 標本検定は，確率分布に対するパラメト
リックモデルを仮定した t-検定や F -検定とほぼ同等の検出力をもつことが，さまざまな数
値実験から確認されている (Fokianos et al. (2001))．
一方，セミパラメトリック 2 標本検定に経験尤度比検定を用いることもできる．経験尤
度比 `(θ) を
`(θ) =
md
∑
i=1
log
1
+
(d)
ρr(xi ; θ) + 1
mn
∑
(n)
ρr(xj ; θ)
log
j=1
(n)
ρr(xj ; θ) + 1
とする．このとき経験尤度比検定統計量 R を
R = 2 sup{`(θ) − `(θ∗ )}
(4.1)
θ∈Θ
と定義する．帰無仮説のもとで R は自由度 d − 1 のカイ 2 乗分布にしたがう (Keziou and
b とすると，統計量 R
Leoni-Aubin (2008))．例 3.4 で示した凸共役関数による推定量を D
b と表せる．前節の結果より，経験尤度比検定統計量 R は相互情報量の推
は 2(md + mn )D
定量のなかで漸近分散を最小にしている．
¯ f を用いた 2 標本検定の検出力について考察する．
f -ダイバージェンスの最適な推定量 D
(m)
データ数の調和平均 m = 1/(1/mn + 1/md ) に対して，確率密度関数 pn (x) を
(m)
p(m)
n (x) = pd (x)r(x; θm ),
hm
θm = θ∗ + √ ,
m
hm ∈ R d
(4.2)
とおく．また hm → h ∈ Rd , m → ∞ とする．パラメータ θ∗ は r(x; θ∗ ) = 1 を満たして
いる．これは，データ数が増えるにしたがって 2 つの分布 pd , pn が近づくという，いわゆ
る “local alternative” の状況に相当し，非常に近い確率分布どうしを区別するときの精度
を考えていることになる．
定理 4.2
行列 M (θ) を Ed [∇ log r∇ log rT ] ∈ Rd×d とする．ここで関数値は θ = θ∗ で
評価している．確率変数 Y を自由度 d − 1，非心度 hT M (θ∗ )h の非心カイ 2 乗分布にした
密度比によるダイバージェンス推定とその応用
109
¯ f を用いて 2 標本検定を行うとき，local alternative の下で
がう確率変数とする．推定量 D
の検出力は漸近的に Pr(Y ≥ χ2d−1 (1 − α)) で与えられる．また，経験尤度スコア S を用い
た検定も，漸近的に同じ検出力をもつ．
¯f を
定理 4.2 より，密度比を用いたセミパラメトリック 2 標本検定で検定統計量として D
¯f と
用いると，検出力は関数 f の選び方に漸近的に依らないことが分かる．また提案法 D
経験尤度スコア S では 2 標本検定の精度に差はないことになる．
次に，密度比モデル r(x; θ) が真の密度比 pn (r)/pd (r) を含まない状況で，検出力を比較
(m)
する．確率密度関数 pn (x) を
(
)
sm (x) + εm
√
p(m)
(x)
=
p
(x)
r(x;
θ
)
+
d
m
n
m
とおく．ここで sm (x) は Ed [sm ] = 0 を満たす関数であり，また m → ∞ のとき εm → ε ∈ R
とする．さらにパラメータ θm は (4.2) で与えられるとする．このとき，検出力について次
の定理が成り立つ．
定理 4.3
確率変数 Y を自由度 d − 1，非心度 hT M (θ∗ )h − ε2 の非心カイ 2 乗分布とす
る．経験尤度スコア S を用いた検定の検出力は漸近的に Pr(Y ≥ χ2d−1 (1 − α)) で与えられ
¯ f による検定の検出力は漸近的に Pr(Y ≥ χ2 (1 − α) − ε2 )
る．また，定理 4.2 の推定量 D
d−1
となる．
カイ 2 乗分布の非心度について，pn (x) が確率密度となることから hT M (θ∗ )h − ε2 は
(m)
非負となる (Kanamori et al. (2012), 定理 4)．定理 4.3 より，密度比モデルが真の密度比
¯ f を用いるほうが経
を含まないときに，local alternative の状況で 2 標本検定を行うと，D
験尤度スコア S を用いる方法より検出力が高いことが分かる．実際のデータ解析では，厳
¯f
密にはモデルは真の密度比を含まないと考えらる．したがって定理 4.3 は，検定統計量 D
が現実的な状況で既存手法より高い精度を達成することを保証している．
数値例を示す．k 次元ユークリッド空間上の確率密度関数 pd (x), pn (x) に対して，密度比
モデルとして指数モデル


k
k


∑
∑
r(x; θ) = exp θ1 +
αi xi +
βj x2i ,


i=1
θ = (θ1 , α1 , . . . , αk , β1 , . . . , βk )
j=1
を仮定する．確率密度 pn (x) を表に示す分布とし，pd (x) を pd (x) = pn (x/σ)/σ と定め
る．ここで多次元確率分布の各要素は独立として密度関数を定義している．データ数を
mn = md = 100 として，カルバック・ライブラーダイバージェンスの推定量を用いる検定
法 (KL)，相互情報量の推定量を用いる検定法 (MI)，経験尤度スコアを用いる検定法 (S)，
ホテリング検定統計量を用いる方法 (Hote.) を比較する．
110
表2
日本統計学会誌第44巻第1号 2014
2 標本検定の検出力を示す．検定法にはカルバック・ライブラーダイバージェンス (KL)，相互情報量 (MI)，
経験尤度スコア (S)，ホテリング検定 (Hote.) を用いた．
mn = md = 1000
10 次元正規分布
10 次元 t-分布 (自由度 5)
σ
KL
MI
S
Hote.
KL
MI
S
Hote.
0.92
1.000
1.000
1.000
0.052
0.885
0.914
0.657
0.036
0.94
0.975
0.980
0.963
0.056
0.546
0.625
0.289
0.051
0.97
0.353
0.387
0.266
0.054
0.142
0.173
0.087
0.058
1.00
0.050
0.049
0.063
0.047
0.098
0.067
0.153
0.048
1.03
0.443
0.379
0.529
0.044
0.286
0.168
0.481
0.046
1.06
0.979
0.971
0.987
0.061
0.739
0.574
0.864
0.049
1.08
1.000
1.000
1.000
0.058
0.932
0.850
0.979
0.054
有意水準 5% の検定を行ったときの検出力を表 2 に示す．これより，ホテリング検定統
計量を用いる方法では，検出力が非常に低いことがか分かる．これは，分布の期待値が同
じであるためである．経験尤度スコアを用いる検定法では，とくに標本が t-分布にしたが
うとき，検出力に若干のバイアスがあることが分かる．最適な f -ダイバージェンス推定量
を用いる方法では，標本が正規分布にしたがうときには，KL と MI とも検出力が高くなっ
ている．一方，標本が t-分布にしたがうときには，相互情報量を用いる方法が優れている
と考えられる．カルバック・ライブラーダイバージェンスを用いる方法では，帰無仮説が
正しいときの検出力が有意水準 5% よりかなり大きくなっている．一方，相互情報量を用
いる方法では，帰無仮説が正しいときに検出力は有意水準 5% に近く，また σ 6= 1 のとき
の検出力が高くなっている．検出力に関する 1 次の漸近論では，最適な推定量を使う限り
検出力は f -ダイバージェンスには依存しないことが結論されていた．一方，数値的にはカ
ルバック・ライブラーダイバージェンスと相互情報量で検出力に差がある問題設定も存在
する．
5.
まとめ
本稿では密度比による f -ダイバージェンスの推定法を提案し，漸近分散の意味で最適な
推定法を導出した．また f -ダイバージェンスの推定量を 2 標本検定に応用し，漸近計算に
より検出力を導出した．1 次の漸近論では，検定の検出力は f -ダイバージェンスの f の選
び方には依存しないが，数値例では関数 f が検出力に影響していることが観測された．高
次漸近論などを適用するなどして，検定における関数 f の選択について考察を深める必要
がある．また，より広いクラスのダイバージェンスに対する推定量を構成し，その統計的
性質を理論的に追求することは重要な課題である．
密度比によるダイバージェンス推定とその応用
111
謝辞
本研究は JSPS 科研費 24500340 の助成を受けた．
参考文献
Ali, S. M. and Silvey, S. D. (1966). A general class of coeﬃcients of divergence of one distribution from another,
J. R. Stat. Soc. Ser. B , 28(1), 131–142.
Cover, T. M. and Thomas, J. A. (2006). Elements of Information Theory (Wiley Series in Telecommunications
and Signal Processing), Wiley-Interscience.
Csisz´ar, I. (1967). Information-type measures of diﬀerence of probability distributions and indirect observation,
Stud. Sci. Math. Hung., 2, 229–318.
Fokianos, K., Kedem, B., Qin, J. and Short, D. A. (2001). A semiparametric approach to the one-way layout,
Technometrics, 43, 56–64.
Kanamori, T., Suzuki, T. and Sugiyama, M. (2012). f -divergence estimation and two-sample homogeneity test
under semiparametric density-ratio models, IEEE Trans. Inf. Theory, 58(2), 708–720.
Keziou, A. and Leoni-Aubin, S. (2005). Test of homogeneity in semiparametric two-sample density ratio models,
Comptes Rendus Mathematique, 340(12), 905–910.
Keziou, A. and Leoni-Aubin, S. (2008). On empirical likelihood for semiparametric two-sample density ratio
models, J. Stat. Plann. Inference, 138(4), 915–928.
Kullback, S. and Leibler, R. A. (1951). On information and suﬃciency, Ann. Math. Stat., 22, 79–86.
Nguyen, X., Wainwright, M. J. and Jordan, M. I. (2010). Estimating divergence functionals and the likelihood
ratio by convex risk minimization, IEEE Trans. Inf. Theory, 56(11), 5847–5861.
Pinsker, M. S. (1964). Information and Information Stability of Random Variables and Processes, Holden-Day.
Qin, J. (1998). Inferences for case-control and semiparametric two-sample density ratio models, Biometrika,
85(3), 619–639.
Rockafellar, R. T. (1970). Convex Analysis, Princeton University Press.
Shimodaira, H. (2000). Improving predictive inference under covariate shift by weighting the log-likelihood
function, J. Stat. Plann. Inference, 90(2), 227–244.
Sugiyama, M. and Kawanabe, M. (2012). Machine Learning in Non-Stationary Environments: Introduction to
Covariate Shift Adaptation, Adaptive computation and machine learning, MIT Press.
Sugiyama, M., Suzuki, T. and Kanamori, T. (2012). Density Ratio Estimation in Machine Learning, Cambridge
University Press.

Download Report