ロバスト推測の基礎とダイバージェンス型メソッドへの発展

藤澤 洋徳
統計数理研究所
fujisawa@ism.ac.jp

1

Contents
第一部：ロバスト推測の基礎
1. 外れ値とロバスト推測
2. 簡単なロバスト推定
第二部：ダイバージェンス型メソッドへの発展
3. 重み付きスコアに基づいたロバスト法と関連したダイバージェンス
4. γ -ダイバージェンスに基づいたロバスト法の性質

2

1. 外れ値とロバスト推測
外れ値とは何であるのか？
ロバスト推測とは何であるのか？

3

ある実験をしていて，以下のような 10 個のデータ値が得られたとする：

5.6, 5.7, 5.4, 5.5, 5.8, 5.5, 5.3, 5.6, 5.4, 5.2.
このデータの標本平均は 5.5 である：

5.6 + 5.7 + 5.4 + 5.5 + 5.8 + 5.5 + 5.3 + 5.6 + 5.4 + 5.2= 5.5.
10
何らかの原因で，外れ値（異常値，outlier）を観測することがある．
最後のデータ値がタイプされる時に，誤って 5 が二回タイプされてしまったとする．
5.6, 5.7, 5.4, 5.5, 5.8, 5.5, 5.3, 5.6, 5.4, 55.2.
このときの標本平均は 10.5 になる．

5.6 + 5.7 + 5.4 + 5.5 + 5.8 + 5.5 + 5.3 + 5.6 + 5.4 + 55.2
10

= 10.5.
4

5.6, 5.7, 5.4, 5.5, 5.8, 5.5, 5.3, 5.6, 5.4, 5.2.
5.6, 5.7, 5.4, 5.5, 5.8, 5.5, 5.3, 5.6, 5.4, 55.2.
外れ値の混入標本平均： 5.2： 5.5→→55.2.
10.5.

何が起きているのか？

x̄ =x1 + x2 + · · · + xn−1 + xn
n→∞as xn → ∞5

標本平均 10.5 は，本来の意図するところを見ているとは，とても言えない．
データに外れ値が入っている場合には，何らかの対策が必要となりやすい．
外れ値に影響されにくい推定をロバスト推定 (robust estimation) といい，外れ値に影響されにくい統計的推測全体をロバスト推測 (robust inference) という．
注意： なお，ロバスト推測という言葉は，広義の意味では，外れ値に限らず，何らかの悪影響に引きずられにくい統計的推測全体のことを言う．

6ρ̂ = 0.44ρ̂∗ = 0.91

ρ̂rob = 0.85
7

Hampel et al. (1986)

8

2. 簡単なロバスト推定
外れ値にロバストである推定方法について，幾つかの簡単な例を挙げる．

9

平均 μ のロバスト推定
5.6, 5.7, 5.4, 5.5, 5.8, 5.5, 5.3, 5.6, 5.4, 55.2.
まずは平均 μ をロバスト推定することを考えよう．
最も単純で汎用的に用いられている推定値は，中央値 (median) である．
データ x1, . . . , xn を小さい順に並び変えた順序統計値を以下で表す：

x(1) ≤ x(2) ≤ · · · ≤ x(n)このとき，中央値は，次で定義される：

Med({xi}) = x(k)= (x(k−1) + x(k))/2

n = 2k − 1 のとき
n = 2k のとき
10タイプミスのある例：

5.6, 5.7, 5.4, 5.5, 5.8, 5.5, 5.3, 5.6, 5.4, 55.2.
順序統計値に並び変える：

5.3, 5.4, 5.4, 5.5, 5.5, 5.6, 5.6, 5.7, 5.8, 55.2.
中央値は次になる：

Med({xi}) = (5.5 + 5.6)/2 = 5.55
これは平均 μ の妥当な推定値であろう．
11ただし，中央値を使うときには，いくつかの注意が必要である．
外れ値の割合が大きい場合は，中央値は平均の妥当な推定値ではない．

0.0

0.3
0.0

0.1

0.1

0.2

0.2

dlnorm (x)

0.4

0.3

0.5

0.4

0.6

0.5

データの背後にある母集団は対称分布である．

