コインを 3 回投げて 3 回とも表が出たとき、「このコインは絶対に裏が出ない」と結論してよいでしょうか。直感では「まだ決めつけられない」と感じるはずです。この記事は、その「決めつけない」感覚を、確率の言葉で正確に書く方法を扱います。同じ考え方で、機械学習でおなじみの「重みを小さく抑える工夫」(正則化)の正体も分かり、さらに予測が「どのくらい当てになるか」まで数字で出せるようになります。
**最尤推定 (MLE)**は「観測データが最も出やすくなるパラメータ」を選ぶ方法です。データが少ないと極端な答え(確率 1 1 1 や無限大の重み)を平然と返します。
ベイズの定理 は、データを見る前の信念(事前分布)とデータの当てはまり(尤度)を掛けて、データを見た後の信念(事後分布)を作る規則です。MAP 推定 は事後分布の最頻値を取ることです。
これまで「過学習を防ぐおまじない」として導入してきた正則化は、事前分布そのもの です。L2 正則化はガウス事前分布の、L1 正則化はラプラス事前分布(原点でとがった山の形をした分布)の MAP 推定に一致し、正則化係数は λ = σ 2 / τ 2 \lambda = \sigma^2/\tau^2 λ = σ 2 / τ 2 (観測ノイズの分散と事前分散の比)という意味を持ちます。
事後分布を 1 点に潰さず、パラメータについて積分すると事後予測分布 が得られます(新しい入力に対する出力を、1 つの値ではなく分布として返すものです)。その分散は「観測ノイズ」と「パラメータの不確実性」の和に分解され、外挿するほど後者が効きます。
MAP は事後分布のたった 1 点であり、座標変換で不変ではありません。L1 のスパース性も「最頻値(分布の山の頂上の位置)の性質」であって「事後平均の性質」ではありません。
この記事で使う事前知識
最尤推定 : 観測データが最も出やすくなるようにパラメータを選ぶ方法。→ 最尤推定
最小二乗問題 : 予測と観測の差(残差)の 2 乗和を最小にする重みを求める問題。→ 最小二乗問題
交差エントロピー誤差 : ロジスティック回帰で使う、負の対数尤度としての損失。→ 交差エントロピー誤差
凸関数 : 2 点を結ぶ線分が常にグラフの上側にある関数。→ 凸関数
線形回帰と最小二乗法 では、残差平方和 ∑ i ( y i − w T x i ) 2 \sum_i (y_i - \boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i)^2 ∑ i ( y i − w T x i ) 2 を最小にする w \boldsymbol{w} w を求めました(最小二乗問題(定義 3.1)[線形回帰と最小二乗法] )。ロジスティック回帰 では交差エントロピー(交差エントロピー誤差(定義 4.1)[ロジスティック回帰] )を最小にしました。そして過学習を抑えるために、どちらでも罰則項 λ ∥ w ∥ 2 \lambda \lVert \boldsymbol{w}\rVert^2 λ ∥ w ∥ 2 を足しました。
ここで素朴な疑問が 3 つ残ります。
なぜ残差の 2 乗 なのですか。絶対値でも 4 乗でもよいはずです。
なぜ罰則が ∥ w ∥ 2 \lVert \boldsymbol{w}\rVert^2 ∥ w ∥ 2 なのですか。λ \lambda λ はどうやって決めるのですか。単位は何ですか。
モデルが返すのは y ^ = w ^ T x \hat{y} = \hat{\boldsymbol{w}}^{\mathsf{T}}\boldsymbol{x} y ^ = w ^ T x という 1 つの数だけです。「この予測はどのくらい信用してよいか」はどこに書いてあるのですか。
この章の主張は、3 つとも確率モデルを明示すれば同時に答えが出る 、というものです。1 の答えは「観測ノイズをガウス分布と仮定したから」、2 の答えは「重みの事前分布をガウス分布と仮定したから、λ \lambda λ はノイズ分散と事前分散の比」、3 の答えは「点推定を捨てて事後分布のまま積分すればよい」です。
歴史的には、この見方は機械学習より 250 年ほど古いものです。トーマス・ベイズの遺稿(1763 年)と、それを独立に一般の形で展開したラプラスの 1774 年の論文は、「結果から原因の確率を測る」問題として条件付き確率の反転を扱いました。一方、20 世紀初頭の R. A. フィッシャーは事前分布を持ち込むことを嫌い、尤度だけを使う最尤推定を統計学の中心に据えました。現代の機械学習はこの両方を、目的に応じて使い分けています。損失関数の設計は前者の言葉で、汎化性能の議論は後者の言葉で語られることが多い、という具合です。
観測データを D \mathcal{D} D 、モデルのパラメータを θ \theta θ (ベクトルのときは w \boldsymbol{w} w )と書きます。確率モデルとは、θ \theta θ を決めるごとに D \mathcal{D} D の生成される確率(密度)p ( D ∣ θ ) p(\mathcal{D} \mid \theta) p ( D ∣ θ ) が定まる仕組みのことです。
定義 2.1 (尤度・事前分布・事後分布 ) ↗
確率モデル p ( D ∣ θ ) p(\mathcal{D} \mid \theta) p ( D ∣ θ ) において、データ D \mathcal{D} D を固定し、θ \theta θ の関数と見た もの
L ( θ ) = p ( D ∣ θ ) L(\theta) = p(\mathcal{D} \mid \theta) L ( θ ) = p ( D ∣ θ ) を尤度関数 といい、その対数 ℓ ( θ ) = log L ( θ ) \ell(\theta) = \log L(\theta) ℓ ( θ ) = log L ( θ ) を対数尤度 といいます。とくに D = ( z 1 , … , z N ) \mathcal{D} = (z_1,\ldots,z_N) D = ( z 1 , … , z N ) が独立同分布のとき、
L ( θ ) = ∏ i = 1 N p ( z i ∣ θ ) , ℓ ( θ ) = ∑ i = 1 N log p ( z i ∣ θ ) L(\theta) = \prod_{i=1}^{N} p(z_i \mid \theta), \qquad \ell(\theta) = \sum_{i=1}^{N} \log p(z_i \mid \theta) L ( θ ) = i = 1 ∏ N p ( z i ∣ θ ) , ℓ ( θ ) = i = 1 ∑ N log p ( z i ∣ θ ) となります。
パラメータ θ \theta θ 自身を確率変数と見なし、データを観測する前の分布 p ( θ ) p(\theta) p ( θ ) を事前分布 、観測後の条件付き分布 p ( θ ∣ D ) p(\theta \mid \mathcal{D}) p ( θ ∣ D ) を事後分布 といいます。また
p ( D ) = ∫ p ( D ∣ θ ) p ( θ ) d θ p(\mathcal{D}) = \int p(\mathcal{D} \mid \theta)\, p(\theta)\, d\theta p ( D ) = ∫ p ( D ∣ θ ) p ( θ ) d θ を周辺尤度 (またはエビデンス)といいます。θ \theta θ が離散なら積分は和に置き換えます。
尤度は「θ \theta θ についての確率分布」ではないことに注意してください。∫ L ( θ ) d θ \int L(\theta)\,d\theta ∫ L ( θ ) d θ は一般に 1 1 1 になりません。L L L はあくまで「この θ \theta θ を信じたとき、手元のデータはどれくらい出やすかったか」を測る物差しです。
定理 2.2 (ベイズの定理 ) ✓
θ \theta θ と D \mathcal{D} D の同時密度 p ( θ , D ) p(\theta, \mathcal{D}) p ( θ , D ) が存在し、p ( D ) > 0 p(\mathcal{D}) > 0 p ( D ) > 0 であるとする。このとき
p ( θ ∣ D ) = p ( D ∣ θ ) p ( θ ) p ( D ) = p ( D ∣ θ ) p ( θ ) ∫ p ( D ∣ θ ′ ) p ( θ ′ ) d θ ′ p(\theta \mid \mathcal{D}) = \frac{p(\mathcal{D} \mid \theta)\, p(\theta)}{p(\mathcal{D})}
= \frac{p(\mathcal{D} \mid \theta)\, p(\theta)}{\displaystyle\int p(\mathcal{D} \mid \theta')\, p(\theta')\, d\theta'} p ( θ ∣ D ) = p ( D ) p ( D ∣ θ ) p ( θ ) = ∫ p ( D ∣ θ ′ ) p ( θ ′ ) d θ ′ p ( D ∣ θ ) p ( θ ) が成り立つ。とくに D \mathcal{D} D を固定すれば分母は θ \theta θ に依存しない定数なので、
p ( θ ∣ D ) ∝ p ( D ∣ θ ) p ( θ ) p(\theta \mid \mathcal{D}) \propto p(\mathcal{D} \mid \theta)\, p(\theta) p ( θ ∣ D ) ∝ p ( D ∣ θ ) p ( θ ) すなわち「事後分布 ∝ \propto ∝ 尤度 × \times × 事前分布」である。
証明(定理 2.2) 条件付き密度の定義から、p ( D ) > 0 p(\mathcal{D}) > 0 p ( D ) > 0 のとき
p ( θ ∣ D ) = p ( θ , D ) p ( D ) p(\theta \mid \mathcal{D}) = \frac{p(\theta, \mathcal{D})}{p(\mathcal{D})} p ( θ ∣ D ) = p ( D ) p ( θ , D ) です。同じ定義を逆向きに使うと、p ( θ ) > 0 p(\theta) > 0 p ( θ ) > 0 となる θ \theta θ について p ( D ∣ θ ) = p ( θ , D ) / p ( θ ) p(\mathcal{D} \mid \theta) = p(\theta, \mathcal{D}) / p(\theta) p ( D ∣ θ ) = p ( θ , D ) / p ( θ ) 、すなわち p ( θ , D ) = p ( D ∣ θ ) p ( θ ) p(\theta, \mathcal{D}) = p(\mathcal{D} \mid \theta)\, p(\theta) p ( θ , D ) = p ( D ∣ θ ) p ( θ ) です(p ( θ ) = 0 p(\theta) = 0 p ( θ ) = 0 の点では両辺とも 0 0 0 なのでこの等式は全体で成り立ちます)。これを最初の式の分子に代入すれば第 1 の等号が得られます。
第 2 の等号は、定義 2.1 の周辺尤度の定義そのものです。実際、同時密度を θ \theta θ について積分すれば D \mathcal{D} D の周辺密度になり、
p ( D ) = ∫ p ( θ ′ , D ) d θ ′ = ∫ p ( D ∣ θ ′ ) p ( θ ′ ) d θ ′ p(\mathcal{D}) = \int p(\theta', \mathcal{D})\, d\theta' = \int p(\mathcal{D} \mid \theta')\, p(\theta')\, d\theta' p ( D ) = ∫ p ( θ ′ , D ) d θ ′ = ∫ p ( D ∣ θ ′ ) p ( θ ′ ) d θ ′ となります。最後の比例関係は、分母が θ \theta θ を含まないことから直ちに従います。
∎
比例関係のほうが実用上は重要です。事後分布を求める作業の大半は、「尤度 × \times × 事前分布」を θ \theta θ の関数として書き下し、θ \theta θ を含まない因子をすべて捨てて 、残った形が何という分布かを見抜くことに尽きます。
例 2.3 (事前分布が効く例:まれな病気の検査 )
有病率 0.1 % 0.1\% 0.1% の病気があり、検査の感度(病気の人が陽性になる確率)が 99 % 99\% 99% 、特異度(健康な人が陰性になる確率)が 95 % 95\% 95% だとします。ある人の検査が陽性でした。この人が病気である確率はいくらでしょうか。
θ ∈ { 病気 , 健康 } \theta \in \{\text{病気}, \text{健康}\} θ ∈ { 病気 , 健康 } 、D = 陽性 \mathcal{D} = \text{陽性} D = 陽性 とします。事前分布は p ( 病気 ) = 0.001 p(\text{病気}) = 0.001 p ( 病気 ) = 0.001 、p ( 健康 ) = 0.999 p(\text{健康}) = 0.999 p ( 健康 ) = 0.999 。尤度は p ( 陽性 ∣ 病気 ) = 0.99 p(\text{陽性} \mid \text{病気}) = 0.99 p ( 陽性 ∣ 病気 ) = 0.99 、p ( 陽性 ∣ 健康 ) = 1 − 0.95 = 0.05 p(\text{陽性} \mid \text{健康}) = 1 - 0.95 = 0.05 p ( 陽性 ∣ 健康 ) = 1 − 0.95 = 0.05 です。定理 2.2 より
p ( 陽性 ) = 0.99 × 0.001 + 0.05 × 0.999 = 0.00099 + 0.04995 = 0.05094 , p ( 病気 ∣ 陽性 ) = 0.00099 0.05094 = 0.01943 … ≈ 1.9 % . \begin{aligned}
p(\text{陽性}) &= 0.99 \times 0.001 + 0.05 \times 0.999 = 0.00099 + 0.04995 = 0.05094,\\
p(\text{病気} \mid \text{陽性}) &= \frac{0.00099}{0.05094} = 0.01943\ldots \approx 1.9\%.
\end{aligned} p ( 陽性 ) p ( 病気 ∣ 陽性 ) = 0.99 × 0.001 + 0.05 × 0.999 = 0.00099 + 0.04995 = 0.05094 , = 0.05094 0.00099 = 0.01943 … ≈ 1.9%. 検査の精度が高いにもかかわらず、陽性でも病気である確率は 2 % 2\% 2% 程度です。尤度比は 0.99 / 0.05 = 19.8 0.99/0.05 = 19.8 0.99/0.05 = 19.8 倍ありますが、事前オッズが 1 : 999 1:999 1 : 999 と極端に小さいので、事後オッズは 19.8 / 999 ≈ 1 / 50.5 19.8/999 \approx 1/50.5 19.8/999 ≈ 1/50.5 にしかならないためです。
「尤度だけを見て θ \theta θ を選ぶ」(この場合は p ( 陽性 ∣ θ ) p(\text{陽性}\mid\theta) p ( 陽性 ∣ θ ) が大きい「病気」を選ぶ)という判断が、事前分布を無視したときにどれほど誤るかを示す例です。
定義 3.1 (最尤推定量 ) ↗
パラメータ空間 Θ \Theta Θ 上の確率モデル p ( D ∣ θ ) p(\mathcal{D} \mid \theta) p ( D ∣ θ ) に対し、尤度関数 L ( θ ) = p ( D ∣ θ ) L(\theta) = p(\mathcal{D}\mid\theta) L ( θ ) = p ( D ∣ θ ) を最大にする θ \theta θ
θ ^ M L = arg max θ ∈ Θ p ( D ∣ θ ) = arg max θ ∈ Θ log p ( D ∣ θ ) \hat{\theta}_{\mathrm{ML}} = \operatorname*{arg\,max}_{\theta \in \Theta} p(\mathcal{D} \mid \theta)
= \operatorname*{arg\,max}_{\theta \in \Theta} \log p(\mathcal{D} \mid \theta) θ ^ ML = θ ∈ Θ arg max p ( D ∣ θ ) = θ ∈ Θ arg max log p ( D ∣ θ ) を最尤推定量 といいます。log \log log は狭義単調増加なので、L L L を最大にする点と ℓ = log L \ell = \log L ℓ = log L を最大にする点は一致します。
対数を取る理由は 2 つあります。独立同分布の積が和に変わって微分しやすくなること、そして N N N が大きいときに L ( θ ) L(\theta) L ( θ ) がアンダーフローするほど小さくなるのを避けられることです。
例 3.2 (コイン投げの最尤推定 )
表の出る確率が θ ∈ [ 0 , 1 ] \theta \in [0,1] θ ∈ [ 0 , 1 ] のコインを n n n 回投げ、k k k 回表が出たとします。各回が独立なので
L ( θ ) = θ k ( 1 − θ ) n − k , ℓ ( θ ) = k log θ + ( n − k ) log ( 1 − θ ) ( 0 < θ < 1 ) . L(\theta) = \theta^{k} (1-\theta)^{n-k}, \qquad
\ell(\theta) = k \log \theta + (n-k)\log(1-\theta) \quad (0 < \theta < 1). L ( θ ) = θ k ( 1 − θ ) n − k , ℓ ( θ ) = k log θ + ( n − k ) log ( 1 − θ ) ( 0 < θ < 1 ) . 微分して
ℓ ′ ( θ ) = k θ − n − k 1 − θ = k ( 1 − θ ) − ( n − k ) θ θ ( 1 − θ ) = k − n θ θ ( 1 − θ ) . \ell'(\theta) = \frac{k}{\theta} - \frac{n-k}{1-\theta}
= \frac{k(1-\theta) - (n-k)\theta}{\theta(1-\theta)}
= \frac{k - n\theta}{\theta(1-\theta)}. ℓ ′ ( θ ) = θ k − 1 − θ n − k = θ ( 1 − θ ) k ( 1 − θ ) − ( n − k ) θ = θ ( 1 − θ ) k − n θ . 分母は 0 < θ < 1 0 < \theta < 1 0 < θ < 1 で正なので、ℓ ′ ( θ ) = 0 ⟺ θ = k / n \ell'(\theta) = 0 \iff \theta = k/n ℓ ′ ( θ ) = 0 ⟺ θ = k / n です。さらに 0 < k < n 0 < k < n 0 < k < n のとき
ℓ ′ ′ ( θ ) = − k θ 2 − n − k ( 1 − θ ) 2 < 0 \ell''(\theta) = -\frac{k}{\theta^{2}} - \frac{n-k}{(1-\theta)^{2}} < 0 ℓ ′′ ( θ ) = − θ 2 k − ( 1 − θ ) 2 n − k < 0 なので ℓ \ell ℓ は狭義凹であり、θ = k / n \theta = k/n θ = k / n が唯一の最大点です。よって θ ^ M L = k / n \hat{\theta}_{\mathrm{ML}} = k/n θ ^ ML = k / n 、つまり単なる標本比率です。
端の場合 を省略しないでおきます。k = n k = n k = n のとき ℓ ( θ ) = n log θ \ell(\theta) = n\log\theta ℓ ( θ ) = n log θ は ( 0 , 1 ) (0,1) ( 0 , 1 ) で狭義単調増加なので、最大は端点 θ = 1 \theta = 1 θ = 1 で達成されます(L ( 1 ) = 1 L(1) = 1 L ( 1 ) = 1 )。同様に k = 0 k = 0 k = 0 なら θ ^ M L = 0 \hat{\theta}_{\mathrm{ML}} = 0 θ ^ ML = 0 です。つまり3 回投げて 3 回表なら、最尤推定は「このコインは絶対に裏が出ない」と断言します 。データと矛盾はしていませんが、賭けの根拠にはできません。この病理が次節の動機です。
次の命題は、線形回帰 で天下り的に採用した「残差の 2 乗和」が、実はガウスノイズの仮定と同じものであることを示します。
命題 3.3 (ガウス雑音の下で最尤推定は最小二乗法 ) ✓
入力 x 1 , … , x N ∈ R d \boldsymbol{x}_1,\ldots,\boldsymbol{x}_N \in \mathbb{R}^{d} x 1 , … , x N ∈ R d は固定された既知の値とし、特徴写像 ϕ : R d → R M \boldsymbol{\phi} : \mathbb{R}^{d} \to \mathbb{R}^{M} ϕ : R d → R M を用いて計画行列 Φ ∈ R N × M \Phi \in \mathbb{R}^{N \times M} Φ ∈ R N × M を、その第 i i i 行が ϕ ( x i ) T \boldsymbol{\phi}(\boldsymbol{x}_i)^{\mathsf{T}} ϕ ( x i ) T であるものとして定める。既知の定数 σ 2 > 0 \sigma^2 > 0 σ 2 > 0 に対し、出力が
y i = w T ϕ ( x i ) + ε i , ε 1 , … , ε N は独立で ε i ∼ N ( 0 , σ 2 ) y_i = \boldsymbol{w}^{\mathsf{T}} \boldsymbol{\phi}(\boldsymbol{x}_i) + \varepsilon_i, \qquad
\varepsilon_1,\ldots,\varepsilon_N \ \text{は独立で} \ \varepsilon_i \sim \mathcal{N}(0, \sigma^2) y i = w T ϕ ( x i ) + ε i , ε 1 , … , ε N は独立で ε i ∼ N ( 0 , σ 2 ) に従って生成されるとする。y = ( y 1 , … , y N ) T \boldsymbol{y} = (y_1,\ldots,y_N)^{\mathsf{T}} y = ( y 1 , … , y N ) T とおくと、w ∈ R M \boldsymbol{w} \in \mathbb{R}^{M} w ∈ R M の最尤推定量は
w ^ M L = arg min w ∈ R M ∥ y − Φ w ∥ 2 \hat{\boldsymbol{w}}_{\mathrm{ML}} = \operatorname*{arg\,min}_{\boldsymbol{w} \in \mathbb{R}^{M}} \lVert \boldsymbol{y} - \Phi \boldsymbol{w} \rVert^{2} w ^ ML = w ∈ R M arg min ∥ y − Φ w ∥ 2 の解全体と一致する。とくに Φ \Phi Φ が列フルランクならば w ^ M L = ( Φ T Φ ) − 1 Φ T y \hat{\boldsymbol{w}}_{\mathrm{ML}} = (\Phi^{\mathsf{T}}\Phi)^{-1}\Phi^{\mathsf{T}}\boldsymbol{y} w ^ ML = ( Φ T Φ ) − 1 Φ T y でただ 1 つに定まる。
証明(命題 3.3) ε i ∼ N ( 0 , σ 2 ) \varepsilon_i \sim \mathcal{N}(0,\sigma^2) ε i ∼ N ( 0 , σ 2 ) より、w \boldsymbol{w} w を与えたときの y i y_i y i の密度は
p ( y i ∣ w ) = 1 2 π σ 2 exp ( − ( y i − w T ϕ ( x i ) ) 2 2 σ 2 ) p(y_i \mid \boldsymbol{w}) = \frac{1}{\sqrt{2\pi\sigma^{2}}} \exp\!\left( -\frac{\bigl(y_i - \boldsymbol{w}^{\mathsf{T}}\boldsymbol{\phi}(\boldsymbol{x}_i)\bigr)^{2}}{2\sigma^{2}} \right) p ( y i ∣ w ) = 2 π σ 2 1 exp ( − 2 σ 2 ( y i − w T ϕ ( x i ) ) 2 ) です。ε i \varepsilon_i ε i が独立なので y 1 , … , y N y_1,\ldots,y_N y 1 , … , y N も(w \boldsymbol{w} w を与えたとき)独立で、定義 2.1 より尤度は積になります。対数を取ると
ℓ ( w ) = ∑ i = 1 N log p ( y i ∣ w ) = − N 2 log ( 2 π σ 2 ) − 1 2 σ 2 ∑ i = 1 N ( y i − w T ϕ ( x i ) ) 2 . \ell(\boldsymbol{w}) = \sum_{i=1}^{N} \log p(y_i \mid \boldsymbol{w})
= -\frac{N}{2}\log(2\pi\sigma^{2}) \;-\; \frac{1}{2\sigma^{2}} \sum_{i=1}^{N} \bigl(y_i - \boldsymbol{w}^{\mathsf{T}}\boldsymbol{\phi}(\boldsymbol{x}_i)\bigr)^{2}. ℓ ( w ) = i = 1 ∑ N log p ( y i ∣ w ) = − 2 N log ( 2 π σ 2 ) − 2 σ 2 1 i = 1 ∑ N ( y i − w T ϕ ( x i ) ) 2 . 右辺第 1 項は w \boldsymbol{w} w を含まない定数です。また ∑ i ( y i − w T ϕ ( x i ) ) 2 = ∥ y − Φ w ∥ 2 \sum_i (y_i - \boldsymbol{w}^{\mathsf{T}}\boldsymbol{\phi}(\boldsymbol{x}_i))^{2} = \lVert \boldsymbol{y} - \Phi\boldsymbol{w}\rVert^{2} ∑ i ( y i − w T ϕ ( x i ) ) 2 = ∥ y − Φ w ∥ 2 は、計画行列の定義から Φ w \Phi\boldsymbol{w} Φ w の第 i i i 成分が ϕ ( x i ) T w \boldsymbol{\phi}(\boldsymbol{x}_i)^{\mathsf{T}}\boldsymbol{w} ϕ ( x i ) T w であることによります。したがって
ℓ ( w ) = c o n s t − 1 2 σ 2 ∥ y − Φ w ∥ 2 . \ell(\boldsymbol{w}) = \mathrm{const} - \frac{1}{2\sigma^{2}} \lVert \boldsymbol{y} - \Phi\boldsymbol{w}\rVert^{2}. ℓ ( w ) = const − 2 σ 2 1 ∥ y − Φ w ∥ 2 . 係数 − 1 / ( 2 σ 2 ) -1/(2\sigma^{2}) − 1/ ( 2 σ 2 ) は負の定数なので、ℓ \ell ℓ を最大化することと ∥ y − Φ w ∥ 2 \lVert \boldsymbol{y} - \Phi\boldsymbol{w}\rVert^{2} ∥ y − Φ w ∥ 2 を最小化することは同値です。
後半は正規方程式です。J ( w ) = ∥ y − Φ w ∥ 2 J(\boldsymbol{w}) = \lVert \boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} J ( w ) = ∥ y − Φ w ∥ 2 を展開すると J = y T y − 2 y T Φ w + w T Φ T Φ w J = \boldsymbol{y}^{\mathsf{T}}\boldsymbol{y} - 2\boldsymbol{y}^{\mathsf{T}}\Phi\boldsymbol{w} + \boldsymbol{w}^{\mathsf{T}}\Phi^{\mathsf{T}}\Phi\boldsymbol{w} J = y T y − 2 y T Φ w + w T Φ T Φ w で、勾配は ∇ J = − 2 Φ T y + 2 Φ T Φ w \nabla J = -2\Phi^{\mathsf{T}}\boldsymbol{y} + 2\Phi^{\mathsf{T}}\Phi\boldsymbol{w} ∇ J = − 2 Φ T y + 2 Φ T Φ w です。J J J は凸(ヘッセ行列 2 Φ T Φ 2\Phi^{\mathsf{T}}\Phi 2 Φ T Φ が半正定値)なので ∇ J = 0 \nabla J = \boldsymbol{0} ∇ J = 0 が最小の必要十分条件で、Φ T Φ w = Φ T y \Phi^{\mathsf{T}}\Phi \boldsymbol{w} = \Phi^{\mathsf{T}}\boldsymbol{y} Φ T Φ w = Φ T y を得ます。Φ \Phi Φ が列フルランクなら Φ T Φ \Phi^{\mathsf{T}}\Phi Φ T Φ は正定値(v ≠ 0 \boldsymbol{v} \ne \boldsymbol{0} v = 0 に対し v T Φ T Φ v = ∥ Φ v ∥ 2 > 0 \boldsymbol{v}^{\mathsf{T}}\Phi^{\mathsf{T}}\Phi\boldsymbol{v} = \lVert \Phi\boldsymbol{v}\rVert^{2} > 0 v T Φ T Φ v = ∥ Φ v ∥ 2 > 0 )で可逆なので、解は一意に定まります。
∎
同じ論法で、ノイズをラプラス分布 p ( ε ) ∝ e − ∣ ε ∣ / c p(\varepsilon) \propto e^{-|\varepsilon|/c} p ( ε ) ∝ e − ∣ ε ∣/ c とすれば最尤推定は残差の絶対値和 の最小化(最小絶対偏差回帰)になります。「なぜ 2 乗か」の答えは「ガウスを仮定したから」であり、逆にいえば外れ値が多いデータでは 2 乗以外を選ぶ理由がある、ということです。
定義 4.1 (MAP 推定量 ) ↗
事前分布 p ( θ ) p(\theta) p ( θ ) を持つモデルにおいて、事後分布 p ( θ ∣ D ) p(\theta \mid \mathcal{D}) p ( θ ∣ D ) を最大にする点
θ ^ M A P = arg max θ ∈ Θ p ( θ ∣ D ) \hat{\theta}_{\mathrm{MAP}} = \operatorname*{arg\,max}_{\theta \in \Theta} p(\theta \mid \mathcal{D}) θ ^ MAP = θ ∈ Θ arg max p ( θ ∣ D ) を最大事後確率推定量 (MAP 推定量)といいます。
命題 4.2 (MAP の目的関数 ) ✓
p ( D ) > 0 p(\mathcal{D}) > 0 p ( D ) > 0 のとき、
θ ^ M A P = arg max θ ∈ Θ { log p ( D ∣ θ ) ⏟ 対数尤度 + log p ( θ ) ⏟ 対数事前分布 } \hat{\theta}_{\mathrm{MAP}} = \operatorname*{arg\,max}_{\theta \in \Theta} \Bigl\{ \underbrace{\log p(\mathcal{D} \mid \theta)}_{\text{対数尤度}} + \underbrace{\log p(\theta)}_{\text{対数事前分布}} \Bigr\} θ ^ MAP = θ ∈ Θ arg max { 対数尤度 log p ( D ∣ θ ) + 対数事前分布 log p ( θ ) } が成り立つ。とくに事前分布が Θ \Theta Θ 上で定数(一様)ならば θ ^ M A P = θ ^ M L \hat{\theta}_{\mathrm{MAP}} = \hat{\theta}_{\mathrm{ML}} θ ^ MAP = θ ^ ML である。
証明(命題 4.2) 定理 2.2 より p ( θ ∣ D ) = p ( D ∣ θ ) p ( θ ) / p ( D ) p(\theta\mid\mathcal{D}) = p(\mathcal{D}\mid\theta)p(\theta)/p(\mathcal{D}) p ( θ ∣ D ) = p ( D ∣ θ ) p ( θ ) / p ( D ) です。分母 p ( D ) p(\mathcal{D}) p ( D ) は θ \theta θ に依存しない正の定数なので、p ( θ ∣ D ) p(\theta\mid\mathcal{D}) p ( θ ∣ D ) を最大にする θ \theta θ と p ( D ∣ θ ) p ( θ ) p(\mathcal{D}\mid\theta)p(\theta) p ( D ∣ θ ) p ( θ ) を最大にする θ \theta θ は同じです。対数は狭義単調増加なので、さらに log ( p ( D ∣ θ ) p ( θ ) ) = log p ( D ∣ θ ) + log p ( θ ) \log\bigl(p(\mathcal{D}\mid\theta)p(\theta)\bigr) = \log p(\mathcal{D}\mid\theta) + \log p(\theta) log ( p ( D ∣ θ ) p ( θ ) ) = log p ( D ∣ θ ) + log p ( θ ) を最大にする θ \theta θ とも同じです。
事前分布が定数 c > 0 c > 0 c > 0 なら log p ( θ ) = log c \log p(\theta) = \log c log p ( θ ) = log c も定数なので、目的関数は対数尤度と定数だけ違い、最大点は 定義 3.1 の最尤推定量に一致します。
∎
この命題が、この章でいちばん使う道具です。MAP 推定は「対数尤度 + + + 対数事前分布」の最大化 であり、機械学習の言葉に翻訳すれば「損失関数 + + + 罰則項」の最小化にほかなりません。第 5 節でこの翻訳を厳密に実行します。
定義 4.3 (共役事前分布 ) ↗
尤度の族 { p ( D ∣ θ ) } θ ∈ Θ \{p(\mathcal{D}\mid\theta)\}_{\theta\in\Theta} { p ( D ∣ θ ) } θ ∈ Θ に対し、事前分布の族 P \mathcal{P} P が共役 であるとは、任意の p ( θ ) ∈ P p(\theta) \in \mathcal{P} p ( θ ) ∈ P と任意の観測 D \mathcal{D} D に対して事後分布 p ( θ ∣ D ) p(\theta\mid\mathcal{D}) p ( θ ∣ D ) もまた P \mathcal{P} P に属することをいいます。
共役性は数学的な必然ではなく、計算の便宜です。事後分布が同じ族に留まってくれれば、積分を実行せずに「パラメータの更新則」だけで事後分布が書けます。
次の例では、コインの表の確率 θ \theta θ の事前分布としてベータ分布を使います。ベータ分布は区間 ( 0 , 1 ) (0,1) ( 0 , 1 ) 上の分布で、密度が θ a − 1 ( 1 − θ ) b − 1 \theta^{a-1}(1-\theta)^{b-1} θ a − 1 ( 1 − θ ) b − 1 に比例するものです(a , b a,b a , b で山の位置や幅が変わります)。密度は全体の積分が 1 1 1 になるよう定数を決めれば形から一意に定まるので、「θ \theta θ に依存する部分の形が同じなら同じ分布」と読み取れます。この例のねらいは、そう読み取ることで積分を一度も計算せずに事後分布が求まる、と確かめることです。
例 4.4 (ベータ事前分布とコイン投げ )
例 3.2 と同じ設定で、事前分布としてベータ分布 B e t a ( a , b ) \mathrm{Beta}(a,b) Beta ( a , b ) (a , b > 0 a,b > 0 a , b > 0 )
p ( θ ) = 1 B ( a , b ) θ a − 1 ( 1 − θ ) b − 1 , 0 < θ < 1 p(\theta) = \frac{1}{B(a,b)}\, \theta^{a-1}(1-\theta)^{b-1}, \qquad 0 < \theta < 1 p ( θ ) = B ( a , b ) 1 θ a − 1 ( 1 − θ ) b − 1 , 0 < θ < 1 を取ります。B ( a , b ) B(a,b) B ( a , b ) は θ \theta θ に依存しない正規化定数です。定理 2.2 の比例形を使うと
p ( θ ∣ D ) ∝ θ k ( 1 − θ ) n − k ⏟ 尤度 ⋅ θ a − 1 ( 1 − θ ) b − 1 ⏟ 事前分布 = θ ( k + a ) − 1 ( 1 − θ ) ( n − k + b ) − 1 . p(\theta \mid \mathcal{D}) \;\propto\; \underbrace{\theta^{k}(1-\theta)^{n-k}}_{\text{尤度}} \cdot \underbrace{\theta^{a-1}(1-\theta)^{b-1}}_{\text{事前分布}}
= \theta^{(k+a)-1} (1-\theta)^{(n-k+b)-1}. p ( θ ∣ D ) ∝ 尤度 θ k ( 1 − θ ) n − k ⋅ 事前分布 θ a − 1 ( 1 − θ ) b − 1 = θ ( k + a ) − 1 ( 1 − θ ) ( n − k + b ) − 1 . 右辺は B e t a ( k + a , n − k + b ) \mathrm{Beta}(k+a,\, n-k+b) Beta ( k + a , n − k + b ) の密度の θ \theta θ 依存部分そのものです。密度は正規化定数まで込めて一意なので、事後分布は B e t a ( k + a , n − k + b ) \mathrm{Beta}(k+a,\, n-k+b) Beta ( k + a , n − k + b ) です。つまりベータ分布族はベルヌーイ/二項尤度の共役事前分布であり、更新則は「a a a に表の回数を、b b b に裏の回数を足す」だけです。a , b a, b a , b は疑似観測回数 と読めます。
具体的な数値 を最後まで追います。n = 3 n = 3 n = 3 、k = 3 k = 3 k = 3 (3 回投げて 3 回とも表)、事前分布は B e t a ( 2 , 2 ) \mathrm{Beta}(2,2) Beta ( 2 , 2 ) (「0.5 0.5 0.5 のあたりが怪しい」という穏やかな信念、表裏 2 回ずつを見たのと同じ重み)とします。事後分布は B e t a ( 5 , 2 ) \mathrm{Beta}(5, 2) Beta ( 5 , 2 ) で、密度は
p ( θ ∣ D ) = θ 4 ( 1 − θ ) B ( 5 , 2 ) , B ( 5 , 2 ) = Γ ( 5 ) Γ ( 2 ) Γ ( 7 ) = 4 ! ⋅ 1 ! 6 ! = 24 720 = 1 30 , p(\theta\mid\mathcal{D}) = \frac{\theta^{4}(1-\theta)}{B(5,2)}, \qquad
B(5,2) = \frac{\Gamma(5)\Gamma(2)}{\Gamma(7)} = \frac{4! \cdot 1!}{6!} = \frac{24}{720} = \frac{1}{30}, p ( θ ∣ D ) = B ( 5 , 2 ) θ 4 ( 1 − θ ) , B ( 5 , 2 ) = Γ ( 7 ) Γ ( 5 ) Γ ( 2 ) = 6 ! 4 ! ⋅ 1 ! = 720 24 = 30 1 , すなわち p ( θ ∣ D ) = 30 θ 4 ( 1 − θ ) p(\theta\mid\mathcal{D}) = 30\,\theta^{4}(1-\theta) p ( θ ∣ D ) = 30 θ 4 ( 1 − θ ) です。最頻値は
d d θ [ 4 log θ + log ( 1 − θ ) ] = 4 θ − 1 1 − θ = 4 − 5 θ θ ( 1 − θ ) = 0 ⟺ θ = 4 5 , \frac{d}{d\theta}\bigl[4\log\theta + \log(1-\theta)\bigr] = \frac{4}{\theta} - \frac{1}{1-\theta} = \frac{4 - 5\theta}{\theta(1-\theta)} = 0
\iff \theta = \frac{4}{5}, d θ d [ 4 log θ + log ( 1 − θ ) ] = θ 4 − 1 − θ 1 = θ ( 1 − θ ) 4 − 5 θ = 0 ⟺ θ = 5 4 , 二階微分は − 4 / θ 2 − 1 / ( 1 − θ ) 2 < 0 -4/\theta^{2} - 1/(1-\theta)^{2} < 0 − 4/ θ 2 − 1/ ( 1 − θ ) 2 < 0 なのでこれが唯一の最大点です。よって θ ^ M A P = 4 / 5 = 0.8 \hat{\theta}_{\mathrm{MAP}} = 4/5 = 0.8 θ ^ MAP = 4/5 = 0.8 。一方、事後平均は B e t a ( α , β ) \mathrm{Beta}(\alpha,\beta) Beta ( α , β ) の平均公式 α / ( α + β ) \alpha/(\alpha+\beta) α / ( α + β ) より 5 / 7 ≈ 0.714 5/7 \approx 0.714 5/7 ≈ 0.714 です。
まとめると、同じデータに対して
θ ^ M L = 1 , θ ^ M A P = 0.8 , E [ θ ∣ D ] = 5 7 ≈ 0.714 \hat{\theta}_{\mathrm{ML}} = 1, \qquad \hat{\theta}_{\mathrm{MAP}} = 0.8, \qquad \mathbb{E}[\theta \mid \mathcal{D}] = \frac{5}{7} \approx 0.714 θ ^ ML = 1 , θ ^ MAP = 0.8 , E [ θ ∣ D ] = 7 5 ≈ 0.714 という 3 つの答えが出ます。最尤推定の「絶対に裏は出ない」という断言は消え、しかも MAP と事後平均も一致しません。3 つ目の値は「次の 1 回が表である確率」でもあります(第 6 節の事後予測分布)。
0 0.5 1 MAP 0.8 MLE 1 事後平均 5/7 事前分布 Beta(2,2) 尤度(規格化) 事後分布 Beta(5,2) 3 回投げて 3 回表だったときの事前分布 Beta(2,2)、尤度、事後分布 Beta(5,2)。尤度は θ の密度ではないので高さは見やすいように規格化してあります。
図で見ると、事後分布が「事前分布と尤度の綱引きの結果」であることがよくわかります。尤度は θ = 1 \theta = 1 θ = 1 に向かって単調に増えていますが、事前分布が θ = 1 \theta = 1 θ = 1 で 0 0 0 に落ちるので、積は 0.8 0.8 0.8 付近で山を作ります。データが増えれば尤度の山が鋭くなり、事前分布の影響は相対的に小さくなっていきます。
次の注意のねらいは、MAP が「パラメータの名前の付け方」で変わってしまうことを、同じ事後分布を使った具体的な計算で確かめることです。θ \theta θ を θ 2 \theta^{2} θ 2 に取り替えるだけで、山の頂上の位置が動きます。
ここが本章の中心です。命題 4.2 の「対数尤度 + + + 対数事前分布」に、命題 3.3 のガウス雑音モデルとガウス事前分布を代入するだけで、リッジ回帰が出てきます。
次の定理の考え方は一言でいえば、「対数をとると、ガウス分布の指数部分がそのまま 2 乗の罰則になる」ということです。ガウス分布の密度は exp ( − ( 2 乗 ) ) \exp(-(\text{2 乗})) exp ( − ( 2 乗 )) の形なので、その対数は 2 乗の項になり、MAP の目的関数「対数尤度 + + + 対数事前分布」に足すと、損失に罰則を加えたものに見えます。
定理 5.1 (L2 正則化はガウス事前分布の MAP 推定 ) ✓
命題 3.3 と同じ設定(入力は固定、Φ ∈ R N × M \Phi \in \mathbb{R}^{N\times M} Φ ∈ R N × M は計画行列、σ 2 > 0 \sigma^2 > 0 σ 2 > 0 は既知の雑音分散、ε i \varepsilon_i ε i は独立に N ( 0 , σ 2 ) \mathcal{N}(0,\sigma^2) N ( 0 , σ 2 ) )に加えて、重みの事前分布を
w ∼ N ( 0 , τ 2 I M ) , τ 2 > 0 は既知 \boldsymbol{w} \sim \mathcal{N}(\boldsymbol{0},\, \tau^{2} I_M), \qquad \tau^{2} > 0 \ \text{は既知} w ∼ N ( 0 , τ 2 I M ) , τ 2 > 0 は既知 とし、w \boldsymbol{w} w と ( ε 1 , … , ε N ) (\varepsilon_1,\ldots,\varepsilon_N) ( ε 1 , … , ε N ) は独立とする。このとき MAP 推定量は
w ^ M A P = arg min w ∈ R M { ∥ y − Φ w ∥ 2 + λ ∥ w ∥ 2 } , λ = σ 2 τ 2 \hat{\boldsymbol{w}}_{\mathrm{MAP}}
= \operatorname*{arg\,min}_{\boldsymbol{w} \in \mathbb{R}^{M}}
\Bigl\{ \lVert \boldsymbol{y} - \Phi\boldsymbol{w}\rVert^{2} + \lambda \lVert \boldsymbol{w}\rVert^{2} \Bigr\},
\qquad \lambda = \frac{\sigma^{2}}{\tau^{2}} w ^ MAP = w ∈ R M arg min { ∥ y − Φ w ∥ 2 + λ ∥ w ∥ 2 } , λ = τ 2 σ 2 で与えられ、これはリッジ回帰(L2 正則化つき最小二乗法)の解に一致する。さらにこの最小化問題の解は一意で、
w ^ M A P = ( Φ T Φ + λ I M ) − 1 Φ T y \hat{\boldsymbol{w}}_{\mathrm{MAP}} = \bigl(\Phi^{\mathsf{T}}\Phi + \lambda I_M\bigr)^{-1} \Phi^{\mathsf{T}} \boldsymbol{y} w ^ MAP = ( Φ T Φ + λ I M ) − 1 Φ T y である。
証明(定理 5.1) 第 1 段:目的関数を書き下す。 命題 3.3 の証明で計算したとおり、対数尤度は
log p ( y ∣ w ) = − N 2 log ( 2 π σ 2 ) − 1 2 σ 2 ∥ y − Φ w ∥ 2 \log p(\boldsymbol{y}\mid\boldsymbol{w}) = -\frac{N}{2}\log(2\pi\sigma^{2}) - \frac{1}{2\sigma^{2}}\lVert \boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} log p ( y ∣ w ) = − 2 N log ( 2 π σ 2 ) − 2 σ 2 1 ∥ y − Φ w ∥ 2 です(w \boldsymbol{w} w と ε i \varepsilon_i ε i の独立性から、w \boldsymbol{w} w を与えたときの y \boldsymbol{y} y の条件付き分布は元のモデルのままであることを使いました)。また M M M 次元ガウス事前分布 N ( 0 , τ 2 I M ) \mathcal{N}(\boldsymbol{0},\tau^{2}I_M) N ( 0 , τ 2 I M ) の密度は p ( w ) = ( 2 π τ 2 ) − M / 2 exp ( − ∥ w ∥ 2 / ( 2 τ 2 ) ) p(\boldsymbol{w}) = (2\pi\tau^{2})^{-M/2}\exp\bigl(-\lVert\boldsymbol{w}\rVert^{2}/(2\tau^{2})\bigr) p ( w ) = ( 2 π τ 2 ) − M /2 exp ( − ∥ w ∥ 2 / ( 2 τ 2 ) ) なので
log p ( w ) = − M 2 log ( 2 π τ 2 ) − 1 2 τ 2 ∥ w ∥ 2 . \log p(\boldsymbol{w}) = -\frac{M}{2}\log(2\pi\tau^{2}) - \frac{1}{2\tau^{2}}\lVert \boldsymbol{w}\rVert^{2}. log p ( w ) = − 2 M log ( 2 π τ 2 ) − 2 τ 2 1 ∥ w ∥ 2 . 第 2 段:定数を落として整理する。 命題 4.2 より w ^ M A P \hat{\boldsymbol{w}}_{\mathrm{MAP}} w ^ MAP は
log p ( y ∣ w ) + log p ( w ) = c o n s t − 1 2 σ 2 ∥ y − Φ w ∥ 2 − 1 2 τ 2 ∥ w ∥ 2 \log p(\boldsymbol{y}\mid\boldsymbol{w}) + \log p(\boldsymbol{w})
= \mathrm{const} - \frac{1}{2\sigma^{2}}\lVert \boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} - \frac{1}{2\tau^{2}}\lVert \boldsymbol{w}\rVert^{2} log p ( y ∣ w ) + log p ( w ) = const − 2 σ 2 1 ∥ y − Φ w ∥ 2 − 2 τ 2 1 ∥ w ∥ 2 の最大点です。c o n s t \mathrm{const} const は w \boldsymbol{w} w を含まないので落とせます。符号を反転すると最小化問題
w ^ M A P = arg min w { 1 2 σ 2 ∥ y − Φ w ∥ 2 + 1 2 τ 2 ∥ w ∥ 2 } \hat{\boldsymbol{w}}_{\mathrm{MAP}} = \operatorname*{arg\,min}_{\boldsymbol{w}}
\left\{ \frac{1}{2\sigma^{2}}\lVert \boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} + \frac{1}{2\tau^{2}}\lVert \boldsymbol{w}\rVert^{2} \right\} w ^ MAP = w arg min { 2 σ 2 1 ∥ y − Φ w ∥ 2 + 2 τ 2 1 ∥ w ∥ 2 } になります。目的関数全体を正の定数 2 σ 2 2\sigma^{2} 2 σ 2 倍しても最小点は変わらないので
w ^ M A P = arg min w { ∥ y − Φ w ∥ 2 + σ 2 τ 2 ∥ w ∥ 2 } \hat{\boldsymbol{w}}_{\mathrm{MAP}} = \operatorname*{arg\,min}_{\boldsymbol{w}}
\left\{ \lVert \boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} + \frac{\sigma^{2}}{\tau^{2}} \lVert \boldsymbol{w}\rVert^{2} \right\} w ^ MAP = w arg min { ∥ y − Φ w ∥ 2 + τ 2 σ 2 ∥ w ∥ 2 } を得ます。λ = σ 2 / τ 2 \lambda = \sigma^{2}/\tau^{2} λ = σ 2 / τ 2 とおけば主張の形です。
第 3 段:閉じた形の解。 J ( w ) = ∥ y − Φ w ∥ 2 + λ ∥ w ∥ 2 J(\boldsymbol{w}) = \lVert \boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} + \lambda\lVert\boldsymbol{w}\rVert^{2} J ( w ) = ∥ y − Φ w ∥ 2 + λ ∥ w ∥ 2 を展開すると
J ( w ) = y T y − 2 y T Φ w + w T ( Φ T Φ + λ I M ) w , J(\boldsymbol{w}) = \boldsymbol{y}^{\mathsf{T}}\boldsymbol{y} - 2\boldsymbol{y}^{\mathsf{T}}\Phi\boldsymbol{w} + \boldsymbol{w}^{\mathsf{T}}\bigl(\Phi^{\mathsf{T}}\Phi + \lambda I_M\bigr)\boldsymbol{w}, J ( w ) = y T y − 2 y T Φ w + w T ( Φ T Φ + λ I M ) w , 勾配は ∇ J ( w ) = − 2 Φ T y + 2 ( Φ T Φ + λ I M ) w \nabla J(\boldsymbol{w}) = -2\Phi^{\mathsf{T}}\boldsymbol{y} + 2(\Phi^{\mathsf{T}}\Phi+\lambda I_M)\boldsymbol{w} ∇ J ( w ) = − 2 Φ T y + 2 ( Φ T Φ + λ I M ) w です。ここで A = Φ T Φ + λ I M A = \Phi^{\mathsf{T}}\Phi + \lambda I_M A = Φ T Φ + λ I M は正定値です。実際、v ≠ 0 \boldsymbol{v}\ne\boldsymbol{0} v = 0 に対し
v T A v = ∥ Φ v ∥ 2 + λ ∥ v ∥ 2 ≥ λ ∥ v ∥ 2 > 0 \boldsymbol{v}^{\mathsf{T}} A \boldsymbol{v} = \lVert \Phi\boldsymbol{v}\rVert^{2} + \lambda\lVert\boldsymbol{v}\rVert^{2} \ge \lambda \lVert \boldsymbol{v}\rVert^{2} > 0 v T A v = ∥ Φ v ∥ 2 + λ ∥ v ∥ 2 ≥ λ ∥ v ∥ 2 > 0 となります(λ = σ 2 / τ 2 > 0 \lambda = \sigma^2/\tau^2 > 0 λ = σ 2 / τ 2 > 0 を使いました)。したがって J J J のヘッセ行列 2 A 2A 2 A は正定値で J J J は狭義凸、∇ J ( w ) = 0 \nabla J(\boldsymbol{w}) = \boldsymbol{0} ∇ J ( w ) = 0 すなわち A w = Φ T y A\boldsymbol{w} = \Phi^{\mathsf{T}}\boldsymbol{y} A w = Φ T y が唯一の最小点を与えます。A A A は正定値ゆえ可逆なので w ^ M A P = A − 1 Φ T y \hat{\boldsymbol{w}}_{\mathrm{MAP}} = A^{-1}\Phi^{\mathsf{T}}\boldsymbol{y} w ^ MAP = A − 1 Φ T y です。
∎
λ = σ 2 / τ 2 \lambda = \sigma^{2}/\tau^{2} λ = σ 2 / τ 2 という表式は、正則化係数に明確な意味を与えます。λ \lambda λ は「観測がどれだけ雑か」と「重みがどれだけ大きくてよいか」の比 です。雑音が大きい(σ 2 \sigma^2 σ 2 大)ほどデータを信用せず罰則を強め、事前の許容幅が広い(τ 2 \tau^2 τ 2 大)ほど罰則を緩めます。τ 2 → ∞ \tau^{2}\to\infty τ 2 → ∞ (何も知らない)とすれば λ → 0 \lambda\to 0 λ → 0 で最尤推定に戻り、これは 命題 4.2 の後半(一様事前分布なら MLE)とも整合します。
例 5.2 (リッジ回帰は固有値の小さい方向を強く縮める )
Φ T Φ \Phi^{\mathsf{T}}\Phi Φ T Φ は対称半正定値なので、スペクトル定理 (系 4.3[スペクトル定理] )により直交行列 U U U と対角行列 Λ = diag ( d 1 , … , d M ) \Lambda = \operatorname{diag}(d_1,\ldots,d_M) Λ = diag ( d 1 , … , d M ) (d j ≥ 0 d_j \ge 0 d j ≥ 0 )を用いて Φ T Φ = U Λ U T \Phi^{\mathsf{T}}\Phi = U\Lambda U^{\mathsf{T}} Φ T Φ = U Λ U T と書けます。U T U = I U^{\mathsf{T}}U = I U T U = I より Φ T Φ + λ I = U ( Λ + λ I ) U T \Phi^{\mathsf{T}}\Phi + \lambda I = U(\Lambda+\lambda I)U^{\mathsf{T}} Φ T Φ + λ I = U ( Λ + λ I ) U T なので、定理 5.1 の解は
w ^ M A P = U ( Λ + λ I ) − 1 U T Φ T y . \hat{\boldsymbol{w}}_{\mathrm{MAP}} = U(\Lambda+\lambda I)^{-1}U^{\mathsf{T}}\Phi^{\mathsf{T}}\boldsymbol{y}. w ^ MAP = U ( Λ + λ I ) − 1 U T Φ T y . z = U T Φ T y \boldsymbol{z} = U^{\mathsf{T}}\Phi^{\mathsf{T}}\boldsymbol{y} z = U T Φ T y とおき、固有ベクトル基底での成分を比べます。d j > 0 d_j > 0 d j > 0 のとき最尤解の第 j j j 成分は z j / d j z_j/d_j z j / d j 、リッジ解の第 j j j 成分は z j / ( d j + λ ) z_j/(d_j+\lambda) z j / ( d j + λ ) なので、その比は
( U T w ^ M A P ) j ( U T w ^ M L ) j = d j d j + λ . \frac{(U^{\mathsf{T}}\hat{\boldsymbol{w}}_{\mathrm{MAP}})_j}{(U^{\mathsf{T}}\hat{\boldsymbol{w}}_{\mathrm{ML}})_j} = \frac{d_j}{d_j + \lambda}. ( U T w ^ ML ) j ( U T w ^ MAP ) j = d j + λ d j . たとえば λ = 1 \lambda = 1 λ = 1 、d 1 = 100 d_1 = 100 d 1 = 100 、d 2 = 0.01 d_2 = 0.01 d 2 = 0.01 なら、縮小率は第 1 方向で 100 / 101 ≈ 0.990 100/101 \approx 0.990 100/101 ≈ 0.990 、第 2 方向で 0.01 / 1.01 ≈ 0.0099 0.01/1.01 \approx 0.0099 0.01/1.01 ≈ 0.0099 です。データの分散が大きい方向はほとんど手つかず、分散がほとんどない方向は 1 / 100 1/100 1/100 に潰されます 。
d j d_j d j は Φ T Φ \Phi^{\mathsf{T}}\Phi Φ T Φ の固有値、つまり主成分分析 でいうところの各主成分方向(定義 3.2[主成分分析] )のデータの散らばりです。「データが何も語っていない方向では事前分布が勝つ」という 定理 5.1 のベイズ的な読みが、そのまま数式に現れています。d j = 0 d_j = 0 d j = 0 (その方向にはデータが皆無)なら成分は 0 0 0 、すなわち事前分布の平均そのものになります。
同じ計算をラプラス事前分布で行うと L1 正則化(ラッソ)が出ます。
命題 5.3 (L1 正則化はラプラス事前分布の MAP 推定 ) ✓
定理 5.1 と同じ雑音モデルの下で、重みの事前分布を、各成分が独立に平均 0 0 0 ・尺度 b > 0 b > 0 b > 0 のラプラス分布
p ( w ) = ∏ j = 1 M 1 2 b exp ( − ∣ w j ∣ b ) = 1 ( 2 b ) M exp ( − ∥ w ∥ 1 b ) p(\boldsymbol{w}) = \prod_{j=1}^{M} \frac{1}{2b}\exp\!\left(-\frac{|w_j|}{b}\right)
= \frac{1}{(2b)^{M}} \exp\!\left(-\frac{\lVert \boldsymbol{w}\rVert_{1}}{b}\right) p ( w ) = j = 1 ∏ M 2 b 1 exp ( − b ∣ w j ∣ ) = ( 2 b ) M 1 exp ( − b ∥ w ∥ 1 ) に従うものとする(∥ w ∥ 1 = ∑ j ∣ w j ∣ \lVert\boldsymbol{w}\rVert_1 = \sum_j |w_j| ∥ w ∥ 1 = ∑ j ∣ w j ∣ )。このとき
w ^ M A P = arg min w ∈ R M { ∥ y − Φ w ∥ 2 + λ ∥ w ∥ 1 } , λ = 2 σ 2 b \hat{\boldsymbol{w}}_{\mathrm{MAP}} = \operatorname*{arg\,min}_{\boldsymbol{w}\in\mathbb{R}^{M}}
\Bigl\{ \lVert \boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} + \lambda \lVert \boldsymbol{w}\rVert_{1} \Bigr\},
\qquad \lambda = \frac{2\sigma^{2}}{b} w ^ MAP = w ∈ R M arg min { ∥ y − Φ w ∥ 2 + λ ∥ w ∥ 1 } , λ = b 2 σ 2 であり、これはラッソ(L1 正則化つき最小二乗法)の解に一致する。
証明(命題 5.3) 対数事前分布は
log p ( w ) = − M log ( 2 b ) − 1 b ∥ w ∥ 1 \log p(\boldsymbol{w}) = -M\log(2b) - \frac{1}{b}\lVert\boldsymbol{w}\rVert_{1} log p ( w ) = − M log ( 2 b ) − b 1 ∥ w ∥ 1 です。命題 4.2 と 定理 5.1 の証明第 1 段の対数尤度を合わせると、w ^ M A P \hat{\boldsymbol{w}}_{\mathrm{MAP}} w ^ MAP は
c o n s t − 1 2 σ 2 ∥ y − Φ w ∥ 2 − 1 b ∥ w ∥ 1 \mathrm{const} - \frac{1}{2\sigma^{2}}\lVert\boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} - \frac{1}{b}\lVert\boldsymbol{w}\rVert_{1} const − 2 σ 2 1 ∥ y − Φ w ∥ 2 − b 1 ∥ w ∥ 1 の最大点です。定数を落とし、符号を反転し、全体を正の定数 2 σ 2 2\sigma^{2} 2 σ 2 倍すると
w ^ M A P = arg min w { ∥ y − Φ w ∥ 2 + 2 σ 2 b ∥ w ∥ 1 } \hat{\boldsymbol{w}}_{\mathrm{MAP}} = \operatorname*{arg\,min}_{\boldsymbol{w}}
\Bigl\{ \lVert\boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} + \frac{2\sigma^{2}}{b}\lVert\boldsymbol{w}\rVert_{1} \Bigr\} w ^ MAP = w arg min { ∥ y − Φ w ∥ 2 + b 2 σ 2 ∥ w ∥ 1 } を得ます。λ = 2 σ 2 / b \lambda = 2\sigma^{2}/b λ = 2 σ 2 / b とおけば主張の形です。
なお目的関数は凸(2 乗項は凸、∥ ⋅ ∥ 1 \lVert\cdot\rVert_1 ∥ ⋅ ∥ 1 はノルムなので凸)ですが ∥ ⋅ ∥ 1 \lVert\cdot\rVert_1 ∥ ⋅ ∥ 1 が w j = 0 w_j = 0 w j = 0 で微分可能でないため、定理 5.1 のような閉じた形の解は一般には得られません。Φ T Φ = I \Phi^{\mathsf{T}}\Phi = I Φ T Φ = I (列が正規直交)の特別な場合は成分ごとに分離でき、解はソフトしきい値関数で書けます(演習 7.3 )。
∎
対応関係を表にまとめます。
罰則項 対応する事前分布 事前密度(θ \theta θ 依存部分) 正則化係数 解の特徴 なし 一様(非正則な事前分布) ∝ 1 \propto 1 ∝ 1 λ = 0 \lambda = 0 λ = 0 最尤推定に一致 λ ∥ w ∥ 2 2 \lambda \lVert \boldsymbol{w}\rVert_{2}^{2} λ ∥ w ∥ 2 2 ガウス N ( 0 , τ 2 I ) \mathcal{N}(\boldsymbol{0},\tau^{2}I) N ( 0 , τ 2 I ) exp ( − ∥ w ∥ 2 2 / ( 2 τ 2 ) ) \exp\bigl(-\lVert\boldsymbol{w}\rVert_2^{2}/(2\tau^{2})\bigr) exp ( − ∥ w ∥ 2 2 / ( 2 τ 2 ) ) λ = σ 2 / τ 2 \lambda = \sigma^{2}/\tau^{2} λ = σ 2 / τ 2 全成分を滑らかに縮小 λ ∥ w ∥ 1 \lambda \lVert \boldsymbol{w}\rVert_{1} λ ∥ w ∥ 1 ラプラス(尺度 b b b ) exp ( − ∥ w ∥ 1 / b ) \exp\bigl(-\lVert\boldsymbol{w}\rVert_1/b\bigr) exp ( − ∥ w ∥ 1 / b ) λ = 2 σ 2 / b \lambda = 2\sigma^{2}/b λ = 2 σ 2 / b 一部の成分が厳密に 0 0 0 λ ∥ w − μ ∥ 2 2 \lambda\lVert\boldsymbol{w}-\boldsymbol{\mu}\rVert_2^{2} λ ∥ w − μ ∥ 2 2 ガウス N ( μ , τ 2 I ) \mathcal{N}(\boldsymbol{\mu},\tau^{2}I) N ( μ , τ 2 I ) exp ( − ∥ w − μ ∥ 2 2 / ( 2 τ 2 ) ) \exp\bigl(-\lVert\boldsymbol{w}-\boldsymbol{\mu}\rVert_2^{2}/(2\tau^{2})\bigr) exp ( − ∥ w − μ ∥ 2 2 / ( 2 τ 2 ) ) λ = σ 2 / τ 2 \lambda = \sigma^{2}/\tau^{2} λ = σ 2 / τ 2 既知の値 μ \boldsymbol{\mu} μ に引き寄せる
4 行目は、事前学習済みモデルからの微調整で「元の重みから離れすぎない」ようにする正則化が、そのまま「元の重みを平均とするガウス事前分布」であることを示しています。正則化を発明するとは、事前分布を設計することです。
-3 0 3 w ラプラス(尖った山) ガウス(滑らかな山) 分散を 1 にそろえたガウス分布(実線)とラプラス分布(尺度 b = 1/√2、破線)の密度。ラプラス分布は原点で尖り、裾は重い
MLE も MAP も、最後に θ \theta θ を 1 点に潰します。潰さずに事後分布のまま持ち歩くとどうなるかを見ます。線形回帰+ガウス事前分布は、この計算が手で最後まで実行できる数少ない例です。
次の定理の証明の方針は 3 段です。(1) ベイズの定理の比例形で「尤度 × \times × 事前分布」を書く。(2) 指数の中を w \boldsymbol{w} w の 2 次式に整理する(長い行列の展開ですが、w \boldsymbol{w} w について 2 次の項と 1 次の項を集めるだけの作業です)。(3) 付録の補題で、その形の密度がどのガウス分布かを読み取る。本質は (3) で、2 次式の形をしていればガウス分布だ、という点です。
定理 6.1 (ベイズ線形回帰の事後分布 ) ✓
命題 3.3 の雑音モデル(入力は固定、Φ ∈ R N × M \Phi\in\mathbb{R}^{N\times M} Φ ∈ R N × M 、σ 2 > 0 \sigma^{2}>0 σ 2 > 0 は既知)の下で、w \boldsymbol{w} w の事前分布を N ( m 0 , S 0 ) \mathcal{N}(\boldsymbol{m}_0, S_0) N ( m 0 , S 0 ) (S 0 S_0 S 0 は M M M 次の対称正定値行列、m 0 ∈ R M \boldsymbol{m}_0\in\mathbb{R}^M m 0 ∈ R M )とし、w \boldsymbol{w} w と雑音は独立とする。このとき事後分布は再びガウス分布
p ( w ∣ y ) = N ( w ∣ m N , S N ) p(\boldsymbol{w}\mid\boldsymbol{y}) = \mathcal{N}(\boldsymbol{w} \mid \boldsymbol{m}_N, S_N) p ( w ∣ y ) = N ( w ∣ m N , S N ) であり、
S N = ( S 0 − 1 + 1 σ 2 Φ T Φ ) − 1 , m N = S N ( S 0 − 1 m 0 + 1 σ 2 Φ T y ) S_N = \left( S_0^{-1} + \frac{1}{\sigma^{2}}\Phi^{\mathsf{T}}\Phi \right)^{-1},
\qquad
\boldsymbol{m}_N = S_N\left( S_0^{-1}\boldsymbol{m}_0 + \frac{1}{\sigma^{2}}\Phi^{\mathsf{T}}\boldsymbol{y} \right) S N = ( S 0 − 1 + σ 2 1 Φ T Φ ) − 1 , m N = S N ( S 0 − 1 m 0 + σ 2 1 Φ T y ) で与えられる。とくにガウス分布は共役事前分布である(定義 4.3 )。
証明(定理 6.1) 第 1 段:逆行列の存在。 A = S 0 − 1 + σ − 2 Φ T Φ A = S_0^{-1} + \sigma^{-2}\Phi^{\mathsf{T}}\Phi A = S 0 − 1 + σ − 2 Φ T Φ とおきます。S 0 S_0 S 0 が対称正定値なら S 0 − 1 S_0^{-1} S 0 − 1 も対称正定値です。また任意の v \boldsymbol{v} v に対し v T Φ T Φ v = ∥ Φ v ∥ 2 ≥ 0 \boldsymbol{v}^{\mathsf{T}}\Phi^{\mathsf{T}}\Phi\boldsymbol{v} = \lVert\Phi\boldsymbol{v}\rVert^{2}\ge 0 v T Φ T Φ v = ∥ Φ v ∥ 2 ≥ 0 なので Φ T Φ \Phi^{\mathsf{T}}\Phi Φ T Φ は半正定値です。よって v ≠ 0 \boldsymbol{v}\ne\boldsymbol{0} v = 0 のとき
v T A v = v T S 0 − 1 v + 1 σ 2 ∥ Φ v ∥ 2 > 0 \boldsymbol{v}^{\mathsf{T}}A\boldsymbol{v} = \boldsymbol{v}^{\mathsf{T}}S_0^{-1}\boldsymbol{v} + \frac{1}{\sigma^{2}}\lVert\Phi\boldsymbol{v}\rVert^{2} > 0 v T A v = v T S 0 − 1 v + σ 2 1 ∥ Φ v ∥ 2 > 0 であり A A A は正定値、したがって可逆です。S N = A − 1 S_N = A^{-1} S N = A − 1 と書きます。
第 2 段:指数部を w \boldsymbol{w} w の 2 次式に整理する。 定理 2.2 の比例形より
p ( w ∣ y ) ∝ p ( y ∣ w ) p ( w ) ∝ exp ( − 1 2 σ 2 ∥ y − Φ w ∥ 2 − 1 2 ( w − m 0 ) T S 0 − 1 ( w − m 0 ) ) p(\boldsymbol{w}\mid\boldsymbol{y}) \propto p(\boldsymbol{y}\mid\boldsymbol{w})\,p(\boldsymbol{w})
\propto \exp\left( -\frac{1}{2\sigma^{2}}\lVert\boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} - \frac{1}{2}(\boldsymbol{w}-\boldsymbol{m}_0)^{\mathsf{T}}S_0^{-1}(\boldsymbol{w}-\boldsymbol{m}_0) \right) p ( w ∣ y ) ∝ p ( y ∣ w ) p ( w ) ∝ exp ( − 2 σ 2 1 ∥ y − Φ w ∥ 2 − 2 1 ( w − m 0 ) T S 0 − 1 ( w − m 0 ) ) です(w \boldsymbol{w} w を含まない正規化定数はすべて比例記号に吸収しました)。指数の中を展開します。
∥ y − Φ w ∥ 2 = y T y − 2 y T Φ w + w T Φ T Φ w , \lVert\boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} = \boldsymbol{y}^{\mathsf{T}}\boldsymbol{y} - 2\boldsymbol{y}^{\mathsf{T}}\Phi\boldsymbol{w} + \boldsymbol{w}^{\mathsf{T}}\Phi^{\mathsf{T}}\Phi\boldsymbol{w}, ∥ y − Φ w ∥ 2 = y T y − 2 y T Φ w + w T Φ T Φ w , ( w − m 0 ) T S 0 − 1 ( w − m 0 ) = w T S 0 − 1 w − 2 m 0 T S 0 − 1 w + m 0 T S 0 − 1 m 0 (\boldsymbol{w}-\boldsymbol{m}_0)^{\mathsf{T}}S_0^{-1}(\boldsymbol{w}-\boldsymbol{m}_0) = \boldsymbol{w}^{\mathsf{T}}S_0^{-1}\boldsymbol{w} - 2\boldsymbol{m}_0^{\mathsf{T}}S_0^{-1}\boldsymbol{w} + \boldsymbol{m}_0^{\mathsf{T}}S_0^{-1}\boldsymbol{m}_0 ( w − m 0 ) T S 0 − 1 ( w − m 0 ) = w T S 0 − 1 w − 2 m 0 T S 0 − 1 w + m 0 T S 0 − 1 m 0 (S 0 − 1 S_0^{-1} S 0 − 1 の対称性から交差項 2 つが等しいことを使いました)。w \boldsymbol{w} w を含まない項をまた比例記号に吸収すると、指数は
− 1 2 w T ( S 0 − 1 + 1 σ 2 Φ T Φ ) ⏟ = A w + ( S 0 − 1 m 0 + 1 σ 2 Φ T y ) T ⏟ = b T w + c o n s t -\frac{1}{2}\boldsymbol{w}^{\mathsf{T}}\underbrace{\left(S_0^{-1}+\frac{1}{\sigma^{2}}\Phi^{\mathsf{T}}\Phi\right)}_{= A}\boldsymbol{w}
\;+\; \underbrace{\left(S_0^{-1}\boldsymbol{m}_0 + \frac{1}{\sigma^{2}}\Phi^{\mathsf{T}}\boldsymbol{y}\right)^{\mathsf{T}}}_{=\ \boldsymbol{b}^{\mathsf{T}}}\boldsymbol{w}
\;+\; \mathrm{const} − 2 1 w T = A ( S 0 − 1 + σ 2 1 Φ T Φ ) w + = b T ( S 0 − 1 m 0 + σ 2 1 Φ T y ) T w + const の形になります。
第 3 段:ガウス分布と同定する。 第 2 段より、事後密度はある定数 C > 0 C > 0 C > 0 を用いて p ( w ∣ y ) = C exp ( − 1 2 w T A w + b T w ) p(\boldsymbol{w}\mid\boldsymbol{y}) = C\exp\bigl(-\tfrac12\boldsymbol{w}^{\mathsf{T}}A\boldsymbol{w} + \boldsymbol{b}^{\mathsf{T}}\boldsymbol{w}\bigr) p ( w ∣ y ) = C exp ( − 2 1 w T A w + b T w ) と書けます。A A A は第 1 段より対称正定値なので、補題 7.5 が適用でき、この密度は N ( A − 1 b , A − 1 ) \mathcal{N}(A^{-1}\boldsymbol{b},\, A^{-1}) N ( A − 1 b , A − 1 ) の密度に一致します。A − 1 = S N A^{-1}=S_N A − 1 = S N 、A − 1 b = S N ( S 0 − 1 m 0 + σ − 2 Φ T y ) = m N A^{-1}\boldsymbol{b} = S_N(S_0^{-1}\boldsymbol{m}_0 + \sigma^{-2}\Phi^{\mathsf{T}}\boldsymbol{y}) = \boldsymbol{m}_N A − 1 b = S N ( S 0 − 1 m 0 + σ − 2 Φ T y ) = m N なので主張が従います。
∎
系 6.2 (ガウス事前分布の下では MAP は事後平均であり、リッジ解に一致する ) ✓
定理 6.1 の設定で m 0 = 0 \boldsymbol{m}_0 = \boldsymbol{0} m 0 = 0 、S 0 = τ 2 I M S_0 = \tau^{2}I_M S 0 = τ 2 I M とすると、
w ^ M A P = E [ w ∣ y ] = m N = ( Φ T Φ + λ I M ) − 1 Φ T y , λ = σ 2 τ 2 \hat{\boldsymbol{w}}_{\mathrm{MAP}} = \mathbb{E}[\boldsymbol{w}\mid\boldsymbol{y}] = \boldsymbol{m}_N = \bigl(\Phi^{\mathsf{T}}\Phi + \lambda I_M\bigr)^{-1}\Phi^{\mathsf{T}}\boldsymbol{y}, \qquad \lambda = \frac{\sigma^{2}}{\tau^{2}} w ^ MAP = E [ w ∣ y ] = m N = ( Φ T Φ + λ I M ) − 1 Φ T y , λ = τ 2 σ 2 が成り立つ。
証明(系 6.2) ガウス分布の密度 ∝ exp ( − 1 2 ( w − m N ) T S N − 1 ( w − m N ) ) \propto \exp\bigl(-\tfrac12(\boldsymbol{w}-\boldsymbol{m}_N)^{\mathsf{T}}S_N^{-1}(\boldsymbol{w}-\boldsymbol{m}_N)\bigr) ∝ exp ( − 2 1 ( w − m N ) T S N − 1 ( w − m N ) ) は、指数が w = m N \boldsymbol{w}=\boldsymbol{m}_N w = m N で最大値 0 0 0 を取り、S N − 1 S_N^{-1} S N − 1 が正定値なので他の点では真に負です。よって最頻値は平均 m N \boldsymbol{m}_N m N に一致します。
具体形は代入するだけです。S 0 − 1 = τ − 2 I S_0^{-1} = \tau^{-2}I S 0 − 1 = τ − 2 I より
S N = ( 1 τ 2 I + 1 σ 2 Φ T Φ ) − 1 = σ 2 ( σ 2 τ 2 I + Φ T Φ ) − 1 = σ 2 ( Φ T Φ + λ I ) − 1 , S_N = \left(\frac{1}{\tau^{2}}I + \frac{1}{\sigma^{2}}\Phi^{\mathsf{T}}\Phi\right)^{-1}
= \sigma^{2}\left(\frac{\sigma^{2}}{\tau^{2}}I + \Phi^{\mathsf{T}}\Phi\right)^{-1}
= \sigma^{2}\bigl(\Phi^{\mathsf{T}}\Phi + \lambda I\bigr)^{-1}, S N = ( τ 2 1 I + σ 2 1 Φ T Φ ) − 1 = σ 2 ( τ 2 σ 2 I + Φ T Φ ) − 1 = σ 2 ( Φ T Φ + λ I ) − 1 , m N = S N ( 0 + 1 σ 2 Φ T y ) = σ 2 ( Φ T Φ + λ I ) − 1 1 σ 2 Φ T y = ( Φ T Φ + λ I ) − 1 Φ T y . \boldsymbol{m}_N = S_N\left(\boldsymbol{0} + \frac{1}{\sigma^{2}}\Phi^{\mathsf{T}}\boldsymbol{y}\right)
= \sigma^{2}\bigl(\Phi^{\mathsf{T}}\Phi+\lambda I\bigr)^{-1}\frac{1}{\sigma^{2}}\Phi^{\mathsf{T}}\boldsymbol{y}
= \bigl(\Phi^{\mathsf{T}}\Phi+\lambda I\bigr)^{-1}\Phi^{\mathsf{T}}\boldsymbol{y}. m N = S N ( 0 + σ 2 1 Φ T y ) = σ 2 ( Φ T Φ + λ I ) − 1 σ 2 1 Φ T y = ( Φ T Φ + λ I ) − 1 Φ T y . これは 定理 5.1 で別の道筋(目的関数の最小化)から得た解と一致しています。
∎
事後共分散 S N = σ 2 ( Φ T Φ + λ I ) − 1 S_N = \sigma^{2}(\Phi^{\mathsf{T}}\Phi+\lambda I)^{-1} S N = σ 2 ( Φ T Φ + λ I ) − 1 は、点推定だけでは決して得られない情報です。これを使って予測に不確実性を持たせます。
定義 6.3 (事後予測分布 ) ↗
新しい入力 x ∗ \boldsymbol{x}_{*} x ∗ に対する出力 y ∗ y_{*} y ∗ の事後予測分布 とは、パラメータを事後分布で積分消去した分布
p ( y ∗ ∣ x ∗ , D ) = ∫ p ( y ∗ ∣ x ∗ , w ) p ( w ∣ D ) d w p(y_{*}\mid \boldsymbol{x}_{*}, \mathcal{D}) = \int p(y_{*}\mid \boldsymbol{x}_{*}, \boldsymbol{w})\, p(\boldsymbol{w}\mid\mathcal{D})\, d\boldsymbol{w} p ( y ∗ ∣ x ∗ , D ) = ∫ p ( y ∗ ∣ x ∗ , w ) p ( w ∣ D ) d w のことです。w ^ \hat{\boldsymbol{w}} w ^ を 1 つ選んで p ( y ∗ ∣ x ∗ , w ^ ) p(y_{*}\mid\boldsymbol{x}_{*},\hat{\boldsymbol{w}}) p ( y ∗ ∣ x ∗ , w ^ ) を使うプラグイン予測 と対比されます。
次の定理は、「新しい入力での予測のばらつき」を、観測ノイズの分だけの部分とパラメータが分からない分の部分に分けて数式で書くものです。証明では、ガウス分布に従う量の線形結合や和がまたガウス分布になることを使います。特性関数(確率分布を E [ e i t X ] \mathbb{E}[e^{itX}] E [ e i tX ] という関数で表したもの)は、その確認に使う道具です。
定理 6.4 (ベイズ線形回帰の事後予測分布 ) ✓
定理 6.1 の設定の下で、新しい入力 x ∗ \boldsymbol{x}_{*} x ∗ の特徴ベクトルを ϕ ∗ = ϕ ( x ∗ ) \boldsymbol{\phi}_{*} = \boldsymbol{\phi}(\boldsymbol{x}_{*}) ϕ ∗ = ϕ ( x ∗ ) とし、y ∗ = w T ϕ ∗ + ε ∗ y_{*} = \boldsymbol{w}^{\mathsf{T}}\boldsymbol{\phi}_{*} + \varepsilon_{*} y ∗ = w T ϕ ∗ + ε ∗ 、ε ∗ ∼ N ( 0 , σ 2 ) \varepsilon_{*}\sim\mathcal{N}(0,\sigma^{2}) ε ∗ ∼ N ( 0 , σ 2 ) は w \boldsymbol{w} w および既存の雑音と独立とする。このとき
p ( y ∗ ∣ x ∗ , D ) = N ( y ∗ ∣ m N T ϕ ∗ , σ 2 + ϕ ∗ T S N ϕ ∗ ) p(y_{*}\mid\boldsymbol{x}_{*},\mathcal{D}) = \mathcal{N}\bigl(y_{*} \,\big|\, \boldsymbol{m}_N^{\mathsf{T}}\boldsymbol{\phi}_{*},\; \sigma^{2} + \boldsymbol{\phi}_{*}^{\mathsf{T}}S_N\boldsymbol{\phi}_{*}\bigr) p ( y ∗ ∣ x ∗ , D ) = N ( y ∗ m N T ϕ ∗ , σ 2 + ϕ ∗ T S N ϕ ∗ ) である。すなわち予測分散は
σ 2 ⏟ 観測ノイズ(偶然的不確実性) + ϕ ∗ T S N ϕ ∗ ⏟ パラメータの不確実性(認識的不確実性) \underbrace{\sigma^{2}}_{\text{観測ノイズ(偶然的不確実性)}} \;+\; \underbrace{\boldsymbol{\phi}_{*}^{\mathsf{T}}S_N\boldsymbol{\phi}_{*}}_{\text{パラメータの不確実性(認識的不確実性)}} 観測ノイズ(偶然的不確実性) σ 2 + パラメータの不確実性(認識的不確実性) ϕ ∗ T S N ϕ ∗ と分解される。
証明(定理 6.4) D \mathcal{D} D を与えたとき、定理 6.1 より w ∼ N ( m N , S N ) \boldsymbol{w}\sim\mathcal{N}(\boldsymbol{m}_N,S_N) w ∼ N ( m N , S N ) です。多変量ガウス分布の線形像はガウス分布です。実際、定ベクトル a \boldsymbol{a} a に対する特性関数は
E [ e i t a T w ] = exp ( i t a T m N − t 2 2 a T S N a ) \mathbb{E}\bigl[e^{it\boldsymbol{a}^{\mathsf{T}}\boldsymbol{w}}\bigr]
= \exp\!\left(it\,\boldsymbol{a}^{\mathsf{T}}\boldsymbol{m}_N - \frac{t^{2}}{2}\boldsymbol{a}^{\mathsf{T}}S_N\boldsymbol{a}\right) E [ e i t a T w ] = exp ( i t a T m N − 2 t 2 a T S N a ) であり、これは N ( a T m N , a T S N a ) \mathcal{N}(\boldsymbol{a}^{\mathsf{T}}\boldsymbol{m}_N,\, \boldsymbol{a}^{\mathsf{T}}S_N\boldsymbol{a}) N ( a T m N , a T S N a ) の特性関数そのものです。a = ϕ ∗ \boldsymbol{a}=\boldsymbol{\phi}_{*} a = ϕ ∗ と取って
w T ϕ ∗ ∣ D ∼ N ( m N T ϕ ∗ , ϕ ∗ T S N ϕ ∗ ) . \boldsymbol{w}^{\mathsf{T}}\boldsymbol{\phi}_{*} \mid \mathcal{D} \;\sim\; \mathcal{N}\bigl(\boldsymbol{m}_N^{\mathsf{T}}\boldsymbol{\phi}_{*},\; \boldsymbol{\phi}_{*}^{\mathsf{T}}S_N\boldsymbol{\phi}_{*}\bigr). w T ϕ ∗ ∣ D ∼ N ( m N T ϕ ∗ , ϕ ∗ T S N ϕ ∗ ) . 仮定より ε ∗ ∼ N ( 0 , σ 2 ) \varepsilon_{*}\sim\mathcal{N}(0,\sigma^{2}) ε ∗ ∼ N ( 0 , σ 2 ) はこれと独立なので、独立なガウス確率変数の和は平均と分散をそれぞれ足したガウス分布に従います(特性関数が積になることから直ちに従います)。よって
y ∗ = w T ϕ ∗ + ε ∗ ∣ D ∼ N ( m N T ϕ ∗ , σ 2 + ϕ ∗ T S N ϕ ∗ ) y_{*} = \boldsymbol{w}^{\mathsf{T}}\boldsymbol{\phi}_{*} + \varepsilon_{*} \mid \mathcal{D}
\;\sim\; \mathcal{N}\bigl(\boldsymbol{m}_N^{\mathsf{T}}\boldsymbol{\phi}_{*},\; \sigma^{2}+\boldsymbol{\phi}_{*}^{\mathsf{T}}S_N\boldsymbol{\phi}_{*}\bigr) y ∗ = w T ϕ ∗ + ε ∗ ∣ D ∼ N ( m N T ϕ ∗ , σ 2 + ϕ ∗ T S N ϕ ∗ ) です。これは 定義 6.3 の積分を実行したものにほかなりません。
∎
この分解が実務上いちばん重要な結論です。σ 2 \sigma^{2} σ 2 はデータをいくら集めても減りません (コインの表裏は永遠に当てられません)。一方 ϕ ∗ T S N ϕ ∗ \boldsymbol{\phi}_{*}^{\mathsf{T}}S_N\boldsymbol{\phi}_{*} ϕ ∗ T S N ϕ ∗ は、S N − 1 = S 0 − 1 + σ − 2 Φ T Φ S_N^{-1} = S_0^{-1}+\sigma^{-2}\Phi^{\mathsf{T}}\Phi S N − 1 = S 0 − 1 + σ − 2 Φ T Φ がデータとともに増えるので N N N が増えれば減っていきます。前者を偶然的(アレアトリック)不確実性、後者を認識的(エピステミック)不確実性と呼び、後者だけが「もっとデータを集めれば解消する」種類の不確実性です。能動学習でどのデータにラベルを付けるかを選ぶとき、指標にするのは後者です。
例 6.5 (数値例:外挿すると分散が増える )
M = 1 M=1 M = 1 、ϕ ( x ) = x \boldsymbol{\phi}(x)=x ϕ ( x ) = x (原点を通る直線 y = w x y = wx y = w x )とし、σ 2 = 1 \sigma^{2}=1 σ 2 = 1 、事前分布は w ∼ N ( 0 , τ 2 ) w\sim\mathcal{N}(0,\tau^{2}) w ∼ N ( 0 , τ 2 ) 、τ 2 = 1 \tau^{2}=1 τ 2 = 1 とします。データは ( x 1 , y 1 ) = ( 1 , 2 ) (x_1,y_1)=(1,2) ( x 1 , y 1 ) = ( 1 , 2 ) 、( x 2 , y 2 ) = ( 2 , 3 ) (x_2,y_2)=(2,3) ( x 2 , y 2 ) = ( 2 , 3 ) の 2 点です。
計画行列は Φ = ( 1 2 ) \Phi = \begin{pmatrix}1\\2\end{pmatrix} Φ = ( 1 2 ) 、y = ( 2 3 ) \boldsymbol{y}=\begin{pmatrix}2\\3\end{pmatrix} y = ( 2 3 ) なので
Φ T Φ = 1 2 + 2 2 = 5 , Φ T y = 1 ⋅ 2 + 2 ⋅ 3 = 8. \Phi^{\mathsf{T}}\Phi = 1^{2}+2^{2} = 5, \qquad \Phi^{\mathsf{T}}\boldsymbol{y} = 1\cdot 2 + 2\cdot 3 = 8. Φ T Φ = 1 2 + 2 2 = 5 , Φ T y = 1 ⋅ 2 + 2 ⋅ 3 = 8. 最尤推定は 命題 3.3 より w ^ M L = 8 / 5 = 1.6 \hat{w}_{\mathrm{ML}} = 8/5 = 1.6 w ^ ML = 8/5 = 1.6 です。事後分布は 定理 6.1 より
S N − 1 = 1 τ 2 + 1 σ 2 Φ T Φ = 1 + 5 = 6 , S N = 1 6 , S_N^{-1} = \frac{1}{\tau^{2}} + \frac{1}{\sigma^{2}}\Phi^{\mathsf{T}}\Phi = 1 + 5 = 6, \qquad S_N = \frac{1}{6}, S N − 1 = τ 2 1 + σ 2 1 Φ T Φ = 1 + 5 = 6 , S N = 6 1 , m N = S N ( 0 + 1 σ 2 Φ T y ) = 1 6 ⋅ 8 = 4 3 ≈ 1.333. m_N = S_N\left(0 + \frac{1}{\sigma^{2}}\Phi^{\mathsf{T}}\boldsymbol{y}\right) = \frac{1}{6}\cdot 8 = \frac{4}{3} \approx 1.333. m N = S N ( 0 + σ 2 1 Φ T y ) = 6 1 ⋅ 8 = 3 4 ≈ 1.333. 検算として 系 6.2 を使うと、λ = σ 2 / τ 2 = 1 \lambda = \sigma^{2}/\tau^{2}=1 λ = σ 2 / τ 2 = 1 で ( 5 + 1 ) − 1 ⋅ 8 = 4 / 3 (5+1)^{-1}\cdot 8 = 4/3 ( 5 + 1 ) − 1 ⋅ 8 = 4/3 となり一致します。
事後予測分布を 定理 6.4 で計算します(ϕ ∗ = x ∗ \boldsymbol{\phi}_{*}=x_{*} ϕ ∗ = x ∗ なので ϕ ∗ T S N ϕ ∗ = x ∗ 2 / 6 \boldsymbol{\phi}_{*}^{\mathsf{T}}S_N\boldsymbol{\phi}_{*} = x_{*}^{2}/6 ϕ ∗ T S N ϕ ∗ = x ∗ 2 /6 )。
x ∗ x_{*} x ∗ 予測平均 4 3 x ∗ \tfrac{4}{3}x_{*} 3 4 x ∗ 認識的分散 x ∗ 2 / 6 x_{*}^{2}/6 x ∗ 2 /6 全分散 標準偏差 95% 予測区間 1 1 1 1.333 1.333 1.333 0.167 0.167 0.167 1.167 1.167 1.167 1.080 1.080 1.080 [ − 0.78 , 3.45 ] [-0.78,\ 3.45] [ − 0.78 , 3.45 ] 3 3 3 4.000 4.000 4.000 1.500 1.500 1.500 2.500 2.500 2.500 1.581 1.581 1.581 [ 0.90 , 7.10 ] [0.90,\ 7.10] [ 0.90 , 7.10 ] 10 10 10 13.333 13.333 13.333 16.667 16.667 16.667 17.667 17.667 17.667 4.203 4.203 4.203 [ 5.09 , 21.57 ] [5.09,\ 21.57] [ 5.09 , 21.57 ]
(95 % 95\% 95% 区間は平均 ± 1.96 × \pm 1.96 \times ± 1.96 × 標準偏差で計算しました。たとえば x ∗ = 10 x_{*}=10 x ∗ = 10 では 13.333 ± 1.96 × 4.203 = 13.333 ± 8.238 13.333 \pm 1.96\times 4.203 = 13.333 \pm 8.238 13.333 ± 1.96 × 4.203 = 13.333 ± 8.238 です。)
もし点推定 w ^ = 4 / 3 \hat{w}=4/3 w ^ = 4/3 を使ったプラグイン予測をしていたら、どの x ∗ x_{*} x ∗ でも分散は σ 2 = 1 \sigma^{2}=1 σ 2 = 1 、区間幅は ± 1.96 \pm 1.96 ± 1.96 の一定です。x ∗ = 10 x_{*}=10 x ∗ = 10 での真の予測区間は幅 ± 8.24 \pm 8.24 ± 8.24 ですから、プラグイン予測は外挿領域で 4 倍以上も自信過剰 になっています。データが x ∈ { 1 , 2 } x\in\{1,2\} x ∈ { 1 , 2 } 付近にしかないのに x = 10 x=10 x = 10 を当てにいっている、という事実は事後共分散 S N S_N S N の中にしか書かれていません。
0 5 10 x 予測平均 プラグイン予測の 95% 区間 事後予測分布の 95% 区間 データ 上の数値例(データ 2 点)の予測平均 4x/3(実線)と 95% 予測区間(帯)。破線は点推定を使ったプラグイン予測の区間で、幅は一定です。データのない遠方ほど帯が扇状に広がります。
最大化 最頻値 θ について積分 MAP 推定:θ を 1 点に決める(=正則化つき学習)
尤度・事前分布から、最尤推定・MAP 推定・事後予測分布へ至る道筋
演習 7.1 易
例 4.4 の設定(n n n 回中 k k k 回表、事前分布 B e t a ( a , b ) \mathrm{Beta}(a,b) Beta ( a , b ) 、事後分布 B e t a ( k + a , n − k + b ) \mathrm{Beta}(k+a,\,n-k+b) Beta ( k + a , n − k + b ) )で、事後平均が最尤推定量と事前平均の凸結合(重み付き平均)として
E [ θ ∣ D ] = γ ⋅ k n + ( 1 − γ ) ⋅ a a + b , γ = n n + a + b \mathbb{E}[\theta\mid\mathcal{D}] = \gamma\cdot\frac{k}{n} + (1-\gamma)\cdot\frac{a}{a+b}, \qquad \gamma = \frac{n}{n+a+b} E [ θ ∣ D ] = γ ⋅ n k + ( 1 − γ ) ⋅ a + b a , γ = n + a + b n と書けることを示してください。また n → ∞ n\to\infty n → ∞ のときの挙動を述べてください。
解答 B e t a ( α , β ) \mathrm{Beta}(\alpha,\beta) Beta ( α , β ) の平均は α / ( α + β ) \alpha/(\alpha+\beta) α / ( α + β ) なので、事後平均は
E [ θ ∣ D ] = k + a ( k + a ) + ( n − k + b ) = k + a n + a + b \mathbb{E}[\theta\mid\mathcal{D}] = \frac{k+a}{(k+a)+(n-k+b)} = \frac{k+a}{n+a+b} E [ θ ∣ D ] = ( k + a ) + ( n − k + b ) k + a = n + a + b k + a です(分母で k k k が打ち消えます)。一方、右辺を計算すると
γ ⋅ k n + ( 1 − γ ) ⋅ a a + b = n n + a + b ⋅ k n + a + b n + a + b ⋅ a a + b = k n + a + b + a n + a + b = k + a n + a + b \gamma\cdot\frac{k}{n} + (1-\gamma)\cdot\frac{a}{a+b}
= \frac{n}{n+a+b}\cdot\frac{k}{n} + \frac{a+b}{n+a+b}\cdot\frac{a}{a+b}
= \frac{k}{n+a+b} + \frac{a}{n+a+b}
= \frac{k+a}{n+a+b} γ ⋅ n k + ( 1 − γ ) ⋅ a + b a = n + a + b n ⋅ n k + n + a + b a + b ⋅ a + b a = n + a + b k + n + a + b a = n + a + b k + a となり一致します。ここで 1 − γ = 1 − n n + a + b = a + b n + a + b 1-\gamma = 1 - \frac{n}{n+a+b} = \frac{a+b}{n+a+b} 1 − γ = 1 − n + a + b n = n + a + b a + b を使いました。γ ∈ ( 0 , 1 ) \gamma\in(0,1) γ ∈ ( 0 , 1 ) かつ γ + ( 1 − γ ) = 1 \gamma+(1-\gamma)=1 γ + ( 1 − γ ) = 1 なので、これは確かに凸結合です。
n → ∞ n\to\infty n → ∞ のとき γ = n n + a + b → 1 \gamma = \frac{n}{n+a+b}\to 1 γ = n + a + b n → 1 なので、事後平均は最尤推定量 k / n k/n k / n に近づきます。事前分布の影響は O ( 1 / n ) O(1/n) O ( 1/ n ) で消えていきます。逆に n = 0 n=0 n = 0 (データなし)なら γ = 0 \gamma=0 γ = 0 で事後平均は事前平均 a / ( a + b ) a/(a+b) a / ( a + b ) そのものです。
a a a と b b b が「疑似観測回数」と呼ばれる理由もここにあります。事前分布 B e t a ( a , b ) \mathrm{Beta}(a,b) Beta ( a , b ) は、実データに先立って表を a a a 回、裏を b b b 回見たのと同じ効き方をします。a = b = 1 a=b=1 a = b = 1 (一様事前分布)のとき事後平均は ( k + 1 ) / ( n + 2 ) (k+1)/(n+2) ( k + 1 ) / ( n + 2 ) で、これはラプラスの継起の法則、機械学習の言葉では加算スムージングです。
演習 7.2 標準
X 1 , … , X n X_1,\ldots,X_n X 1 , … , X n が独立に平均 λ > 0 \lambda > 0 λ > 0 のポアソン分布 p ( x ∣ λ ) = e − λ λ x / x ! p(x\mid\lambda) = e^{-\lambda}\lambda^{x}/x! p ( x ∣ λ ) = e − λ λ x / x ! (x = 0 , 1 , 2 , … x=0,1,2,\ldots x = 0 , 1 , 2 , … )に従うとします。S = ∑ i = 1 n X i S=\sum_{i=1}^{n} X_i S = ∑ i = 1 n X i とおきます。
λ \lambda λ の最尤推定量を求めてください(S = 0 S=0 S = 0 の場合も論じること)。
事前分布としてガンマ分布 p ( λ ) ∝ λ a − 1 e − b λ p(\lambda)\propto \lambda^{a-1}e^{-b\lambda} p ( λ ) ∝ λ a − 1 e − bλ (a > 0 a>0 a > 0 、b > 0 b>0 b > 0 、λ > 0 \lambda>0 λ > 0 )を取ったときの事後分布を求め、ガンマ分布がポアソン尤度の共役事前分布であることを示してください。
MAP 推定量を求め、n → ∞ n\to\infty n → ∞ での挙動を述べてください。
解答 1. 尤度と対数尤度は
L ( λ ) = ∏ i = 1 n e − λ λ x i x i ! = e − n λ λ S ∏ i 1 x i ! , ℓ ( λ ) = − n λ + S log λ − ∑ i log ( x i ! ) . L(\lambda) = \prod_{i=1}^{n} \frac{e^{-\lambda}\lambda^{x_i}}{x_i!}
= e^{-n\lambda}\lambda^{S}\prod_{i}\frac{1}{x_i!},
\qquad
\ell(\lambda) = -n\lambda + S\log\lambda - \sum_i \log(x_i!). L ( λ ) = i = 1 ∏ n x i ! e − λ λ x i = e − nλ λ S i ∏ x i ! 1 , ℓ ( λ ) = − nλ + S log λ − i ∑ log ( x i !) . S ≥ 1 S \ge 1 S ≥ 1 のとき ℓ ′ ( λ ) = − n + S / λ = 0 \ell'(\lambda) = -n + S/\lambda = 0 ℓ ′ ( λ ) = − n + S / λ = 0 より λ = S / n \lambda = S/n λ = S / n 、また ℓ ′ ′ ( λ ) = − S / λ 2 < 0 \ell''(\lambda) = -S/\lambda^{2} < 0 ℓ ′′ ( λ ) = − S / λ 2 < 0 なので ℓ \ell ℓ は狭義凹で、λ ^ M L = S / n = x ˉ \hat{\lambda}_{\mathrm{ML}} = S/n = \bar{x} λ ^ ML = S / n = x ˉ が唯一の最大点です。
S = 0 S=0 S = 0 のときは ℓ ( λ ) = − n λ \ell(\lambda) = -n\lambda ℓ ( λ ) = − nλ が λ > 0 \lambda>0 λ > 0 で狭義単調減少なので、( 0 , ∞ ) (0,\infty) ( 0 , ∞ ) 上に最大点は存在せず、λ ↓ 0 \lambda\downarrow 0 λ ↓ 0 が上限を与えます。パラメータ空間を [ 0 , ∞ ) [0,\infty) [ 0 , ∞ ) に広げれば λ ^ M L = 0 \hat{\lambda}_{\mathrm{ML}}=0 λ ^ ML = 0 で、これも x ˉ = 0 \bar{x}=0 x ˉ = 0 と書けます。例 3.2 の端の場合と同じ病理です。
2. 定理 2.2 の比例形より
p ( λ ∣ D ) ∝ ( e − n λ λ S ) ⋅ ( λ a − 1 e − b λ ) = λ ( S + a ) − 1 e − ( n + b ) λ . p(\lambda\mid\mathcal{D}) \propto \bigl(e^{-n\lambda}\lambda^{S}\bigr)\cdot\bigl(\lambda^{a-1}e^{-b\lambda}\bigr)
= \lambda^{(S+a)-1}\,e^{-(n+b)\lambda}. p ( λ ∣ D ) ∝ ( e − nλ λ S ) ⋅ ( λ a − 1 e − bλ ) = λ ( S + a ) − 1 e − ( n + b ) λ . 右辺は形状 S + a S+a S + a 、率 n + b n+b n + b のガンマ分布 G a m m a ( S + a , n + b ) \mathrm{Gamma}(S+a,\, n+b) Gamma ( S + a , n + b ) の密度の λ \lambda λ 依存部分です。S + a > 0 S+a>0 S + a > 0 、n + b > 0 n+b>0 n + b > 0 なのでこれは正しく正規化可能な密度であり、事後分布は G a m m a ( S + a , n + b ) \mathrm{Gamma}(S+a,\,n+b) Gamma ( S + a , n + b ) です。事後分布が再びガンマ分布族に属するので、定義 4.3 の意味で共役です。更新則は「形状に観測値の総和を足し、率に標本数を足す」となります。
3. 対数事後分布は ( S + a − 1 ) log λ − ( n + b ) λ (S+a-1)\log\lambda - (n+b)\lambda ( S + a − 1 ) log λ − ( n + b ) λ (定数を除く)で、微分すると S + a − 1 λ − ( n + b ) \frac{S+a-1}{\lambda} - (n+b) λ S + a − 1 − ( n + b ) です。S + a > 1 S+a>1 S + a > 1 のとき、これが 0 0 0 になる λ \lambda λ は
λ ^ M A P = S + a − 1 n + b \hat{\lambda}_{\mathrm{MAP}} = \frac{S+a-1}{n+b} λ ^ MAP = n + b S + a − 1 であり、二階微分 − ( S + a − 1 ) / λ 2 < 0 -(S+a-1)/\lambda^{2}<0 − ( S + a − 1 ) / λ 2 < 0 より唯一の最大点です。S + a ≤ 1 S+a\le 1 S + a ≤ 1 のときは対数事後分布が単調減少なので最頻値は λ = 0 \lambda=0 λ = 0 です。
n → ∞ n\to\infty n → ∞ の挙動は、x ˉ = S / n \bar{x}=S/n x ˉ = S / n を使って
λ ^ M A P = n x ˉ + a − 1 n + b = x ˉ ⋅ n n + b + a − 1 n + b \hat{\lambda}_{\mathrm{MAP}} = \frac{n\bar{x}+a-1}{n+b} = \bar{x}\cdot\frac{n}{n+b} + \frac{a-1}{n+b} λ ^ MAP = n + b n x ˉ + a − 1 = x ˉ ⋅ n + b n + n + b a − 1 と書き直せば見えます。x ˉ \bar{x} x ˉ が大数の法則で真の λ 0 \lambda_0 λ 0 に収束するなら、第 1 項は λ 0 \lambda_0 λ 0 に、第 2 項は 0 0 0 に収束するので λ ^ M A P → λ 0 \hat{\lambda}_{\mathrm{MAP}}\to\lambda_0 λ ^ MAP → λ 0 です。事前分布の影響は O ( 1 / n ) O(1/n) O ( 1/ n ) で消え、演習 7.1 と同じ結論になります。
演習 7.3 標準
λ > 0 \lambda>0 λ > 0 、z ∈ R z\in\mathbb{R} z ∈ R を定数とし、1 変数関数
J ( w ) = 1 2 ( w − z ) 2 + λ ∣ w ∣ ( w ∈ R ) J(w) = \frac{1}{2}(w-z)^{2} + \lambda|w| \qquad (w\in\mathbb{R}) J ( w ) = 2 1 ( w − z ) 2 + λ ∣ w ∣ ( w ∈ R ) を考えます。
J J J の最小点が唯一存在し、それがソフトしきい値関数 S λ ( z ) = sign ( z ) max ( ∣ z ∣ − λ , 0 ) S_\lambda(z) = \operatorname{sign}(z)\max(|z|-\lambda,\,0) S λ ( z ) = sign ( z ) max ( ∣ z ∣ − λ , 0 ) で与えられることを示してください。
L1 罰則を L2 罰則に替えた J 2 ( w ) = 1 2 ( w − z ) 2 + λ 2 w 2 J_2(w) = \frac{1}{2}(w-z)^{2} + \frac{\lambda}{2}w^{2} J 2 ( w ) = 2 1 ( w − z ) 2 + 2 λ w 2 の最小点を求め、1 と比較して「なぜ L1 だけが厳密な 0 0 0 を生むのか」を説明してください。
解答 1. 1 2 ( w − z ) 2 \frac12(w-z)^2 2 1 ( w − z ) 2 は狭義凸、λ ∣ w ∣ \lambda|w| λ ∣ w ∣ は凸なので J J J は狭義凸です。したがって最小点は高々 1 つです。また λ ∣ w ∣ ≥ 0 \lambda|w|\ge 0 λ ∣ w ∣ ≥ 0 より J ( w ) ≥ 1 2 ( w − z ) 2 J(w)\ge\frac{1}{2}(w-z)^{2} J ( w ) ≥ 2 1 ( w − z ) 2 であり、∣ w ∣ → ∞ |w|\to\infty ∣ w ∣ → ∞ のとき右辺が ∞ \infty ∞ に発散するので J J J は強制的です。連続な強制的関数は最小値を取るので、最小点はちょうど 1 つ存在します。凸関数では w ⋆ w^{\star} w ⋆ が最小点であることと 0 ∈ ∂ J ( w ⋆ ) 0\in\partial J(w^{\star}) 0 ∈ ∂ J ( w ⋆ ) が同値です。ここで劣微分は
∂ J ( w ) = { w − z + λ s : s ∈ ∂ ∣ w ∣ } , ∂ ∣ w ∣ = { { sign ( w ) } ( w ≠ 0 ) [ − 1 , 1 ] ( w = 0 ) \partial J(w) = \{\,w - z + \lambda s \;:\; s\in\partial|w|\,\},
\qquad
\partial|w| = \begin{cases} \{\operatorname{sign}(w)\} & (w\ne 0)\\ [-1,1] & (w=0)\end{cases} ∂ J ( w ) = { w − z + λ s : s ∈ ∂ ∣ w ∣ } , ∂ ∣ w ∣ = { { sign ( w )} [ − 1 , 1 ] ( w = 0 ) ( w = 0 ) です。場合分けします。
w > 0 w>0 w > 0 :条件は w − z + λ = 0 w-z+\lambda=0 w − z + λ = 0 、すなわち w = z − λ w = z-\lambda w = z − λ 。これが w > 0 w>0 w > 0 を満たすのは z > λ z>\lambda z > λ のときだけです。
w < 0 w<0 w < 0 :条件は w − z − λ = 0 w-z-\lambda=0 w − z − λ = 0 、すなわち w = z + λ w = z+\lambda w = z + λ 。これが w < 0 w<0 w < 0 を満たすのは z < − λ z<-\lambda z < − λ のときだけです。
w = 0 w=0 w = 0 :条件は 0 ∈ { − z + λ s : s ∈ [ − 1 , 1 ] } = [ − z − λ , − z + λ ] 0\in\{-z+\lambda s: s\in[-1,1]\} = [-z-\lambda,\, -z+\lambda] 0 ∈ { − z + λ s : s ∈ [ − 1 , 1 ]} = [ − z − λ , − z + λ ] 、すなわち − z − λ ≤ 0 ≤ − z + λ -z-\lambda\le 0\le -z+\lambda − z − λ ≤ 0 ≤ − z + λ 、整理して ∣ z ∣ ≤ λ |z|\le\lambda ∣ z ∣ ≤ λ 。
3 つの場合は z z z について排他的かつ網羅的で、まとめると
w ⋆ = { z − λ ( z > λ ) 0 ( ∣ z ∣ ≤ λ ) z + λ ( z < − λ ) w^{\star} = \begin{cases} z-\lambda & (z>\lambda)\\ 0 & (|z|\le\lambda)\\ z+\lambda & (z<-\lambda)\end{cases} w ⋆ = ⎩ ⎨ ⎧ z − λ 0 z + λ ( z > λ ) ( ∣ z ∣ ≤ λ ) ( z < − λ ) です。z > λ z>\lambda z > λ のとき sign ( z ) max ( ∣ z ∣ − λ , 0 ) = z − λ \operatorname{sign}(z)\max(|z|-\lambda,0) = z-\lambda sign ( z ) max ( ∣ z ∣ − λ , 0 ) = z − λ 、z < − λ z<-\lambda z < − λ のとき − ( ∣ z ∣ − λ ) = − ( − z − λ ) = z + λ -(|z|-\lambda) = -(-z-\lambda) = z+\lambda − ( ∣ z ∣ − λ ) = − ( − z − λ ) = z + λ 、∣ z ∣ ≤ λ |z|\le\lambda ∣ z ∣ ≤ λ のとき max ( ∣ z ∣ − λ , 0 ) = 0 \max(|z|-\lambda,0)=0 max ( ∣ z ∣ − λ , 0 ) = 0 なので、これは S λ ( z ) S_\lambda(z) S λ ( z ) に一致します。
2. J 2 J_2 J 2 は滑らかで狭義凸なので J 2 ′ ( w ) = ( w − z ) + λ w = 0 J_2'(w) = (w-z)+\lambda w = 0 J 2 ′ ( w ) = ( w − z ) + λ w = 0 、すなわち w ⋆ = z / ( 1 + λ ) w^{\star} = z/(1+\lambda) w ⋆ = z / ( 1 + λ ) が唯一の最小点です。これが 0 0 0 になるのは z = 0 z=0 z = 0 のときだけです。
違いは原点での「傾きの跳び」にあります。λ ∣ w ∣ \lambda|w| λ ∣ w ∣ の劣微分は w = 0 w=0 w = 0 で幅 2 λ 2\lambda 2 λ の区間 [ − λ , λ ] [-\lambda,\lambda] [ − λ , λ ] を持つので、データ由来の勾配 − z -z − z がこの区間に収まっている限り w = 0 w=0 w = 0 が最適であり続けます。∣ z ∣ ≤ λ |z|\le\lambda ∣ z ∣ ≤ λ という「不感帯」が生じるのはこのためです。一方 λ 2 w 2 \frac{\lambda}{2}w^{2} 2 λ w 2 の微分 λ w \lambda w λ w は w = 0 w=0 w = 0 で 0 0 0 になり、罰則が 0 0 0 の近傍で w w w を押し戻す力を失うため、どんなに小さくても z ≠ 0 z\ne0 z = 0 なら w ⋆ ≠ 0 w^\star\ne0 w ⋆ = 0 になります。
事前分布の言葉では、命題 5.3 のラプラス密度が原点に尖った角を持つのに対し、定理 5.1 のガウス密度は原点で滑らかである、という違いに対応します。ただし 注意 5.4 で述べたとおり、これはあくまで最頻値の性質です。
演習 7.4 難
ロジスティック回帰 を扱います。σ ( t ) = 1 / ( 1 + e − t ) \sigma(t)=1/(1+e^{-t}) σ ( t ) = 1/ ( 1 + e − t ) 、データ ( x i , y i ) i = 1 n (\boldsymbol{x}_i,y_i)_{i=1}^{n} ( x i , y i ) i = 1 n (x i ∈ R M \boldsymbol{x}_i\in\mathbb{R}^{M} x i ∈ R M 、y i ∈ { 0 , 1 } y_i\in\{0,1\} y i ∈ { 0 , 1 } )に対しモデルを p ( y = 1 ∣ x , w ) = σ ( w T x ) p(y=1\mid\boldsymbol{x},\boldsymbol{w}) = \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}) p ( y = 1 ∣ x , w ) = σ ( w T x ) とします。
負の対数尤度が交差エントロピー損失 E ( w ) = − ∑ i [ y i log p ^ i + ( 1 − y i ) log ( 1 − p ^ i ) ] E(\boldsymbol{w}) = -\sum_i\bigl[y_i\log\hat{p}_i + (1-y_i)\log(1-\hat{p}_i)\bigr] E ( w ) = − ∑ i [ y i log p ^ i + ( 1 − y i ) log ( 1 − p ^ i ) ] (p ^ i = σ ( w T x i ) \hat p_i = \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i) p ^ i = σ ( w T x i ) )に一致することを示し、さらに t i = ( 2 y i − 1 ) w T x i t_i = (2y_i-1)\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i t i = ( 2 y i − 1 ) w T x i とおくと E ( w ) = ∑ i log ( 1 + e − t i ) E(\boldsymbol{w}) = \sum_i \log(1+e^{-t_i}) E ( w ) = ∑ i log ( 1 + e − t i ) と書けることを示してください。
事前分布 w ∼ N ( 0 , τ 2 I ) \boldsymbol{w}\sim\mathcal{N}(\boldsymbol{0},\tau^{2}I) w ∼ N ( 0 , τ 2 I ) の下で MAP 推定が E ( w ) + 1 2 τ 2 ∥ w ∥ 2 E(\boldsymbol{w}) + \frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2} E ( w ) + 2 τ 2 1 ∥ w ∥ 2 の最小化になることを示してください。
データが線形分離可能、すなわちある w 0 \boldsymbol{w}_0 w 0 (∥ w 0 ∥ = 1 \lVert\boldsymbol{w}_0\rVert=1 ∥ w 0 ∥ = 1 )が存在してすべての i i i で ( 2 y i − 1 ) w 0 T x i > 0 (2y_i-1)\boldsymbol{w}_0^{\mathsf{T}}\boldsymbol{x}_i>0 ( 2 y i − 1 ) w 0 T x i > 0 が成り立つとします。このとき最尤推定量は存在しないが、2 の MAP 推定量はただ 1 つ存在することを示してください。
解答 1. y i ∈ { 0 , 1 } y_i\in\{0,1\} y i ∈ { 0 , 1 } なので p ( y i ∣ x i , w ) = p ^ i y i ( 1 − p ^ i ) 1 − y i p(y_i\mid\boldsymbol{x}_i,\boldsymbol{w}) = \hat p_i^{\,y_i}(1-\hat p_i)^{1-y_i} p ( y i ∣ x i , w ) = p ^ i y i ( 1 − p ^ i ) 1 − y i です(y i = 1 y_i=1 y i = 1 なら p ^ i \hat p_i p ^ i 、y i = 0 y_i=0 y i = 0 なら 1 − p ^ i 1-\hat p_i 1 − p ^ i になります)。各データが独立なので対数尤度は和になり、
log p ( D ∣ w ) = ∑ i [ y i log p ^ i + ( 1 − y i ) log ( 1 − p ^ i ) ] = − E ( w ) . \log p(\mathcal{D}\mid\boldsymbol{w}) = \sum_i \bigl[y_i\log\hat p_i + (1-y_i)\log(1-\hat p_i)\bigr] = -E(\boldsymbol{w}). log p ( D ∣ w ) = i ∑ [ y i log p ^ i + ( 1 − y i ) log ( 1 − p ^ i ) ] = − E ( w ) . 次に 1 − σ ( t ) = 1 − 1 1 + e − t = e − t 1 + e − t = 1 1 + e t = σ ( − t ) 1-\sigma(t) = 1 - \frac{1}{1+e^{-t}} = \frac{e^{-t}}{1+e^{-t}} = \frac{1}{1+e^{t}} = \sigma(-t) 1 − σ ( t ) = 1 − 1 + e − t 1 = 1 + e − t e − t = 1 + e t 1 = σ ( − t ) に注意します。y i = 1 y_i=1 y i = 1 なら p ^ i = σ ( w T x i ) = σ ( t i ) \hat p_i = \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i) = \sigma(t_i) p ^ i = σ ( w T x i ) = σ ( t i ) 、y i = 0 y_i=0 y i = 0 なら 1 − p ^ i = σ ( − w T x i ) = σ ( t i ) 1-\hat p_i = \sigma(-\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i) = \sigma(t_i) 1 − p ^ i = σ ( − w T x i ) = σ ( t i ) (2 y i − 1 = − 1 2y_i-1=-1 2 y i − 1 = − 1 だから)なので、どちらの場合も i i i 番目の項は − log σ ( t i ) -\log\sigma(t_i) − log σ ( t i ) です。− log σ ( t ) = log ( 1 + e − t ) -\log\sigma(t) = \log(1+e^{-t}) − log σ ( t ) = log ( 1 + e − t ) なので E ( w ) = ∑ i log ( 1 + e − t i ) E(\boldsymbol{w}) = \sum_i\log(1+e^{-t_i}) E ( w ) = ∑ i log ( 1 + e − t i ) を得ます。
2. 命題 4.2 より MAP は log p ( D ∣ w ) + log p ( w ) \log p(\mathcal{D}\mid\boldsymbol{w}) + \log p(\boldsymbol{w}) log p ( D ∣ w ) + log p ( w ) の最大点です。1 より第 1 項は − E ( w ) -E(\boldsymbol{w}) − E ( w ) 、第 2 項は 定理 5.1 の証明第 1 段と同じ計算で c o n s t − 1 2 τ 2 ∥ w ∥ 2 \mathrm{const} - \frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2} const − 2 τ 2 1 ∥ w ∥ 2 です。符号を反転して定数を落とせば、F ( w ) = E ( w ) + 1 2 τ 2 ∥ w ∥ 2 F(\boldsymbol{w}) = E(\boldsymbol{w}) + \frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2} F ( w ) = E ( w ) + 2 τ 2 1 ∥ w ∥ 2 の最小化になります。これは本文と同じ E ( w ) + λ ∥ w ∥ 2 E(\boldsymbol{w}) + \lambda\lVert\boldsymbol{w}\rVert^{2} E ( w ) + λ ∥ w ∥ 2 という書き方をすれば、正則化係数 λ = 1 / ( 2 τ 2 ) \lambda = 1/(2\tau^{2}) λ = 1/ ( 2 τ 2 ) の L2 正則化つきロジスティック回帰にほかなりません。
3. 最尤推定量の非存在。 γ = min i ( 2 y i − 1 ) w 0 T x i \gamma = \min_i (2y_i-1)\boldsymbol{w}_0^{\mathsf{T}}\boldsymbol{x}_i γ = min i ( 2 y i − 1 ) w 0 T x i とおくと、有限個の正数の最小値なので γ > 0 \gamma>0 γ > 0 です。w = c w 0 \boldsymbol{w}=c\boldsymbol{w}_0 w = c w 0 (c > 0 c>0 c > 0 )と取ると t i = c ( 2 y i − 1 ) w 0 T x i ≥ c γ t_i = c(2y_i-1)\boldsymbol{w}_0^{\mathsf{T}}\boldsymbol{x}_i \ge c\gamma t i = c ( 2 y i − 1 ) w 0 T x i ≥ c γ なので、log ( 1 + e − t ) \log(1+e^{-t}) log ( 1 + e − t ) が t t t について減少することから
E ( c w 0 ) = ∑ i = 1 n log ( 1 + e − t i ) ≤ n log ( 1 + e − c γ ) → c → ∞ 0. E(c\boldsymbol{w}_0) = \sum_{i=1}^n \log(1+e^{-t_i}) \le n\log\bigl(1+e^{-c\gamma}\bigr) \xrightarrow[c\to\infty]{} 0. E ( c w 0 ) = i = 1 ∑ n log ( 1 + e − t i ) ≤ n log ( 1 + e − c γ ) c → ∞ 0. 一方、任意の有限な w \boldsymbol{w} w に対し e − t i > 0 e^{-t_i}>0 e − t i > 0 なので log ( 1 + e − t i ) > 0 \log(1+e^{-t_i})>0 log ( 1 + e − t i ) > 0 、したがって E ( w ) > 0 E(\boldsymbol{w})>0 E ( w ) > 0 です。よって inf w E = 0 \inf_{\boldsymbol{w}} E = 0 inf w E = 0 ですが、この下限を達成する w \boldsymbol{w} w は存在しません。E E E の最小化は尤度の最大化と同値なので、最尤推定量は存在しません(数値的には重みが発散します)。これは 定理 6.3[ロジスティック回帰] を、いま用意した事前分布の言葉で言い直したものです。
MAP 推定量の存在。 E ( w ) ≥ 0 E(\boldsymbol{w})\ge 0 E ( w ) ≥ 0 より
F ( w ) = E ( w ) + 1 2 τ 2 ∥ w ∥ 2 ≥ 1 2 τ 2 ∥ w ∥ 2 F(\boldsymbol{w}) = E(\boldsymbol{w}) + \frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2} \ge \frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2} F ( w ) = E ( w ) + 2 τ 2 1 ∥ w ∥ 2 ≥ 2 τ 2 1 ∥ w ∥ 2 なので、∥ w ∥ → ∞ \lVert\boldsymbol{w}\rVert\to\infty ∥ w ∥ → ∞ のとき F ( w ) → ∞ F(\boldsymbol{w})\to\infty F ( w ) → ∞ (強制的)です。いま F ( 0 ) = E ( 0 ) = n log 2 F(\boldsymbol{0}) = E(\boldsymbol{0}) = n\log 2 F ( 0 ) = E ( 0 ) = n log 2 なので、劣位集合 C = { w : F ( w ) ≤ n log 2 } C = \{\boldsymbol{w} : F(\boldsymbol{w})\le n\log 2\} C = { w : F ( w ) ≤ n log 2 } は空でなく、その上では 1 2 τ 2 ∥ w ∥ 2 ≤ n log 2 \frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2}\le n\log 2 2 τ 2 1 ∥ w ∥ 2 ≤ n log 2 すなわち ∥ w ∥ ≤ τ 2 n log 2 \lVert\boldsymbol{w}\rVert\le\tau\sqrt{2n\log 2} ∥ w ∥ ≤ τ 2 n log 2 となるので有界です。F F F は連続なので C C C は閉、よって C C C はコンパクトです。ワイエルシュトラスの定理より F F F は C C C 上で最小値を取り、C C C の外では F > n log 2 ≥ min C F F>n\log 2\ge \min_C F F > n log 2 ≥ min C F なので、それは R M \mathbb{R}^{M} R M 全体での最小値です。
一意性。 g ( t ) = log ( 1 + e − t ) g(t)=\log(1+e^{-t}) g ( t ) = log ( 1 + e − t ) は g ′ ′ ( t ) = e − t ( 1 + e − t ) 2 > 0 g''(t) = \frac{e^{-t}}{(1+e^{-t})^{2}}>0 g ′′ ( t ) = ( 1 + e − t ) 2 e − t > 0 より凸で、t i t_i t i は w \boldsymbol{w} w の線形関数なので、凸関数と線形写像の合成である g ( t i ( w ) ) g(t_i(\boldsymbol{w})) g ( t i ( w )) は凸、その和 E E E も凸です。1 2 τ 2 ∥ w ∥ 2 \frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2} 2 τ 2 1 ∥ w ∥ 2 は τ 2 > 0 \tau^2 > 0 τ 2 > 0 より狭義凸なので、F F F は狭義凸です。狭義凸関数の最小点は高々 1 つなので、MAP 推定量はただ 1 つに定まります。
解釈。 分離可能なデータでは「境界からできるだけ遠ざけたい」という力が無限に働き、最尤推定は破綻します。ガウス事前分布はこの力に 1 2 τ 2 ∥ w ∥ 2 \frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2} 2 τ 2 1 ∥ w ∥ 2 という有限の対価を課し、釣り合いの位置で解を止めます。正則化が「数値的な安定化テクニック」ではなく「モデルの一部」であることが、この例にはっきり現れています。
C. M. Bishop, Pattern Recognition and Machine Learning , Springer, 2006 — 第 1 章(確率とベイズの立場)、第 3 章(ベイズ線形回帰と事後予測分布)。
K. P. Murphy, Probabilistic Machine Learning: An Introduction , MIT Press, 2022 — 第 4 章(最尤推定・MAP 推定・共役事前分布)。著者サイトで公開されています: probml.github.io/pml-book
A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, D. B. Rubin, Bayesian Data Analysis , 3rd ed., CRC Press, 2013 — 第 1〜3 章(ベイズ推論の枠組みと共役事前分布)。
T. Hastie, R. Tibshirani, J. Friedman, The Elements of Statistical Learning , 2nd ed., Springer, 2009 — 第 3 章(リッジ回帰の縮小効果とラッソ)。著者サイトで公開されています: hastie.su.domains/ElemStatLearn
R. Tibshirani, “Regression Shrinkage and Selection via the Lasso”, Journal of the Royal Statistical Society, Series B 58 (1996), 267–288. — ラッソの原論文。ラプラス事前分布による解釈も述べられています。
Y. Gal and Z. Ghahramani, “Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning”, ICML 2016 . arXiv:1506.02142
この補題は何のためにあるか。 定理 6.1 の証明では、事後密度の対数が「w \boldsymbol{w} w の 2 次形式 + + + 1 次形式 + + + 定数」の形になることまでを示しました。そこから「だからガウス分布である」と言うには、その形の関数が実際にどのガウス分布の密度なのかを特定する必要があります。行列版の平方完成がその作業です。
補題 7.5 (2 次形式で書かれた密度の同定 ) ✓
A A A を M M M 次の対称正定値行列、b ∈ R M \boldsymbol{b}\in\mathbb{R}^{M} b ∈ R M 、C > 0 C>0 C > 0 を定数とする。R M \mathbb{R}^{M} R M 上の確率密度 q q q が
q ( w ) = C exp ( − 1 2 w T A w + b T w ) ( w ∈ R M ) q(\boldsymbol{w}) = C \exp\!\left( -\frac{1}{2}\boldsymbol{w}^{\mathsf{T}}A\boldsymbol{w} + \boldsymbol{b}^{\mathsf{T}}\boldsymbol{w} \right)
\qquad (\boldsymbol{w}\in\mathbb{R}^{M}) q ( w ) = C exp ( − 2 1 w T A w + b T w ) ( w ∈ R M ) を満たすならば、q q q は多変量正規分布 N ( A − 1 b , A − 1 ) \mathcal{N}(A^{-1}\boldsymbol{b},\, A^{-1}) N ( A − 1 b , A − 1 ) の密度である。
証明(補題 7.5) A A A は正定値なので可逆で、A − 1 A^{-1} A − 1 も対称正定値です(A = A T A = A^{\mathsf{T}} A = A T より ( A − 1 ) T = ( A T ) − 1 = A − 1 (A^{-1})^{\mathsf{T}} = (A^{\mathsf{T}})^{-1} = A^{-1} ( A − 1 ) T = ( A T ) − 1 = A − 1 、また A A A の固有値がすべて正なら A − 1 A^{-1} A − 1 の固有値もすべて正)。μ = A − 1 b \boldsymbol{\mu} = A^{-1}\boldsymbol{b} μ = A − 1 b とおきます。
まず恒等式を確かめます。展開すると
− 1 2 ( w − μ ) T A ( w − μ ) = − 1 2 ( w T A w − w T A μ − μ T A w + μ T A μ ) -\frac{1}{2}(\boldsymbol{w}-\boldsymbol{\mu})^{\mathsf{T}}A(\boldsymbol{w}-\boldsymbol{\mu})
= -\frac{1}{2}\Bigl( \boldsymbol{w}^{\mathsf{T}}A\boldsymbol{w} - \boldsymbol{w}^{\mathsf{T}}A\boldsymbol{\mu} - \boldsymbol{\mu}^{\mathsf{T}}A\boldsymbol{w} + \boldsymbol{\mu}^{\mathsf{T}}A\boldsymbol{\mu} \Bigr) − 2 1 ( w − μ ) T A ( w − μ ) = − 2 1 ( w T A w − w T A μ − μ T A w + μ T A μ ) ですが、A μ = A A − 1 b = b A\boldsymbol{\mu} = AA^{-1}\boldsymbol{b} = \boldsymbol{b} A μ = A A − 1 b = b なので w T A μ = w T b = b T w \boldsymbol{w}^{\mathsf{T}}A\boldsymbol{\mu} = \boldsymbol{w}^{\mathsf{T}}\boldsymbol{b} = \boldsymbol{b}^{\mathsf{T}}\boldsymbol{w} w T A μ = w T b = b T w (スカラーなので転置しても同じ)、同様に μ T A w = ( A μ ) T w = b T w \boldsymbol{\mu}^{\mathsf{T}}A\boldsymbol{w} = (A\boldsymbol{\mu})^{\mathsf{T}}\boldsymbol{w} = \boldsymbol{b}^{\mathsf{T}}\boldsymbol{w} μ T A w = ( A μ ) T w = b T w (A A A の対称性を使いました)、そして μ T A μ = b T A − 1 b \boldsymbol{\mu}^{\mathsf{T}}A\boldsymbol{\mu} = \boldsymbol{b}^{\mathsf{T}}A^{-1}\boldsymbol{b} μ T A μ = b T A − 1 b です。したがって
− 1 2 ( w − μ ) T A ( w − μ ) = − 1 2 w T A w + b T w − 1 2 b T A − 1 b . -\frac{1}{2}(\boldsymbol{w}-\boldsymbol{\mu})^{\mathsf{T}}A(\boldsymbol{w}-\boldsymbol{\mu})
= -\frac{1}{2}\boldsymbol{w}^{\mathsf{T}}A\boldsymbol{w} + \boldsymbol{b}^{\mathsf{T}}\boldsymbol{w} - \frac{1}{2}\boldsymbol{b}^{\mathsf{T}}A^{-1}\boldsymbol{b}. − 2 1 ( w − μ ) T A ( w − μ ) = − 2 1 w T A w + b T w − 2 1 b T A − 1 b . これを q q q の式に代入すると、K = C exp ( 1 2 b T A − 1 b ) > 0 K = C\exp\bigl(\tfrac12\boldsymbol{b}^{\mathsf{T}}A^{-1}\boldsymbol{b}\bigr) > 0 K = C exp ( 2 1 b T A − 1 b ) > 0 とおいて
q ( w ) = K exp ( − 1 2 ( w − μ ) T A ( w − μ ) ) . q(\boldsymbol{w}) = K \exp\!\left( -\frac{1}{2}(\boldsymbol{w}-\boldsymbol{\mu})^{\mathsf{T}}A(\boldsymbol{w}-\boldsymbol{\mu}) \right). q ( w ) = K exp ( − 2 1 ( w − μ ) T A ( w − μ ) ) . 一方、Σ = A − 1 \Sigma = A^{-1} Σ = A − 1 とした多変量正規分布 N ( μ , Σ ) \mathcal{N}(\boldsymbol{\mu},\Sigma) N ( μ , Σ ) の密度は
g ( w ) = 1 ( 2 π ) M / 2 ( det Σ ) 1 / 2 exp ( − 1 2 ( w − μ ) T Σ − 1 ( w − μ ) ) g(\boldsymbol{w}) = \frac{1}{(2\pi)^{M/2}(\det\Sigma)^{1/2}} \exp\!\left( -\frac{1}{2}(\boldsymbol{w}-\boldsymbol{\mu})^{\mathsf{T}}\Sigma^{-1}(\boldsymbol{w}-\boldsymbol{\mu}) \right) g ( w ) = ( 2 π ) M /2 ( det Σ ) 1/2 1 exp ( − 2 1 ( w − μ ) T Σ − 1 ( w − μ ) ) であり、Σ − 1 = A \Sigma^{-1}=A Σ − 1 = A なので指数部は q q q のそれと完全に一致します。つまり q = ( K / K ′ ) g q = (K/K')\,g q = ( K / K ′ ) g (K ′ K' K ′ は g g g の正規化定数)という比例関係が全点で成り立ちます。
最後に比例定数が 1 1 1 であることを言います。q q q も g g g も確率密度なので ∫ q = ∫ g = 1 \int q = \int g = 1 ∫ q = ∫ g = 1 です。q = ( K / K ′ ) g q = (K/K')g q = ( K / K ′ ) g の両辺を積分すると 1 = ( K / K ′ ) ⋅ 1 1 = (K/K')\cdot 1 1 = ( K / K ′ ) ⋅ 1 、すなわち K = K ′ K = K' K = K ′ を得ます。よって q = g q = g q = g 、つまり q q q は N ( μ , A − 1 ) = N ( A − 1 b , A − 1 ) \mathcal{N}(\boldsymbol{\mu},A^{-1}) = \mathcal{N}(A^{-1}\boldsymbol{b},\,A^{-1}) N ( μ , A − 1 ) = N ( A − 1 b , A − 1 ) の密度です。
∎
1 次元で確かめておきます。 M = 1 M=1 M = 1 、A = a > 0 A=a>0 A = a > 0 、b ∈ R b\in\mathbb{R} b ∈ R なら q ( w ) = C exp ( − 1 2 a w 2 + b w ) q(w) = C\exp(-\tfrac12 a w^{2}+bw) q ( w ) = C exp ( − 2 1 a w 2 + b w ) で、補題は「これは平均 b / a b/a b / a 、分散 1 / a 1/a 1/ a の正規分布」と言っています。実際 − 1 2 a w 2 + b w = − a 2 ( w − b a ) 2 + b 2 2 a -\tfrac12 a w^2 + bw = -\tfrac{a}{2}\bigl(w - \tfrac{b}{a}\bigr)^{2} + \tfrac{b^{2}}{2a} − 2 1 a w 2 + b w = − 2 a ( w − a b ) 2 + 2 a b 2 なので、確かに平均 b / a b/a b / a 、分散 1 / a 1/a 1/ a です。a a a は分散の逆数、すなわち精度です。定理 6.1 の S N − 1 = S 0 − 1 + σ − 2 Φ T Φ S_N^{-1} = S_0^{-1}+\sigma^{-2}\Phi^{\mathsf{T}}\Phi S N − 1 = S 0 − 1 + σ − 2 Φ T Φ が「精度は足し算になる」と読めるのは、このためです。