−4−20x2

4

02

4x68

Median12

標準偏差 σ のロバスト推定
5.6, 5.7, 5.4, 5.5, 5.8, 5.5, 5.3, 5.6, 5.4, 5.2.
5.6, 5.7, 5.4, 5.5, 5.8, 5.5, 5.3, 5.6, 5.4, 55.2.

S2 =1nn−1i=1S 2 = 0.185(xi − x̄)2
S 2 = 15.7

σ̂rob(= MADN) = 0.222
13平均 μ の他の簡単なロバスト推定法
刈り込み平均 (trimmed mean)x(1) ≤ · · · ≤ x[nα] ≤ · · · ≤ xn−[nα] ≤ · · · ≤ x(n)上側 100α%と下側 100α%のデータを使わない平均．

1μ̂ =n − 2[nα]n−[nα]

x(i)i=[nα]+1

ただし [a] は a を超えない整数．

14重み付きスコアに基づいたロバスト推定
母集団が N (μ, σ 2 ) に従う場合を考える．
簡単のために標準偏差 σ は既知であるとする．
平均パラメータ μ を推定する問題を考えよう．

最尤推定量 尤度方程式n(xi − μ) = 0(μ = x̄)i=115平均 μ をロバスト推定することを考えよう．
正規分布の密度関数を φ(x; μ, σ 2 ) で表す．
2φ(x; μ, σ ) = √12πσ 2exp −(x − μ)2

2σ2

外れ値の典型的な性質の一つは，その生起確率が小さいことである．
外れ値 x∗ に対して φ(x∗; μ, σ 2 ) は小さくなるだろう．

16

尤度方程式において，外れ値の寄与を小さくするために，平均 μ の推定値 μ̂ を，次の推定方程式の解として与えることにしよう：nφ(xi; μ, σ 2)(xi − μ) = 0.

i=1

最初の例で，タイプミスがあった場合に，上記の推定方法を適用してみる．ただし，標準偏差 σ は，外れ値を除外したデータの標本標準偏差の近似値 0.158 であったと仮定してみる．そのとき，推定値は，μ̂ = 5.517 となった．

17

疑問
もしも σ 2 が未知だったら？
もしも正規分布（対称分布）でなくて一般の分布 f (x; θ) だったら？
第二部で関連方法が提案される．

18

3. 重み付きスコアに基づいたロバスト法と関連したダイバージェンス
本節と次節では，データ x に対して，何らかのパラメトリックモデル fθ (x) = f (x; θ)が想定されている．

Basu et al. (1998, Biometrika)
Eguchi and Kano (1998, unpublished)
Jones et al. (2001, Biometrika)

19平均パラメータ μ の重み付き推定方程式を思い出そう：nφ(xi; μ, σ 2 )(xi − μ) = 0i=1

コアはスコア関数だと考えよう：ns(xi; θ) = 0s(x; θ) =i=1ddθlog f (x; θ)

重み付き推定方程式を考えてみよう：nf (xi; θ)β s(xi; θ) = 0

β>0i=1

正規分布における平均パラメータの推定の場合は，このままで良かったが，一般のパラメータを推定する場合には，もう一工夫が必要である．
20重み付き推定方程式：nf (xi; θ)β s(xi; θ) = 0.

i=1推定方程式の不偏性を思い出す：βEfθ f (x; θ) s(x; θ) = 0 とは限らない．
）（正規分布での平均 μ の推定の時は自動的に満たされていた．
核の関数を修正する：

ψ(x; θ) = f (x; θ)β s(x; θ) − Efθ f (x; θ)β s(x; θ)

これによって推定方程式の不偏性が満たされる： Efθ [ψ(x; θ)] = 0.
21M 推定：10=ψ(xi; θ)nni=1

重み付きスコアに基づいた推定方程式：n 
1
f (xi; θ)β s(xi; θ) − Efθ f (x; θ)β s(x; θ)0 =n1=ni=1nf (xi; θ)β s(xi; θ) −f (x; θ)1+β s(x; θ)dxi=122尤度方程式と KL ダイバージェンスの関係を思い出そう：n1θ̂ = arg max
log f (xi; θ)θ ni=1

これを微分すると尤度方程式が得られる：

0=n1

ns(xi; θ)i=1

右辺を積分して戻せば元に戻る：

1log f (xi; θ)nni=1

これの極限にマイナスをつけたものを考える（相互エントロピー）．
23相互エントロピー：

dKL(g, fθ ) = −Eg [log f (x; θ)] = −g(x) log f (x; θ)dx

KL ダイバージェンスは次のように表現できる：

DKL(g, f ) = Eg logg(X)= Eg [log g(X)] − Eg [log f (X)]f (X)= −dKL(g, g) + dKL (g, f )ダイバージェンスの性質 (≈ 距離)

DKL(g, f ) ≥ 0.

DKL(g, f ) = 0 ⇔ g = f.
24

β -ダイバージェンス
重み付きスコアに基づく推定方程式から，相互エントロピーとダイバージェンスを作る：

0=n1

nf (xi; θ)β s(xi; θ) −f (x; θ)1+β s(x; θ)dxi=1

KL のときと同様に積分して極限を取ってマイナスをつける．
β -相互エントロピー：

1dβ (g, fθ ) = −β1g(x)f (x; θ) dx +1+ββf (x; θ)1+β dx

25β -ダイバージェンス (β -divergence, density power divergence)：Dβ (g, f ) = −dβ (g, g) + dβ (g, f )1
1=g(x)1+β dx −β(1 + β)1+β1f (x)1+βdx+1+βg(x)f (x)β dx

極限を考えると KL ダイバージェンスになる：

lim Dβ (g, f ) = DKL(g, f ).

β→026

γ -ダイバージェンス
重み付きスコアに基づく推定方程式を見直す：
n1
0=f (xi; θ)γ s(xi; θ) − Efθ [f (x; θ)γ s(x; θ)]ni=1

重みを 1 にする：

0=n1nf (xi; θ)γ s(xi; θ)

i=1−Efθ [f (x; θ)γ s(x; θ)]n1

nf (xi; θ)γi=1Efθ [f (x; θ)γ ]

先ほどと同じように行うと，対応する γ -相互エントロピーが得られる：

1dγ (g, fθ ) = − logγγg(x)f (x; θ) dx +

11+γ

log

f (x; θ)1+γ dx
27

γ -ダイバージェンス：
Dγ (g, f ) = −dγ (g, g) + dγ (g, f )111+γ
log g(x)logdx −=γ(1 + γ)1+γ1log f (x)1+γ dx+1+γg(x)f (x)γ dx

極限を考えると KL ダイバージェンスになる：

lim Dγ (g, f ) = DKL(g, f ).

γ→028

二つの相互エントロピーの違い
β -相互エントロピー：
dβ (g, f ) = −1ββg(x)f (x) dx +11+βf (x)1+β dx

γ -相互エントロピー：
1dγ (g, f ) = − logγγg(x)f (x) dx +11+γlogf (x)1+γ dx

注意： 二つの相互エントロピーの違いは「log」が付いているかどうかである．
これだけの違いだが大きな性質の違いも生み出す．
29

4. γ -ダイバージェンスに基づいたロバスト法の性質
Fujisawa and Eguchi (2008). Robust parameter estimation witha small bias against heavy contamination. Journal of Multivariate Analysis, Vol.99, 2053-2081.

30外れ値とは何か
データを発生している分布（汚染された分布）

g(x) = (1 − ε)f (x) + εδ(x)f : 目的分布 (= fθ∗ )δ : 外れ値の分布（汚染分布）
ε: 外れ値の割合
外れ値の分布が目的分布の裾にある：

νf =f (x)δ(x)dx ≈ 0f (x)γ0 δ(x)dx

1/γ0


≈0

γ0 > γ > 031

唯一つの本質的な仮定
外れ値の分布が目的とする分布の裾にある：

νf =f (x)γ0 δ(x)dx

1/γ0

≈0

外れ値を x∗ として，外れ値の分布を δ(x) = δx∗ (x) とすると，

νf = f (x∗) ≈ 0となるので，想定された条件は，外れ値の生起確率が小さいことを意味する．

32

Review: 経験推定
δa (x)

ディラック関数f (x)δa(x)dx = f (a)経験密度関数(厳密には経験分布関数で議論すべき)

1
ḡ(x) =
δxi (x) → g(x)
nni=1

経験推定

Eg [h(X)] =g(x)h(x)dx ←n1ḡ(x)h(x)dx =h(xi)
ni=1

33

γ -相互エントロピー：
1dγ (g, f ) = − log
γg(x)f (x)γ dx +

1
log
1+γ

f (x)1+γ dx

γ 相互エントロピーの経験推定
11γlogdγ (ḡ, f ) = − log ḡ(x)f (x) dx +γ
⎞ 1+γ
⎛1n1

= − log ⎝γ
ni=1f (xi)γ⎠+11+γlog

f (x)1+γ dx
f (x)1+γ dx

34

(γ -) ロバスト推定量θ̂γ= arg min dγ (ḡ, fθ )θ= arg min −dγ (ḡ, ḡ) + dγ (ḡ, fθ ) = arg min Dγ (ḡ, fθ )θ→θθγ∗ = arg min dγ (g, fθ )θ

潜在的なバイアス
θγ∗ − θ ∗(≈ 0)= O (εν γ )35

ピタゴリアン関係
Dγ (g, f ) = − dγ (g, g) + dγ (g, f )Dγ (g, fθ ) = Dγ (g, f ) + Dγ (f, fθ ) + O(εν γ )­  ­    Æ ­  36

潜在的なバイアス (ε = 0.05)
g = (1 − ε)N (0, 1) + εN (6, 1)
μσ0.3

0.6

0.25

0.5

0.2

0.4

0.15

0.3

0.1

0.2

0.05

0.1
0.20.4fθ = N (μ, σ 2 )

←0.6

0.8

1

θγ∗ :(m)

θβ
0.2

0.4

0.6

0.8

1

x 軸:
y 軸:

実線 dγ (g, f ):点線 dβ (g, f )γ or β
潜在バイアス37

潜在的なバイアス (ε = 0.2)
g = (1 − ε)N (0, 1) + εN (6, 1)
μfθ = N (μ, σ 2 )←σ1.2θγ∗ :

1.5

1

1.25

0.8

10.60.75

0.4

0.50.2

0.25
0.2

0.4

0.60.8

1(m)

θβ0.2

0.4

0.60.81x 軸:
y 軸:

実線 dγ (g, f ):点線 dβ (g, f )γ or β
潜在バイアス38

更新アルゴリズム
パラメトリック分布を N (μ, σ 2) としたときは，次の繰り返しアルゴリズムの収束値によって推定値が得られます： θ = (μ, σ 2)．
(a)wi= f (xi; θ (a))γμ(a+1) =(σ 2)(a+1) =n


f (xi; θ (a) )γi=1(a)

w i xii=1⎧n⎨⎩n

i=1(a)⎫⎬

wi x2i − (μ(a+1))2 (1 + γ)⎭

このアルゴリズムは次のような単調性をもちます：dγ (ḡ, fθ(a) ) ≥ dγ (ḡ, fθ(a+1) ) ≥ · · · ≥ dγ (ḡ, fθ̂ )
γ

39

漸近分散
漸近正規性：

√ dn θ̂γ − θγ∗ → N 0, Σg (θγ∗ )Jg (θ) = Egddθ

ξ(x; θ)Σg (θ) = Jg (θ)−1Ig (θ)Jg (θ)
Ig (θ) = Eg ξ(x; θ)ξ(x; θ)−1!

漸近分散の関係：Σg (θγ∗ ) =1
Σf (θ ∗) + O(εν γ )1−ε

結論： 外れ値は自動的に無視されている．
40

ある種のロバスト性を生じる相互エントロピーの唯一性
外れ値の割合が大きい場合にもバイアスが小さくなるロバスト推定をもたらす相互
エントロピー d(g, f ) は，幾つかの条件の下では，本質的に唯一つである：

d(g, f ) = φ(dγ (g, f )).
ただし φ(u) は適当な単調増加関数である．

41

THANK YOU

Hironori Fujisawa
Institute of Statistical Mathematics
42