Skip to content

確率論とベイズ統計:最尤推定・MAP 推定と正則化の正体

Prerequisite:主成分分析:分散最大化はなぜ固有値問題になるのか

Raw

This content is not available in your language yet.

  • **最尤推定(MLE)**は「観測データが最も出やすくなるパラメータ」を選ぶ方法です。データが少ないと極端な答え(確率 11 や無限大の重み)を平然と返します。
  • ベイズの定理は、データを見る前の信念(事前分布)とデータの当てはまり(尤度)を掛けて、データを見た後の信念(事後分布)を作る規則です。MAP 推定は事後分布の最頻値を取ることです。
  • これまで「過学習を防ぐおまじない」として導入してきた正則化は、事前分布そのものです。L2 正則化はガウス事前分布の、L1 正則化はラプラス事前分布の MAP 推定に一致し、正則化係数は λ=σ2/τ2\lambda = \sigma^2/\tau^2(観測ノイズの分散と事前分散の比)という意味を持ちます。
  • 事後分布を 1 点に潰さず、パラメータについて積分すると事後予測分布が得られます。その分散は「観測ノイズ」と「パラメータの不確実性」の和に分解され、外挿するほど後者が効きます。
  • MAP は事後分布のたった 1 点であり、座標変換で不変ではありません。L1 のスパース性も「最頻値の性質」であって「事後平均の性質」ではありません。

1. 動機:点で答えるか、分布で答えるか

Section titled “1. 動機:点で答えるか、分布で答えるか”

線形回帰と最小二乗法では、残差平方和 i(yiwTxi)2\sum_i (y_i - \boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i)^2 を最小にする w\boldsymbol{w} を求めました(最小二乗問題(Definition 3.1)[Linear Regression and Least Squares])。ロジスティック回帰では交差エントロピー(交差エントロピー誤差(Definition 4.1)[Logistic Regression])を最小にしました。そして過学習を抑えるために、どちらでも罰則項 λw2\lambda \lVert \boldsymbol{w}\rVert^2 を足しました。

ここで素朴な疑問が 3 つ残ります。

  1. なぜ残差の 2 乗なのですか。絶対値でも 4 乗でもよいはずです。
  2. なぜ罰則が w2\lVert \boldsymbol{w}\rVert^2 なのですか。λ\lambda はどうやって決めるのですか。単位は何ですか。
  3. モデルが返すのは y^=w^Tx\hat{y} = \hat{\boldsymbol{w}}^{\mathsf{T}}\boldsymbol{x} という 1 つの数だけです。「この予測はどのくらい信用してよいか」はどこに書いてあるのですか。

この章の主張は、3 つとも確率モデルを明示すれば同時に答えが出る、というものです。1 の答えは「観測ノイズをガウス分布と仮定したから」、2 の答えは「重みの事前分布をガウス分布と仮定したから、λ\lambda はノイズ分散と事前分散の比」、3 の答えは「点推定を捨てて事後分布のまま積分すればよい」です。

歴史的には、この見方は機械学習より 250 年ほど古いものです。トーマス・ベイズの遺稿(1763 年)と、それを独立に一般の形で展開したラプラスの 1774 年の論文は、「結果から原因の確率を測る」問題として条件付き確率の反転を扱いました。一方、20 世紀初頭の R. A. フィッシャーは事前分布を持ち込むことを嫌い、尤度だけを使う最尤推定を統計学の中心に据えました。現代の機械学習はこの両方を、目的に応じて使い分けています。損失関数の設計は前者の言葉で、汎化性能の議論は後者の言葉で語られることが多い、という具合です。

2. 準備:尤度・事前分布・ベイズの定理

Section titled “2. 準備:尤度・事前分布・ベイズの定理”

観測データを D\mathcal{D}、モデルのパラメータを θ\theta(ベクトルのときは w\boldsymbol{w})と書きます。確率モデルとは、θ\theta を決めるごとに D\mathcal{D} の生成される確率(密度)p(Dθ)p(\mathcal{D} \mid \theta) が定まる仕組みのことです。

Definition 2.1尤度・事前分布・事後分布

確率モデル p(Dθ)p(\mathcal{D} \mid \theta) において、データ D\mathcal{D} を固定し、θ\theta の関数と見たもの

L(θ)=p(Dθ)L(\theta) = p(\mathcal{D} \mid \theta)

尤度関数といい、その対数 (θ)=logL(θ)\ell(\theta) = \log L(\theta)対数尤度といいます。とくに D=(z1,,zN)\mathcal{D} = (z_1,\ldots,z_N) が独立同分布のとき、

L(θ)=i=1Np(ziθ),(θ)=i=1Nlogp(ziθ)L(\theta) = \prod_{i=1}^{N} p(z_i \mid \theta), \qquad \ell(\theta) = \sum_{i=1}^{N} \log p(z_i \mid \theta)

となります。

パラメータ θ\theta 自身を確率変数と見なし、データを観測する前の分布 p(θ)p(\theta)事前分布、観測後の条件付き分布 p(θD)p(\theta \mid \mathcal{D})事後分布といいます。また

p(D)=p(Dθ)p(θ)dθp(\mathcal{D}) = \int p(\mathcal{D} \mid \theta)\, p(\theta)\, d\theta

周辺尤度(またはエビデンス)といいます。θ\theta が離散なら積分は和に置き換えます。

尤度は「θ\theta についての確率分布」ではないことに注意してください。L(θ)dθ\int L(\theta)\,d\theta は一般に 11 になりません。LL はあくまで「この θ\theta を信じたとき、手元のデータはどれくらい出やすかったか」を測る物差しです。

Theorem 2.2ベイズの定理

θ\thetaD\mathcal{D} の同時密度 p(θ,D)p(\theta, \mathcal{D}) が存在し、p(D)>0p(\mathcal{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'}

が成り立つ。とくに D\mathcal{D} を固定すれば分母は θ\theta に依存しない定数なので、

p(θD)p(Dθ)p(θ)p(\theta \mid \mathcal{D}) \propto p(\mathcal{D} \mid \theta)\, p(\theta)

すなわち「事後分布 \propto 尤度 ×\times 事前分布」である。

Proof(Theorem 2.2)

条件付き密度の定義から、p(D)>0p(\mathcal{D}) > 0 のとき

p(θD)=p(θ,D)p(D)p(\theta \mid \mathcal{D}) = \frac{p(\theta, \mathcal{D})}{p(\mathcal{D})}

です。同じ定義を逆向きに使うと、p(θ)>0p(\theta) > 0 となる θ\theta について p(Dθ)=p(θ,D)/p(θ)p(\mathcal{D} \mid \theta) = p(\theta, \mathcal{D}) / p(\theta)、すなわち p(θ,D)=p(Dθ)p(θ)p(\theta, \mathcal{D}) = p(\mathcal{D} \mid \theta)\, p(\theta) です(p(θ)=0p(\theta) = 0 の点では両辺とも 00 なのでこの等式は全体で成り立ちます)。これを最初の式の分子に代入すれば第 1 の等号が得られます。

第 2 の等号は、Definition 2.1 の周辺尤度の定義そのものです。実際、同時密度を θ\theta について積分すれば D\mathcal{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'

となります。最後の比例関係は、分母が θ\theta を含まないことから直ちに従います。

比例関係のほうが実用上は重要です。事後分布を求める作業の大半は、「尤度 ×\times 事前分布」を θ\theta の関数として書き下し、θ\theta を含まない因子をすべて捨てて、残った形が何という分布かを見抜くことに尽きます。

Example 2.3事前分布が効く例:まれな病気の検査

有病率 0.1%0.1\% の病気があり、検査の感度(病気の人が陽性になる確率)が 99%99\%、特異度(健康な人が陰性になる確率)が 95%95\% だとします。ある人の検査が陽性でした。この人が病気である確率はいくらでしょうか。

θ{病気,健康}\theta \in \{\text{病気}, \text{健康}\}D=陽性\mathcal{D} = \text{陽性} とします。事前分布は p(病気)=0.001p(\text{病気}) = 0.001p(健康)=0.999p(\text{健康}) = 0.999。尤度は p(陽性病気)=0.99p(\text{陽性} \mid \text{病気}) = 0.99p(陽性健康)=10.95=0.05p(\text{陽性} \mid \text{健康}) = 1 - 0.95 = 0.05 です。Theorem 2.2 より

p(陽性)=0.99×0.001+0.05×0.999=0.00099+0.04995=0.05094,p(病気陽性)=0.000990.05094=0.019431.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}

検査の精度が高いにもかかわらず、陽性でも病気である確率は 2%2\% 程度です。尤度比は 0.99/0.05=19.80.99/0.05 = 19.8 倍ありますが、事前オッズが 1:9991:999 と極端に小さいので、事後オッズは 19.8/9991/50.519.8/999 \approx 1/50.5 にしかならないためです。

「尤度だけを見て θ\theta を選ぶ」(この場合は p(陽性θ)p(\text{陽性}\mid\theta) が大きい「病気」を選ぶ)という判断が、事前分布を無視したときにどれほど誤るかを示す例です。

Remark 2.4頻度論とベイズの分かれ目

θ\theta 自身に確率分布を置く」という一歩は、見た目より大きな一歩です。頻度論の立場では θ\theta は未知だが固定された定数であり、確率は繰り返し試行の相対頻度としてのみ意味を持つので、p(θ)p(\theta) は書けません。ベイズの立場では確率を「信念の度合い」と読み、未知の定数にも分布を置きます。

機械学習の実務では、この対立を教義として扱う必要はほとんどありません。むしろ「事前分布は、モデルに入れたい構造(重みは小さいはず、係数の多くは 00 のはず、関数は滑らかなはず)を確率の言葉で書く道具である」と考えると、この章の内容はすべて技術的な道具として使えます。

3. 最尤推定:データだけを見る

Section titled “3. 最尤推定:データだけを見る”

Definition 3.1最尤推定量

パラメータ空間 Θ\Theta 上の確率モデル p(Dθ)p(\mathcal{D} \mid \theta) に対し、尤度関数 L(θ)=p(Dθ)L(\theta) = p(\mathcal{D}\mid\theta) を最大にする θ\theta

θ^ML=arg maxθΘp(Dθ)=arg maxθΘlogp(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)

最尤推定量といいます。log\log は狭義単調増加なので、LL を最大にする点と =logL\ell = \log L を最大にする点は一致します。

対数を取る理由は 2 つあります。独立同分布の積が和に変わって微分しやすくなること、そして NN が大きいときに L(θ)L(\theta) がアンダーフローするほど小さくなるのを避けられることです。

Example 3.2コイン投げの最尤推定

表の出る確率が θ[0,1]\theta \in [0,1] のコインを nn 回投げ、kk 回表が出たとします。各回が独立なので

L(θ)=θk(1θ)nk,(θ)=klogθ+(nk)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).

微分して

(θ)=kθnk1θ=k(1θ)(nk)θθ(1θ)=knθθ(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)}.

分母は 0<θ<10 < \theta < 1 で正なので、(θ)=0    θ=k/n\ell'(\theta) = 0 \iff \theta = k/n です。さらに 0<k<n0 < k < n のとき

(θ)=kθ2nk(1θ)2<0\ell''(\theta) = -\frac{k}{\theta^{2}} - \frac{n-k}{(1-\theta)^{2}} < 0

なので \ell は狭義凹であり、θ=k/n\theta = k/n が唯一の最大点です。よって θ^ML=k/n\hat{\theta}_{\mathrm{ML}} = k/n、つまり単なる標本比率です。

端の場合を省略しないでおきます。k=nk = n のとき (θ)=nlogθ\ell(\theta) = n\log\theta(0,1)(0,1) で狭義単調増加なので、最大は端点 θ=1\theta = 1 で達成されます(L(1)=1L(1) = 1)。同様に k=0k = 0 なら θ^ML=0\hat{\theta}_{\mathrm{ML}} = 0 です。つまり3 回投げて 3 回表なら、最尤推定は「このコインは絶対に裏が出ない」と断言します。データと矛盾はしていませんが、賭けの根拠にはできません。この病理が次節の動機です。

次の命題は、線形回帰 で天下り的に採用した「残差の 2 乗和」が、実はガウスノイズの仮定と同じものであることを示します。

Proposition 3.3ガウス雑音の下で最尤推定は最小二乗法

入力 x1,,xNRd\boldsymbol{x}_1,\ldots,\boldsymbol{x}_N \in \mathbb{R}^{d} は固定された既知の値とし、特徴写像 ϕ:RdRM\boldsymbol{\phi} : \mathbb{R}^{d} \to \mathbb{R}^{M} を用いて計画行列 ΦRN×M\Phi \in \mathbb{R}^{N \times M} を、その第 ii 行が ϕ(xi)T\boldsymbol{\phi}(\boldsymbol{x}_i)^{\mathsf{T}} であるものとして定める。既知の定数 σ2>0\sigma^2 > 0 に対し、出力が

yi=wTϕ(xi)+εi,ε1,,εN は独立で εiN(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=(y1,,yN)T\boldsymbol{y} = (y_1,\ldots,y_N)^{\mathsf{T}} とおくと、wRM\boldsymbol{w} \in \mathbb{R}^{M} の最尤推定量は

w^ML=arg minwRMyΦw2\hat{\boldsymbol{w}}_{\mathrm{ML}} = \operatorname*{arg\,min}_{\boldsymbol{w} \in \mathbb{R}^{M}} \lVert \boldsymbol{y} - \Phi \boldsymbol{w} \rVert^{2}

の解全体と一致する。とくに Φ\Phi が列フルランクならば w^ML=(ΦTΦ)1ΦTy\hat{\boldsymbol{w}}_{\mathrm{ML}} = (\Phi^{\mathsf{T}}\Phi)^{-1}\Phi^{\mathsf{T}}\boldsymbol{y} でただ 1 つに定まる。

Proof(Proposition 3.3)

εiN(0,σ2)\varepsilon_i \sim \mathcal{N}(0,\sigma^2) より、w\boldsymbol{w} を与えたときの yiy_i の密度は

p(yiw)=12πσ2exp ⁣((yiwTϕ(xi))22σ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)

です。εi\varepsilon_i が独立なので y1,,yNy_1,\ldots,y_N も(w\boldsymbol{w} を与えたとき)独立で、Definition 2.1 より尤度は積になります。対数を取ると

(w)=i=1Nlogp(yiw)=N2log(2πσ2)    12σ2i=1N(yiwTϕ(xi))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}.

右辺第 1 項は w\boldsymbol{w} を含まない定数です。また i(yiwTϕ(xi))2=yΦw2\sum_i (y_i - \boldsymbol{w}^{\mathsf{T}}\boldsymbol{\phi}(\boldsymbol{x}_i))^{2} = \lVert \boldsymbol{y} - \Phi\boldsymbol{w}\rVert^{2} は、計画行列の定義から Φw\Phi\boldsymbol{w} の第 ii 成分が ϕ(xi)Tw\boldsymbol{\phi}(\boldsymbol{x}_i)^{\mathsf{T}}\boldsymbol{w} であることによります。したがって

(w)=const12σ2yΦw2.\ell(\boldsymbol{w}) = \mathrm{const} - \frac{1}{2\sigma^{2}} \lVert \boldsymbol{y} - \Phi\boldsymbol{w}\rVert^{2}.

係数 1/(2σ2)-1/(2\sigma^{2}) は負の定数なので、\ell を最大化することと yΦw2\lVert \boldsymbol{y} - \Phi\boldsymbol{w}\rVert^{2} を最小化することは同値です。

後半は正規方程式です。J(w)=yΦw2J(\boldsymbol{w}) = \lVert \boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} を展開すると J=yTy2yTΦw+wTΦTΦwJ = \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=2ΦTy+2ΦTΦw\nabla J = -2\Phi^{\mathsf{T}}\boldsymbol{y} + 2\Phi^{\mathsf{T}}\Phi\boldsymbol{w} です。JJ は凸(ヘッセ行列 2ΦTΦ2\Phi^{\mathsf{T}}\Phi が半正定値)なので J=0\nabla J = \boldsymbol{0} が最小の必要十分条件で、ΦTΦw=ΦTy\Phi^{\mathsf{T}}\Phi \boldsymbol{w} = \Phi^{\mathsf{T}}\boldsymbol{y} を得ます。Φ\Phi が列フルランクなら ΦTΦ\Phi^{\mathsf{T}}\Phi は正定値(v0\boldsymbol{v} \ne \boldsymbol{0} に対し vTΦTΦv=Φv2>0\boldsymbol{v}^{\mathsf{T}}\Phi^{\mathsf{T}}\Phi\boldsymbol{v} = \lVert \Phi\boldsymbol{v}\rVert^{2} > 0)で可逆なので、解は一意に定まります。

同じ論法で、ノイズをラプラス分布 p(ε)eε/cp(\varepsilon) \propto e^{-|\varepsilon|/c} とすれば最尤推定は残差の絶対値和の最小化(最小絶対偏差回帰)になります。「なぜ 2 乗か」の答えは「ガウスを仮定したから」であり、逆にいえば外れ値が多いデータでは 2 乗以外を選ぶ理由がある、ということです。

4. 最大事後確率推定:事前分布を足す

Section titled “4. 最大事後確率推定:事前分布を足す”

Definition 4.1MAP 推定量

事前分布 p(θ)p(\theta) を持つモデルにおいて、事後分布 p(θD)p(\theta \mid \mathcal{D}) を最大にする点

θ^MAP=arg maxθΘp(θD)\hat{\theta}_{\mathrm{MAP}} = \operatorname*{arg\,max}_{\theta \in \Theta} p(\theta \mid \mathcal{D})

最大事後確率推定量(MAP 推定量)といいます。

Proposition 4.2MAP の目的関数

p(D)>0p(\mathcal{D}) > 0 のとき、

θ^MAP=arg maxθΘ{logp(Dθ)対数尤度+logp(θ)対数事前分布}\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\}

が成り立つ。とくに事前分布が Θ\Theta 上で定数(一様)ならば θ^MAP=θ^ML\hat{\theta}_{\mathrm{MAP}} = \hat{\theta}_{\mathrm{ML}} である。

Proof(Proposition 4.2)

Theorem 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(\mathcal{D})θ\theta に依存しない正の定数なので、p(θD)p(\theta\mid\mathcal{D}) を最大にする θ\thetap(Dθ)p(θ)p(\mathcal{D}\mid\theta)p(\theta) を最大にする θ\theta は同じです。対数は狭義単調増加なので、さらに log(p(Dθ)p(θ))=logp(Dθ)+logp(θ)\log\bigl(p(\mathcal{D}\mid\theta)p(\theta)\bigr) = \log p(\mathcal{D}\mid\theta) + \log p(\theta) を最大にする θ\theta とも同じです。

事前分布が定数 c>0c > 0 なら logp(θ)=logc\log p(\theta) = \log c も定数なので、目的関数は対数尤度と定数だけ違い、最大点は Definition 3.1 の最尤推定量に一致します。

この命題が、この章でいちばん使う道具です。MAP 推定は「対数尤度 ++ 対数事前分布」の最大化であり、機械学習の言葉に翻訳すれば「損失関数 ++ 罰則項」の最小化にほかなりません。第 5 節でこの翻訳を厳密に実行します。

Definition 4.3共役事前分布

尤度の族 {p(Dθ)}θΘ\{p(\mathcal{D}\mid\theta)\}_{\theta\in\Theta} に対し、事前分布の族 P\mathcal{P}共役であるとは、任意の p(θ)Pp(\theta) \in \mathcal{P} と任意の観測 D\mathcal{D} に対して事後分布 p(θD)p(\theta\mid\mathcal{D}) もまた P\mathcal{P} に属することをいいます。

共役性は数学的な必然ではなく、計算の便宜です。事後分布が同じ族に留まってくれれば、積分を実行せずに「パラメータの更新則」だけで事後分布が書けます。

Example 4.4ベータ事前分布とコイン投げ

Example 3.2 と同じ設定で、事前分布としてベータ分布 Beta(a,b)\mathrm{Beta}(a,b)a,b>0a,b > 0

p(θ)=1B(a,b)θa1(1θ)b1,0<θ<1p(\theta) = \frac{1}{B(a,b)}\, \theta^{a-1}(1-\theta)^{b-1}, \qquad 0 < \theta < 1

を取ります。B(a,b)B(a,b)θ\theta に依存しない正規化定数です。Theorem 2.2 の比例形を使うと

p(θD)    θk(1θ)nk尤度θa1(1θ)b1事前分布=θ(k+a)1(1θ)(nk+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}.

右辺は Beta(k+a,nk+b)\mathrm{Beta}(k+a,\, n-k+b) の密度の θ\theta 依存部分そのものです。密度は正規化定数まで込めて一意なので、事後分布は Beta(k+a,nk+b)\mathrm{Beta}(k+a,\, n-k+b) です。つまりベータ分布族はベルヌーイ/二項尤度の共役事前分布であり、更新則は「aa に表の回数を、bb に裏の回数を足す」だけです。a,ba, b疑似観測回数と読めます。

具体的な数値を最後まで追います。n=3n = 3k=3k = 3(3 回投げて 3 回とも表)、事前分布は Beta(2,2)\mathrm{Beta}(2,2)(「0.50.5 のあたりが怪しい」という穏やかな信念、表裏 1 回ずつを見たのと同じ重み)とします。事後分布は Beta(5,2)\mathrm{Beta}(5, 2) で、密度は

p(θD)=θ4(1θ)B(5,2),B(5,2)=Γ(5)Γ(2)Γ(7)=4!1!6!=24720=130,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)=30θ4(1θ)p(\theta\mid\mathcal{D}) = 30\,\theta^{4}(1-\theta) です。最頻値は

ddθ[4logθ+log(1θ)]=4θ11θ=45θθ(1θ)=0    θ=45,\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},

二階微分は 4/θ21/(1θ)2<0-4/\theta^{2} - 1/(1-\theta)^{2} < 0 なのでこれが唯一の最大点です。よって θ^MAP=4/5=0.8\hat{\theta}_{\mathrm{MAP}} = 4/5 = 0.8。一方、事後平均は Beta(α,β)\mathrm{Beta}(\alpha,\beta) の平均公式 α/(α+β)\alpha/(\alpha+\beta) より 5/70.7145/7 \approx 0.714 です。

まとめると、同じデータに対して

θ^ML=1,θ^MAP=0.8,E[θD]=570.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

という 3 つの答えが出ます。最尤推定の「絶対に裏は出ない」という断言は消え、しかも MAP と事後平均も一致しません。3 つ目の値は「次の 1 回が表である確率」でもあります(第 6 節の事後予測分布)。

00.51MAP 0.8MLE 1事後平均 5/7事前分布 Beta(2,2)尤度(規格化)事後分布 Beta(5,2)
3 回投げて 3 回表だったときの事前分布 Beta(2,2)、尤度、事後分布 Beta(5,2)。尤度は θ の密度ではないので高さは見やすいように規格化してあります。

図で見ると、事後分布が「事前分布と尤度の綱引きの結果」であることがよくわかります。尤度は θ=1\theta = 1 に向かって単調に増えていますが、事前分布が θ=1\theta = 100 に落ちるので、積は 0.80.8 付近で山を作ります。データが増えれば尤度の山が鋭くなり、事前分布の影響は相対的に小さくなっていきます。

Remark 4.5MAP は座標変換で不変ではない

MAP は「いちばんありそうな値」という直感的な説明をされますが、これは密度の最頻値であって、パラメータの取り方に依存します。

Example 4.4 の事後分布 p(θD)=30θ4(1θ)p(\theta\mid\mathcal{D}) = 30\theta^{4}(1-\theta) で、パラメータを η=θ2\eta = \theta^{2} に取り替えてみます。θ=η\theta = \sqrt{\eta}dθ/dη=1/(2η)d\theta/d\eta = 1/(2\sqrt{\eta}) なので、変数変換の公式より η\eta の密度は

g(η)=p(ηD)12η=30η2(1η)12η=15(η3/2η2).g(\eta) = p(\sqrt{\eta}\mid\mathcal{D}) \cdot \frac{1}{2\sqrt{\eta}} = 30\,\eta^{2}\,(1-\sqrt{\eta}) \cdot \frac{1}{2\sqrt{\eta}} = 15\bigl(\eta^{3/2} - \eta^{2}\bigr).

その最頻値は g(η)=15(32η1/22η)=0g'(\eta) = 15\bigl(\tfrac{3}{2}\eta^{1/2} - 2\eta\bigr) = 0 から、η>0\eta > 0η1/2\eta^{1/2} で割って 32=2η1/2\tfrac{3}{2} = 2\eta^{1/2}、すなわち η=9/16\eta = 9/16 です。これを θ\theta に戻すと 9/16=3/4=0.75\sqrt{9/16} = 3/4 = 0.75 で、θ\theta 座標での MAP 0.80.8 とは一致しません

一方、事後分布そのものは変数変換の公式で正しく移り変わりますし、事後平均も「何を損失関数と考えるか」を決めれば意味が定まります(2 乗損失の下でのベイズ推定量が事後平均です)。MAP は計算が軽い代わりに、この種の恣意性を抱えていることを覚えておいてください。

5. 正則化の正体:罰則項は事前分布である

Section titled “5. 正則化の正体:罰則項は事前分布である”

ここが本章の中心です。Proposition 4.2 の「対数尤度 ++ 対数事前分布」に、Proposition 3.3 のガウス雑音モデルとガウス事前分布を代入するだけで、リッジ回帰が出てきます。

Theorem 5.1L2 正則化はガウス事前分布の MAP 推定

Proposition 3.3 と同じ設定(入力は固定、ΦRN×M\Phi \in \mathbb{R}^{N\times M} は計画行列、σ2>0\sigma^2 > 0 は既知の雑音分散、εi\varepsilon_i は独立に N(0,σ2)\mathcal{N}(0,\sigma^2))に加えて、重みの事前分布を

wN(0,τ2IM),τ2>0 は既知\boldsymbol{w} \sim \mathcal{N}(\boldsymbol{0},\, \tau^{2} I_M), \qquad \tau^{2} > 0 \ \text{は既知}

とし、w\boldsymbol{w}(ε1,,εN)(\varepsilon_1,\ldots,\varepsilon_N) は独立とする。このとき MAP 推定量は

w^MAP=arg minwRM{yΦw2+λw2},λ=σ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}}

で与えられ、これはリッジ回帰(L2 正則化つき最小二乗法)の解に一致する。さらにこの最小化問題の解は一意で、

w^MAP=(ΦTΦ+λIM)1ΦTy\hat{\boldsymbol{w}}_{\mathrm{MAP}} = \bigl(\Phi^{\mathsf{T}}\Phi + \lambda I_M\bigr)^{-1} \Phi^{\mathsf{T}} \boldsymbol{y}

である。

Proof(Theorem 5.1)

第 1 段:目的関数を書き下す。 Proposition 3.3 の証明で計算したとおり、対数尤度は

logp(yw)=N2log(2πσ2)12σ2yΦw2\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}

です(w\boldsymbol{w}εi\varepsilon_i の独立性から、w\boldsymbol{w} を与えたときの y\boldsymbol{y} の条件付き分布は元のモデルのままであることを使いました)。また MM 次元ガウス事前分布 N(0,τ2IM)\mathcal{N}(\boldsymbol{0},\tau^{2}I_M) の密度は p(w)=(2πτ2)M/2exp(w2/(2τ2))p(\boldsymbol{w}) = (2\pi\tau^{2})^{-M/2}\exp\bigl(-\lVert\boldsymbol{w}\rVert^{2}/(2\tau^{2})\bigr) なので

logp(w)=M2log(2πτ2)12τ2w2.\log p(\boldsymbol{w}) = -\frac{M}{2}\log(2\pi\tau^{2}) - \frac{1}{2\tau^{2}}\lVert \boldsymbol{w}\rVert^{2}.

第 2 段:定数を落として整理する。 Proposition 4.2 より w^MAP\hat{\boldsymbol{w}}_{\mathrm{MAP}}

logp(yw)+logp(w)=const12σ2yΦw212τ2w2\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}

の最大点です。const\mathrm{const}w\boldsymbol{w} を含まないので落とせます。符号を反転すると最小化問題

w^MAP=arg minw{12σ2yΦw2+12τ2w2}\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\}

になります。目的関数全体を正の定数 2σ22\sigma^{2} 倍しても最小点は変わらないので

w^MAP=arg minw{yΦw2+σ2τ2w2}\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\}

を得ます。λ=σ2/τ2\lambda = \sigma^{2}/\tau^{2} とおけば主張の形です。

第 3 段:閉じた形の解。 J(w)=yΦw2+λw2J(\boldsymbol{w}) = \lVert \boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} + \lambda\lVert\boldsymbol{w}\rVert^{2} を展開すると

J(w)=yTy2yTΦw+wT(ΦTΦ+λIM)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)=2ΦTy+2(ΦTΦ+λIM)w\nabla J(\boldsymbol{w}) = -2\Phi^{\mathsf{T}}\boldsymbol{y} + 2(\Phi^{\mathsf{T}}\Phi+\lambda I_M)\boldsymbol{w} です。ここで A=ΦTΦ+λIMA = \Phi^{\mathsf{T}}\Phi + \lambda I_M は正定値です。実際、v0\boldsymbol{v}\ne\boldsymbol{0} に対し

vTAv=Φv2+λv2λv2>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

となります(λ=σ2/τ2>0\lambda = \sigma^2/\tau^2 > 0 を使いました)。したがって JJ のヘッセ行列 2A2A は正定値で JJ は狭義凸、J(w)=0\nabla J(\boldsymbol{w}) = \boldsymbol{0} すなわち Aw=ΦTyA\boldsymbol{w} = \Phi^{\mathsf{T}}\boldsymbol{y} が唯一の最小点を与えます。AA は正定値ゆえ可逆なので w^MAP=A1ΦTy\hat{\boldsymbol{w}}_{\mathrm{MAP}} = A^{-1}\Phi^{\mathsf{T}}\boldsymbol{y} です。

λ=σ2/τ2\lambda = \sigma^{2}/\tau^{2} という表式は、正則化係数に明確な意味を与えます。λ\lambda は「観測がどれだけ雑か」と「重みがどれだけ大きくてよいか」の比です。雑音が大きい(σ2\sigma^2 大)ほどデータを信用せず罰則を強め、事前の許容幅が広い(τ2\tau^2 大)ほど罰則を緩めます。τ2\tau^{2}\to\infty(何も知らない)とすれば λ0\lambda\to 0 で最尤推定に戻り、これは Proposition 4.2 の後半(一様事前分布なら MLE)とも整合します。

Example 5.2リッジ回帰は固有値の小さい方向を強く縮める

ΦTΦ\Phi^{\mathsf{T}}\Phi は対称半正定値なので、スペクトル定理Corollary 4.3[スペクトル定理])により直交行列 UU と対角行列 Λ=diag(d1,,dM)\Lambda = \operatorname{diag}(d_1,\ldots,d_M)dj0d_j \ge 0)を用いて ΦTΦ=UΛUT\Phi^{\mathsf{T}}\Phi = U\Lambda U^{\mathsf{T}} と書けます。UTU=IU^{\mathsf{T}}U = I より ΦTΦ+λI=U(Λ+λI)UT\Phi^{\mathsf{T}}\Phi + \lambda I = U(\Lambda+\lambda I)U^{\mathsf{T}} なので、Theorem 5.1 の解は

w^MAP=U(Λ+λI)1UTΦTy.\hat{\boldsymbol{w}}_{\mathrm{MAP}} = U(\Lambda+\lambda I)^{-1}U^{\mathsf{T}}\Phi^{\mathsf{T}}\boldsymbol{y}.

z=UTΦTy\boldsymbol{z} = U^{\mathsf{T}}\Phi^{\mathsf{T}}\boldsymbol{y} とおき、固有ベクトル基底での成分を比べます。dj>0d_j > 0 のとき最尤解の第 jj 成分は zj/djz_j/d_j、リッジ解の第 jj 成分は zj/(dj+λ)z_j/(d_j+\lambda) なので、その比は

(UTw^MAP)j(UTw^ML)j=djdj+λ.\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}.

たとえば λ=1\lambda = 1d1=100d_1 = 100d2=0.01d_2 = 0.01 なら、縮小率は第 1 方向で 100/1010.990100/101 \approx 0.990、第 2 方向で 0.01/1.010.00990.01/1.01 \approx 0.0099 です。データの分散が大きい方向はほとんど手つかず、分散がほとんどない方向は 1/1001/100 に潰されます

djd_jΦTΦ\Phi^{\mathsf{T}}\Phi の固有値、つまり主成分分析でいうところの各主成分方向(Definition 3.2[主成分分析])のデータの散らばりです。「データが何も語っていない方向では事前分布が勝つ」という Theorem 5.1 のベイズ的な読みが、そのまま数式に現れています。dj=0d_j = 0(その方向にはデータが皆無)なら成分は 00、すなわち事前分布の平均そのものになります。

同じ計算をラプラス事前分布で行うと L1 正則化(ラッソ)が出ます。

Proposition 5.3L1 正則化はラプラス事前分布の MAP 推定

Theorem 5.1 と同じ雑音モデルの下で、重みの事前分布を、各成分が独立に平均 00・尺度 b>0b > 0 のラプラス分布

p(w)=j=1M12bexp ⁣(wjb)=1(2b)Mexp ⁣(w1b)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)

に従うものとする(w1=jwj\lVert\boldsymbol{w}\rVert_1 = \sum_j |w_j|)。このとき

w^MAP=arg minwRM{yΦw2+λw1},λ=2σ2b\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}

であり、これはラッソ(L1 正則化つき最小二乗法)の解に一致する。

Proof(Proposition 5.3)

対数事前分布は

logp(w)=Mlog(2b)1bw1\log p(\boldsymbol{w}) = -M\log(2b) - \frac{1}{b}\lVert\boldsymbol{w}\rVert_{1}

です。Proposition 4.2Theorem 5.1 の証明第 1 段の対数尤度を合わせると、w^MAP\hat{\boldsymbol{w}}_{\mathrm{MAP}}

const12σ2yΦw21bw1\mathrm{const} - \frac{1}{2\sigma^{2}}\lVert\boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} - \frac{1}{b}\lVert\boldsymbol{w}\rVert_{1}

の最大点です。定数を落とし、符号を反転し、全体を正の定数 2σ22\sigma^{2} 倍すると

w^MAP=arg minw{yΦw2+2σ2bw1}\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\}

を得ます。λ=2σ2/b\lambda = 2\sigma^{2}/b とおけば主張の形です。

なお目的関数は凸(2 乗項は凸、1\lVert\cdot\rVert_1 はノルムなので凸)ですが 1\lVert\cdot\rVert_1wj=0w_j = 0 で微分可能でないため、Theorem 5.1 のような閉じた形の解は一般には得られません。ΦTΦ=I\Phi^{\mathsf{T}}\Phi = I(列が正規直交)の特別な場合は成分ごとに分離でき、解はソフトしきい値関数で書けます(Exercise 7.3)。

対応関係を表にまとめます。

罰則項対応する事前分布事前密度(θ\theta 依存部分)正則化係数解の特徴
なし一様(非正則な事前分布)1\propto 1λ=0\lambda = 0最尤推定に一致
λw22\lambda \lVert \boldsymbol{w}\rVert_{2}^{2}ガウス N(0,τ2I)\mathcal{N}(\boldsymbol{0},\tau^{2}I)exp(w22/(2τ2))\exp\bigl(-\lVert\boldsymbol{w}\rVert_2^{2}/(2\tau^{2})\bigr)λ=σ2/τ2\lambda = \sigma^{2}/\tau^{2}全成分を滑らかに縮小
λw1\lambda \lVert \boldsymbol{w}\rVert_{1}ラプラス(尺度 bbexp(w1/b)\exp\bigl(-\lVert\boldsymbol{w}\rVert_1/b\bigr)λ=2σ2/b\lambda = 2\sigma^{2}/b一部の成分が厳密に 00
λwμ22\lambda\lVert\boldsymbol{w}-\boldsymbol{\mu}\rVert_2^{2}ガウス N(μ,τ2I)\mathcal{N}(\boldsymbol{\mu},\tau^{2}I)exp(wμ22/(2τ2))\exp\bigl(-\lVert\boldsymbol{w}-\boldsymbol{\mu}\rVert_2^{2}/(2\tau^{2})\bigr)λ=σ2/τ2\lambda = \sigma^{2}/\tau^{2}既知の値 μ\boldsymbol{\mu} に引き寄せる

4 行目は、事前学習済みモデルからの微調整で「元の重みから離れすぎない」ようにする正則化が、そのまま「元の重みを平均とするガウス事前分布」であることを示しています。正則化を発明するとは、事前分布を設計することです。

Remark 5.4スパース性は「最頻値」の性質であって「事後分布」の性質ではない

L1 正則化が厳密な 00 を生むのは、ラプラス密度が原点で尖っている(微分不可能な角を持つ)ためです。ところがこれは事後分布の最頻値についての話です。ラプラス事前分布の下でも、事後分布 p(wD)p(\boldsymbol{w}\mid\mathcal{D}) は連続分布なので Pr(wj=0D)=0\Pr(w_j = 0 \mid \mathcal{D}) = 0 であり、事後平均 E[wjD]\mathbb{E}[w_j\mid\mathcal{D}] が厳密に 00 になることもまずありません。

つまり「ラッソはベイズ的にはラプラス事前分布の MAP である」は正しい一方、「ラッソの変数選択はベイズ的な変数選択である」とは言えません。本当に「その変数が不要である確率」を扱いたいなら、wj=0w_j = 0 に正の確率質量を置くスパイク・アンド・スラブ型の事前分布が必要になります。Remark 4.5 と合わせて、MAP は事後分布の要約としてかなり乱暴なものだと理解しておいてください。

6. 不確実性を扱う:ベイズ線形回帰と事後予測分布

Section titled “6. 不確実性を扱う:ベイズ線形回帰と事後予測分布”

MLE も MAP も、最後に θ\theta を 1 点に潰します。潰さずに事後分布のまま持ち歩くとどうなるかを見ます。線形回帰+ガウス事前分布は、この計算が手で最後まで実行できる数少ない例です。

Theorem 6.1ベイズ線形回帰の事後分布

Proposition 3.3 の雑音モデル(入力は固定、ΦRN×M\Phi\in\mathbb{R}^{N\times M}σ2>0\sigma^{2}>0 は既知)の下で、w\boldsymbol{w} の事前分布を N(m0,S0)\mathcal{N}(\boldsymbol{m}_0, S_0)S0S_0MM 次の対称正定値行列、m0RM\boldsymbol{m}_0\in\mathbb{R}^M)とし、w\boldsymbol{w} と雑音は独立とする。このとき事後分布は再びガウス分布

p(wy)=N(wmN,SN)p(\boldsymbol{w}\mid\boldsymbol{y}) = \mathcal{N}(\boldsymbol{w} \mid \boldsymbol{m}_N, S_N)

であり、

SN=(S01+1σ2ΦTΦ)1,mN=SN(S01m0+1σ2ΦTy)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)

で与えられる。とくにガウス分布は共役事前分布である(Definition 4.3)。

Proof(Theorem 6.1)

第 1 段:逆行列の存在。 A=S01+σ2ΦTΦA = S_0^{-1} + \sigma^{-2}\Phi^{\mathsf{T}}\Phi とおきます。S0S_0 が対称正定値なら S01S_0^{-1} も対称正定値です。また任意の v\boldsymbol{v} に対し vTΦTΦv=Φv20\boldsymbol{v}^{\mathsf{T}}\Phi^{\mathsf{T}}\Phi\boldsymbol{v} = \lVert\Phi\boldsymbol{v}\rVert^{2}\ge 0 なので ΦTΦ\Phi^{\mathsf{T}}\Phi は半正定値です。よって v0\boldsymbol{v}\ne\boldsymbol{0} のとき

vTAv=vTS01v+1σ2Φv2>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

であり AA は正定値、したがって可逆です。SN=A1S_N = A^{-1} と書きます。

第 2 段:指数部を w\boldsymbol{w} の 2 次式に整理する。 Theorem 2.2 の比例形より

p(wy)p(yw)p(w)exp(12σ2yΦw212(wm0)TS01(wm0))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)

です(w\boldsymbol{w} を含まない正規化定数はすべて比例記号に吸収しました)。指数の中を展開します。

yΦw2=yTy2yTΦw+wTΦ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},(wm0)TS01(wm0)=wTS01w2m0TS01w+m0TS01m0(\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

S01S_0^{-1} の対称性から交差項 2 つが等しいことを使いました)。w\boldsymbol{w} を含まない項をまた比例記号に吸収すると、指数は

12wT(S01+1σ2ΦTΦ)=Aw  +  (S01m0+1σ2ΦTy)T= bTw  +  const-\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}

の形になります。

第 3 段:ガウス分布と同定する。 第 2 段より、事後密度はある定数 C>0C > 0 を用いて p(wy)=Cexp(12wTAw+bTw)p(\boldsymbol{w}\mid\boldsymbol{y}) = C\exp\bigl(-\tfrac12\boldsymbol{w}^{\mathsf{T}}A\boldsymbol{w} + \boldsymbol{b}^{\mathsf{T}}\boldsymbol{w}\bigr) と書けます。AA は第 1 段より対称正定値なので、Lemma 7.5 が適用でき、この密度は N(A1b,A1)\mathcal{N}(A^{-1}\boldsymbol{b},\, A^{-1}) の密度に一致します。A1=SNA^{-1}=S_NA1b=SN(S01m0+σ2ΦTy)=mNA^{-1}\boldsymbol{b} = S_N(S_0^{-1}\boldsymbol{m}_0 + \sigma^{-2}\Phi^{\mathsf{T}}\boldsymbol{y}) = \boldsymbol{m}_N なので主張が従います。

Corollary 6.2ガウス事前分布の下では MAP は事後平均であり、リッジ解に一致する

Theorem 6.1 の設定で m0=0\boldsymbol{m}_0 = \boldsymbol{0}S0=τ2IMS_0 = \tau^{2}I_M とすると、

w^MAP=E[wy]=mN=(ΦTΦ+λIM)1ΦTy,λ=σ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}}

が成り立つ。

Proof(Corollary 6.2)

ガウス分布の密度 exp(12(wmN)TSN1(wmN))\propto \exp\bigl(-\tfrac12(\boldsymbol{w}-\boldsymbol{m}_N)^{\mathsf{T}}S_N^{-1}(\boldsymbol{w}-\boldsymbol{m}_N)\bigr) は、指数が w=mN\boldsymbol{w}=\boldsymbol{m}_N で最大値 00 を取り、SN1S_N^{-1} が正定値なので他の点では真に負です。よって最頻値は平均 mN\boldsymbol{m}_N に一致します。

具体形は代入するだけです。S01=τ2IS_0^{-1} = \tau^{-2}I より

SN=(1τ2I+1σ2ΦTΦ)1=σ2(σ2τ2I+Φ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},mN=SN(0+1σ2ΦTy)=σ2(ΦTΦ+λI)11σ2ΦTy=(ΦTΦ+λI)1ΦTy.\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}.

これは Theorem 5.1 で別の道筋(目的関数の最小化)から得た解と一致しています。

事後共分散 SN=σ2(ΦTΦ+λI)1S_N = \sigma^{2}(\Phi^{\mathsf{T}}\Phi+\lambda I)^{-1} は、点推定だけでは決して得られない情報です。これを使って予測に不確実性を持たせます。

Definition 6.3事後予測分布

新しい入力 x\boldsymbol{x}_{*} に対する出力 yy_{*}事後予測分布とは、パラメータを事後分布で積分消去した分布

p(yx,D)=p(yx,w)p(wD)dwp(y_{*}\mid \boldsymbol{x}_{*}, \mathcal{D}) = \int p(y_{*}\mid \boldsymbol{x}_{*}, \boldsymbol{w})\, p(\boldsymbol{w}\mid\mathcal{D})\, d\boldsymbol{w}

のことです。w^\hat{\boldsymbol{w}} を 1 つ選んで p(yx,w^)p(y_{*}\mid\boldsymbol{x}_{*},\hat{\boldsymbol{w}}) を使うプラグイン予測と対比されます。

Theorem 6.4ベイズ線形回帰の事後予測分布

Theorem 6.1 の設定の下で、新しい入力 x\boldsymbol{x}_{*} の特徴ベクトルを ϕ=ϕ(x)\boldsymbol{\phi}_{*} = \boldsymbol{\phi}(\boldsymbol{x}_{*}) とし、y=wTϕ+εy_{*} = \boldsymbol{w}^{\mathsf{T}}\boldsymbol{\phi}_{*} + \varepsilon_{*}εN(0,σ2)\varepsilon_{*}\sim\mathcal{N}(0,\sigma^{2})w\boldsymbol{w} および既存の雑音と独立とする。このとき

p(yx,D)=N(ymNTϕ,  σ2+ϕTSNϕ)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)

である。すなわち予測分散は

σ2観測ノイズ(偶然的不確実性)  +  ϕTSNϕパラメータの不確実性(認識的不確実性)\underbrace{\sigma^{2}}_{\text{観測ノイズ(偶然的不確実性)}} \;+\; \underbrace{\boldsymbol{\phi}_{*}^{\mathsf{T}}S_N\boldsymbol{\phi}_{*}}_{\text{パラメータの不確実性(認識的不確実性)}}

と分解される。

Proof(Theorem 6.4)

D\mathcal{D} を与えたとき、Theorem 6.1 より wN(mN,SN)\boldsymbol{w}\sim\mathcal{N}(\boldsymbol{m}_N,S_N) です。多変量ガウス分布の線形像はガウス分布です。実際、定ベクトル a\boldsymbol{a} に対する特性関数は

E[eitaTw]=exp ⁣(itaTmNt22aTSNa)\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)

であり、これは N(aTmN,aTSNa)\mathcal{N}(\boldsymbol{a}^{\mathsf{T}}\boldsymbol{m}_N,\, \boldsymbol{a}^{\mathsf{T}}S_N\boldsymbol{a}) の特性関数そのものです。a=ϕ\boldsymbol{a}=\boldsymbol{\phi}_{*} と取って

wTϕD    N(mNTϕ,  ϕTSNϕ).\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).

仮定より εN(0,σ2)\varepsilon_{*}\sim\mathcal{N}(0,\sigma^{2}) はこれと独立なので、独立なガウス確率変数の和は平均と分散をそれぞれ足したガウス分布に従います(特性関数が積になることから直ちに従います)。よって

y=wTϕ+εD    N(mNTϕ,  σ2+ϕTSNϕ)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)

です。これは Definition 6.3 の積分を実行したものにほかなりません。

この分解が実務上いちばん重要な結論です。σ2\sigma^{2} はデータをいくら集めても減りません(コインの表裏は永遠に当てられません)。一方 ϕTSNϕ\boldsymbol{\phi}_{*}^{\mathsf{T}}S_N\boldsymbol{\phi}_{*} は、SN1=S01+σ2ΦTΦS_N^{-1} = S_0^{-1}+\sigma^{-2}\Phi^{\mathsf{T}}\Phi がデータとともに増えるので NN が増えれば減っていきます。前者を偶然的(アレアトリック)不確実性、後者を認識的(エピステミック)不確実性と呼び、後者だけが「もっとデータを集めれば解消する」種類の不確実性です。能動学習でどのデータにラベルを付けるかを選ぶとき、指標にするのは後者です。

Example 6.5数値例:外挿すると分散が増える

M=1M=1ϕ(x)=x\boldsymbol{\phi}(x)=x(原点を通る直線 y=wxy = wx)とし、σ2=1\sigma^{2}=1、事前分布は wN(0,τ2)w\sim\mathcal{N}(0,\tau^{2})τ2=1\tau^{2}=1 とします。データは (x1,y1)=(1,2)(x_1,y_1)=(1,2)(x2,y2)=(2,3)(x_2,y_2)=(2,3) の 2 点です。

計画行列は Φ=(12)\Phi = \begin{pmatrix}1\\2\end{pmatrix}y=(23)\boldsymbol{y}=\begin{pmatrix}2\\3\end{pmatrix} なので

ΦTΦ=12+22=5,ΦTy=12+23=8.\Phi^{\mathsf{T}}\Phi = 1^{2}+2^{2} = 5, \qquad \Phi^{\mathsf{T}}\boldsymbol{y} = 1\cdot 2 + 2\cdot 3 = 8.

最尤推定は Proposition 3.3 より w^ML=8/5=1.6\hat{w}_{\mathrm{ML}} = 8/5 = 1.6 です。事後分布は Theorem 6.1 より

SN1=1τ2+1σ2ΦTΦ=1+5=6,SN=16,S_N^{-1} = \frac{1}{\tau^{2}} + \frac{1}{\sigma^{2}}\Phi^{\mathsf{T}}\Phi = 1 + 5 = 6, \qquad S_N = \frac{1}{6},mN=SN(0+1σ2ΦTy)=168=431.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.

検算として Corollary 6.2 を使うと、λ=σ2/τ2=1\lambda = \sigma^{2}/\tau^{2}=1(5+1)18=4/3(5+1)^{-1}\cdot 8 = 4/3 となり一致します。

事後予測分布を Theorem 6.4 で計算します(ϕ=x\boldsymbol{\phi}_{*}=x_{*} なので ϕTSNϕ=x2/6\boldsymbol{\phi}_{*}^{\mathsf{T}}S_N\boldsymbol{\phi}_{*} = x_{*}^{2}/6)。

xx_{*}予測平均 43x\tfrac{4}{3}x_{*}認識的分散 x2/6x_{*}^{2}/6全分散標準偏差95% 予測区間
111.3331.3330.1670.1671.1671.1671.0801.080[0.78, 3.45][-0.78,\ 3.45]
334.0004.0001.5001.5002.5002.5001.5811.581[0.90, 7.10][0.90,\ 7.10]
101013.33313.33316.66716.66717.66717.6674.2034.203[5.09, 21.57][5.09,\ 21.57]

95%95\% 区間は平均 ±1.96×\pm 1.96 \times 標準偏差で計算しました。たとえば x=10x_{*}=10 では 13.333±1.96×4.203=13.333±8.23813.333 \pm 1.96\times 4.203 = 13.333 \pm 8.238 です。)

もし点推定 w^=4/3\hat{w}=4/3 を使ったプラグイン予測をしていたら、どの xx_{*} でも分散は σ2=1\sigma^{2}=1、区間幅は ±1.96\pm 1.96 の一定です。x=10x_{*}=10 での真の予測区間は幅 ±8.24\pm 8.24 ですから、プラグイン予測は外挿領域で 4 倍以上も自信過剰になっています。データが x{1,2}x\in\{1,2\} 付近にしかないのに x=10x=10 を当てにいっている、という事実は事後共分散 SNS_N の中にしか書かれていません。

flowchart TD
D["観測データ"] --> L["尤度:θ ごとのデータの出やすさ"]
P["事前分布:データを見る前の θ の分布"] --> B["ベイズの定理"]
L --> B
L -->|最大化| MLE["最尤推定:θ を 1 点に決める"]
B --> Post["事後分布:データを見た後の θ の分布"]
Post -->|最頻値| MAP["MAP 推定:θ を 1 点に決める(=正則化つき学習)"]
Post -->|θ について積分| Pred["事後予測分布:予測を分布として返す"]
尤度・事前分布から、最尤推定・MAP 推定・事後予測分布へ至る道筋

Remark 6.6現実のモデルではどうするか

Theorem 6.1 が手で解けたのは、尤度がガウス、事前分布もガウス、モデルがパラメータについて線形、という三拍子が揃っていたからです。ニューラルネットワークでは事後分布 p(wD)p(\boldsymbol{w}\mid\mathcal{D}) は数百万次元の、正規化定数すら計算できない分布になります。実務では次の近似が使われます。

  • ラプラス近似w^MAP\hat{\boldsymbol{w}}_{\mathrm{MAP}} のまわりで対数事後分布を 2 次まで展開し、共分散をヘッセ行列の逆行列で近似する。Theorem 6.1 はこの近似が厳密に正しくなる場合にあたります。
  • マルコフ連鎖モンテカルロ(MCMC):事後分布からのサンプル列を作り、期待値をサンプル平均で置き換える。厳密だが計算量が大きい。
  • 変分推論:扱いやすい分布族の中で事後分布に最も近いものを最適化で探す。目的関数(ELBO)の最大化に帰着するので勾配降下法Definition 4.1[勾配降下法])がそのまま使えます。
  • アンサンブル・MC ドロップアウト:初期値や推論時ドロップアウトを変えた複数の予測のばらつきを、認識的不確実性の代用にする。理論的な保証は弱いものの、実装が容易なため広く使われています。

いずれの方法でも、目指しているものは Theorem 6.4 の「σ2\sigma^{2}ϕTSNϕ\boldsymbol{\phi}_{*}^{\mathsf{T}}S_N\boldsymbol{\phi}_{*} の分解」の一般化です。

Exercise 7.1

Example 4.4 の設定(nn 回中 kk 回表、事前分布 Beta(a,b)\mathrm{Beta}(a,b)、事後分布 Beta(k+a,nk+b)\mathrm{Beta}(k+a,\,n-k+b))で、事後平均が最尤推定量と事前平均の凸結合(重み付き平均)として

E[θD]=γkn+(1γ)aa+b,γ=nn+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}

と書けることを示してください。また nn\to\infty のときの挙動を述べてください。

Solution

Beta(α,β)\mathrm{Beta}(\alpha,\beta) の平均は α/(α+β)\alpha/(\alpha+\beta) なので、事後平均は

E[θD]=k+a(k+a)+(nk+b)=k+an+a+b\mathbb{E}[\theta\mid\mathcal{D}] = \frac{k+a}{(k+a)+(n-k+b)} = \frac{k+a}{n+a+b}

です(分母で kk が打ち消えます)。一方、右辺を計算すると

γkn+(1γ)aa+b=nn+a+bkn+a+bn+a+baa+b=kn+a+b+an+a+b=k+an+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}

となり一致します。ここで 1γ=1nn+a+b=a+bn+a+b1-\gamma = 1 - \frac{n}{n+a+b} = \frac{a+b}{n+a+b} を使いました。γ(0,1)\gamma\in(0,1) かつ γ+(1γ)=1\gamma+(1-\gamma)=1 なので、これは確かに凸結合です。

nn\to\infty のとき γ=nn+a+b1\gamma = \frac{n}{n+a+b}\to 1 なので、事後平均は最尤推定量 k/nk/n に近づきます。事前分布の影響は O(1/n)O(1/n) で消えていきます。逆に n=0n=0(データなし)なら γ=0\gamma=0 で事後平均は事前平均 a/(a+b)a/(a+b) そのものです。

aabb が「疑似観測回数」と呼ばれる理由もここにあります。事前分布 Beta(a,b)\mathrm{Beta}(a,b) は、実データに先立って表を aa 回、裏を bb 回見たのと同じ効き方をします。a=b=1a=b=1(一様事前分布)のとき事後平均は (k+1)/(n+2)(k+1)/(n+2) で、これはラプラスの継起の法則、機械学習の言葉では加算スムージングです。

Exercise 7.2標準

X1,,XnX_1,\ldots,X_n が独立に平均 λ>0\lambda > 0 のポアソン分布 p(xλ)=eλλx/x!p(x\mid\lambda) = e^{-\lambda}\lambda^{x}/x!x=0,1,2,x=0,1,2,\ldots)に従うとします。S=i=1nXiS=\sum_{i=1}^{n} X_i とおきます。

  1. λ\lambda の最尤推定量を求めてください(S=0S=0 の場合も論じること)。
  2. 事前分布としてガンマ分布 p(λ)λa1ebλp(\lambda)\propto \lambda^{a-1}e^{-b\lambda}a>0a>0b>0b>0λ>0\lambda>0)を取ったときの事後分布を求め、ガンマ分布がポアソン尤度の共役事前分布であることを示してください。
  3. MAP 推定量を求め、nn\to\infty での挙動を述べてください。
Solution

1. 尤度と対数尤度は

L(λ)=i=1neλλxixi!=enλλSi1xi!,(λ)=nλ+Slogλilog(xi!).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!).

S1S \ge 1 のとき (λ)=n+S/λ=0\ell'(\lambda) = -n + S/\lambda = 0 より λ=S/n\lambda = S/n、また (λ)=S/λ2<0\ell''(\lambda) = -S/\lambda^{2} < 0 なので \ell は狭義凹で、λ^ML=S/n=xˉ\hat{\lambda}_{\mathrm{ML}} = S/n = \bar{x} が唯一の最大点です。

S=0S=0 のときは (λ)=nλ\ell(\lambda) = -n\lambdaλ>0\lambda>0 で狭義単調減少なので、(0,)(0,\infty) 上に最大点は存在せず、λ0\lambda\downarrow 0 が上限を与えます。パラメータ空間を [0,)[0,\infty) に広げれば λ^ML=0\hat{\lambda}_{\mathrm{ML}}=0 で、これも xˉ=0\bar{x}=0 と書けます。Example 3.2 の端の場合と同じ病理です。

2. Theorem 2.2 の比例形より

p(λD)(enλλS)(λa1ebλ)=λ(S+a)1e(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}.

右辺は形状 S+aS+a、率 n+bn+b のガンマ分布 Gamma(S+a,n+b)\mathrm{Gamma}(S+a,\, n+b) の密度の λ\lambda 依存部分です。S+a>0S+a>0n+b>0n+b>0 なのでこれは正しく正規化可能な密度であり、事後分布は Gamma(S+a,n+b)\mathrm{Gamma}(S+a,\,n+b) です。事後分布が再びガンマ分布族に属するので、Definition 4.3 の意味で共役です。更新則は「形状に観測値の総和を足し、率に標本数を足す」となります。

3. 対数事後分布は (S+a1)logλ(n+b)λ(S+a-1)\log\lambda - (n+b)\lambda(定数を除く)で、微分すると S+a1λ(n+b)\frac{S+a-1}{\lambda} - (n+b) です。S+a>1S+a>1 のとき、これが 00 になる λ\lambda

λ^MAP=S+a1n+b\hat{\lambda}_{\mathrm{MAP}} = \frac{S+a-1}{n+b}

であり、二階微分 (S+a1)/λ2<0-(S+a-1)/\lambda^{2}<0 より唯一の最大点です。S+a1S+a\le 1 のときは対数事後分布が単調減少なので最頻値は λ=0\lambda=0 です。

nn\to\infty の挙動は、xˉ=S/n\bar{x}=S/n を使って

λ^MAP=nxˉ+a1n+b=xˉnn+b+a1n+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}

と書き直せば見えます。xˉ\bar{x} が大数の法則で真の λ0\lambda_0 に収束するなら、第 1 項は λ0\lambda_0 に、第 2 項は 00 に収束するので λ^MAPλ0\hat{\lambda}_{\mathrm{MAP}}\to\lambda_0 です。事前分布の影響は O(1/n)O(1/n) で消え、Exercise 7.1 と同じ結論になります。

Exercise 7.3標準

λ>0\lambda>0zRz\in\mathbb{R} を定数とし、1 変数関数

J(w)=12(wz)2+λw(wR)J(w) = \frac{1}{2}(w-z)^{2} + \lambda|w| \qquad (w\in\mathbb{R})

を考えます。

  1. JJ の最小点が唯一存在し、それがソフトしきい値関数 Sλ(z)=sign(z)max(zλ,0)S_\lambda(z) = \operatorname{sign}(z)\max(|z|-\lambda,\,0) で与えられることを示してください。
  2. L1 罰則を L2 罰則に替えた J2(w)=12(wz)2+λ2w2J_2(w) = \frac{1}{2}(w-z)^{2} + \frac{\lambda}{2}w^{2} の最小点を求め、1 と比較して「なぜ L1 だけが厳密な 00 を生むのか」を説明してください。
Solution

1. 12(wz)2\frac12(w-z)^2 は狭義凸、λw\lambda|w| は凸なので JJ は狭義凸です。したがって最小点は高々 1 つです。また λw0\lambda|w|\ge 0 より J(w)12(wz)2J(w)\ge\frac{1}{2}(w-z)^{2} であり、w|w|\to\infty のとき右辺が \infty に発散するので JJ は強制的です。連続な強制的関数は最小値を取るので、最小点はちょうど 1 つ存在します。凸関数では ww^{\star} が最小点であることと 0J(w)0\in\partial J(w^{\star}) が同値です。ここで劣微分は

J(w)={wz+λs  :  sw},w={{sign(w)}(w0)[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}

です。場合分けします。

  • w>0w>0:条件は wz+λ=0w-z+\lambda=0、すなわち w=zλw = z-\lambda。これが w>0w>0 を満たすのは z>λz>\lambda のときだけです。
  • w<0w<0:条件は wzλ=0w-z-\lambda=0、すなわち w=z+λw = z+\lambda。これが w<0w<0 を満たすのは z<λz<-\lambda のときだけです。
  • w=0w=0:条件は 0{z+λs:s[1,1]}=[zλ,z+λ]0\in\{-z+\lambda s: s\in[-1,1]\} = [-z-\lambda,\, -z+\lambda]、すなわち zλ0z+λ-z-\lambda\le 0\le -z+\lambda、整理して zλ|z|\le\lambda

3 つの場合は zz について排他的かつ網羅的で、まとめると

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}

です。z>λz>\lambda のとき sign(z)max(zλ,0)=zλ\operatorname{sign}(z)\max(|z|-\lambda,0) = z-\lambdaz<λz<-\lambda のとき (zλ)=(zλ)=z+λ-(|z|-\lambda) = -(-z-\lambda) = z+\lambdazλ|z|\le\lambda のとき max(zλ,0)=0\max(|z|-\lambda,0)=0 なので、これは Sλ(z)S_\lambda(z) に一致します。

2. J2J_2 は滑らかで狭義凸なので J2(w)=(wz)+λw=0J_2'(w) = (w-z)+\lambda w = 0、すなわち w=z/(1+λ)w^{\star} = z/(1+\lambda) が唯一の最小点です。これが 00 になるのは z=0z=0 のときだけです。

違いは原点での「傾きの跳び」にあります。λw\lambda|w| の劣微分は w=0w=0 で幅 2λ2\lambda の区間 [λ,λ][-\lambda,\lambda] を持つので、データ由来の勾配 z-z がこの区間に収まっている限り w=0w=0 が最適であり続けます。zλ|z|\le\lambda という「不感帯」が生じるのはこのためです。一方 λ2w2\frac{\lambda}{2}w^{2} の微分 λw\lambda ww=0w=000 になり、罰則が 00 の近傍で ww を押し戻す力を失うため、どんなに小さくても z0z\ne0 なら w0w^\star\ne0 になります。

事前分布の言葉では、Proposition 5.3 のラプラス密度が原点に尖った角を持つのに対し、Theorem 5.1 のガウス密度は原点で滑らかである、という違いに対応します。ただし Remark 5.4 で述べたとおり、これはあくまで最頻値の性質です。

Exercise 7.4

ロジスティック回帰を扱います。σ(t)=1/(1+et)\sigma(t)=1/(1+e^{-t})、データ (xi,yi)i=1n(\boldsymbol{x}_i,y_i)_{i=1}^{n}xiRM\boldsymbol{x}_i\in\mathbb{R}^{M}yi{0,1}y_i\in\{0,1\})に対しモデルを p(y=1x,w)=σ(wTx)p(y=1\mid\boldsymbol{x},\boldsymbol{w}) = \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}) とします。

  1. 負の対数尤度が交差エントロピー損失 E(w)=i[yilogp^i+(1yi)log(1p^i)]E(\boldsymbol{w}) = -\sum_i\bigl[y_i\log\hat{p}_i + (1-y_i)\log(1-\hat{p}_i)\bigr]p^i=σ(wTxi)\hat p_i = \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i))に一致することを示し、さらに ti=(2yi1)wTxit_i = (2y_i-1)\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i とおくと E(w)=ilog(1+eti)E(\boldsymbol{w}) = \sum_i \log(1+e^{-t_i}) と書けることを示してください。
  2. 事前分布 wN(0,τ2I)\boldsymbol{w}\sim\mathcal{N}(\boldsymbol{0},\tau^{2}I) の下で MAP 推定が E(w)+12τ2w2E(\boldsymbol{w}) + \frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2} の最小化になることを示してください。
  3. データが線形分離可能、すなわちある w0\boldsymbol{w}_0w0=1\lVert\boldsymbol{w}_0\rVert=1)が存在してすべての ii(2yi1)w0Txi>0(2y_i-1)\boldsymbol{w}_0^{\mathsf{T}}\boldsymbol{x}_i>0 が成り立つとします。このとき最尤推定量は存在しないが、2 の MAP 推定量はただ 1 つ存在することを示してください。
Solution

1. yi{0,1}y_i\in\{0,1\} なので p(yixi,w)=p^iyi(1p^i)1yip(y_i\mid\boldsymbol{x}_i,\boldsymbol{w}) = \hat p_i^{\,y_i}(1-\hat p_i)^{1-y_i} です(yi=1y_i=1 なら p^i\hat p_iyi=0y_i=0 なら 1p^i1-\hat p_i になります)。各データが独立なので対数尤度は和になり、

logp(Dw)=i[yilogp^i+(1yi)log(1p^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}).

次に 1σ(t)=111+et=et1+et=11+et=σ(t)1-\sigma(t) = 1 - \frac{1}{1+e^{-t}} = \frac{e^{-t}}{1+e^{-t}} = \frac{1}{1+e^{t}} = \sigma(-t) に注意します。yi=1y_i=1 なら p^i=σ(wTxi)=σ(ti)\hat p_i = \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i) = \sigma(t_i)yi=0y_i=0 なら 1p^i=σ(wTxi)=σ(ti)1-\hat p_i = \sigma(-\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i) = \sigma(t_i)2yi1=12y_i-1=-1 だから)なので、どちらの場合も ii 番目の項は logσ(ti)-\log\sigma(t_i) です。logσ(t)=log(1+et)-\log\sigma(t) = \log(1+e^{-t}) なので E(w)=ilog(1+eti)E(\boldsymbol{w}) = \sum_i\log(1+e^{-t_i}) を得ます。

2. Proposition 4.2 より MAP は logp(Dw)+logp(w)\log p(\mathcal{D}\mid\boldsymbol{w}) + \log p(\boldsymbol{w}) の最大点です。1 より第 1 項は E(w)-E(\boldsymbol{w})、第 2 項は Theorem 5.1 の証明第 1 段と同じ計算で const12τ2w2\mathrm{const} - \frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2} です。符号を反転して定数を落とせば、F(w)=E(w)+12τ2w2F(\boldsymbol{w}) = E(\boldsymbol{w}) + \frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2} の最小化になります。これは本文と同じ E(w)+λw2E(\boldsymbol{w}) + \lambda\lVert\boldsymbol{w}\rVert^{2} という書き方をすれば、正則化係数 λ=1/(2τ2)\lambda = 1/(2\tau^{2}) の L2 正則化つきロジスティック回帰にほかなりません。

3. 最尤推定量の非存在。 γ=mini(2yi1)w0Txi\gamma = \min_i (2y_i-1)\boldsymbol{w}_0^{\mathsf{T}}\boldsymbol{x}_i とおくと、有限個の正数の最小値なので γ>0\gamma>0 です。w=cw0\boldsymbol{w}=c\boldsymbol{w}_0c>0c>0)と取ると ti=c(2yi1)w0Txicγt_i = c(2y_i-1)\boldsymbol{w}_0^{\mathsf{T}}\boldsymbol{x}_i \ge c\gamma なので、log(1+et)\log(1+e^{-t})tt について減少することから

E(cw0)=i=1nlog(1+eti)nlog(1+ecγ)c0.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.

一方、任意の有限な w\boldsymbol{w} に対し eti>0e^{-t_i}>0 なので log(1+eti)>0\log(1+e^{-t_i})>0、したがって E(w)>0E(\boldsymbol{w})>0 です。よって infwE=0\inf_{\boldsymbol{w}} E = 0 ですが、この下限を達成する w\boldsymbol{w} は存在しません。EE の最小化は尤度の最大化と同値なので、最尤推定量は存在しません(数値的には重みが発散します)。これは Theorem 6.3[Logistic Regression] を、いま用意した事前分布の言葉で言い直したものです。

MAP 推定量の存在。 E(w)0E(\boldsymbol{w})\ge 0 より

F(w)=E(w)+12τ2w212τ2w2F(\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}

なので、w\lVert\boldsymbol{w}\rVert\to\infty のとき F(w)F(\boldsymbol{w})\to\infty(強制的)です。いま F(0)=E(0)=nlog2F(\boldsymbol{0}) = E(\boldsymbol{0}) = n\log 2 なので、劣位集合 C={w:F(w)nlog2}C = \{\boldsymbol{w} : F(\boldsymbol{w})\le n\log 2\} は空でなく、その上では 12τ2w2nlog2\frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2}\le n\log 2 すなわち wτ2nlog2\lVert\boldsymbol{w}\rVert\le\tau\sqrt{2n\log 2} となるので有界です。FF は連続なので CC は閉、よって CC はコンパクトです。ワイエルシュトラスの定理より FFCC 上で最小値を取り、CC の外では F>nlog2minCFF>n\log 2\ge \min_C F なので、それは RM\mathbb{R}^{M} 全体での最小値です。

一意性。 g(t)=log(1+et)g(t)=\log(1+e^{-t})g(t)=et(1+et)2>0g''(t) = \frac{e^{-t}}{(1+e^{-t})^{2}}>0 より凸で、tit_iw\boldsymbol{w} の線形関数なので、凸関数と線形写像の合成である g(ti(w))g(t_i(\boldsymbol{w})) は凸、その和 EE も凸です。12τ2w2\frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2}τ2>0\tau^2 > 0 より狭義凸なので、FF は狭義凸です。狭義凸関数の最小点は高々 1 つなので、MAP 推定量はただ 1 つに定まります。

解釈。 分離可能なデータでは「境界からできるだけ遠ざけたい」という力が無限に働き、最尤推定は破綻します。ガウス事前分布はこの力に 12τ2w2\frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{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

この補題は何のためにあるか。 Theorem 6.1 の証明では、事後密度の対数が「w\boldsymbol{w} の 2 次形式 ++ 1 次形式 ++ 定数」の形になることまでを示しました。そこから「だからガウス分布である」と言うには、その形の関数が実際にどのガウス分布の密度なのかを特定する必要があります。行列版の平方完成がその作業です。

Lemma 7.52 次形式で書かれた密度の同定

AAMM 次の対称正定値行列、bRM\boldsymbol{b}\in\mathbb{R}^{M}C>0C>0 を定数とする。RM\mathbb{R}^{M} 上の確率密度 qq

q(w)=Cexp ⁣(12wTAw+bTw)(wRM)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})

を満たすならば、qq は多変量正規分布 N(A1b,A1)\mathcal{N}(A^{-1}\boldsymbol{b},\, A^{-1}) の密度である。

Proof(Lemma 7.5)

AA は正定値なので可逆で、A1A^{-1} も対称正定値です(A=ATA = A^{\mathsf{T}} より (A1)T=(AT)1=A1(A^{-1})^{\mathsf{T}} = (A^{\mathsf{T}})^{-1} = A^{-1}、また AA の固有値がすべて正なら A1A^{-1} の固有値もすべて正)。μ=A1b\boldsymbol{\mu} = A^{-1}\boldsymbol{b} とおきます。

まず恒等式を確かめます。展開すると

12(wμ)TA(wμ)=12(wTAwwTAμμTAw+μTAμ)-\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)

ですが、Aμ=AA1b=bA\boldsymbol{\mu} = AA^{-1}\boldsymbol{b} = \boldsymbol{b} なので wTAμ=wTb=bTw\boldsymbol{w}^{\mathsf{T}}A\boldsymbol{\mu} = \boldsymbol{w}^{\mathsf{T}}\boldsymbol{b} = \boldsymbol{b}^{\mathsf{T}}\boldsymbol{w}(スカラーなので転置しても同じ)、同様に μTAw=(Aμ)Tw=bTw\boldsymbol{\mu}^{\mathsf{T}}A\boldsymbol{w} = (A\boldsymbol{\mu})^{\mathsf{T}}\boldsymbol{w} = \boldsymbol{b}^{\mathsf{T}}\boldsymbol{w}AA の対称性を使いました)、そして μTAμ=bTA1b\boldsymbol{\mu}^{\mathsf{T}}A\boldsymbol{\mu} = \boldsymbol{b}^{\mathsf{T}}A^{-1}\boldsymbol{b} です。したがって

12(wμ)TA(wμ)=12wTAw+bTw12bTA1b.-\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}.

これを qq の式に代入すると、K=Cexp(12bTA1b)>0K = C\exp\bigl(\tfrac12\boldsymbol{b}^{\mathsf{T}}A^{-1}\boldsymbol{b}\bigr) > 0 とおいて

q(w)=Kexp ⁣(12(wμ)TA(wμ)).q(\boldsymbol{w}) = K \exp\!\left( -\frac{1}{2}(\boldsymbol{w}-\boldsymbol{\mu})^{\mathsf{T}}A(\boldsymbol{w}-\boldsymbol{\mu}) \right).

一方、Σ=A1\Sigma = A^{-1} とした多変量正規分布 N(μ,Σ)\mathcal{N}(\boldsymbol{\mu},\Sigma) の密度は

g(w)=1(2π)M/2(detΣ)1/2exp ⁣(12(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)

であり、Σ1=A\Sigma^{-1}=A なので指数部は qq のそれと完全に一致します。つまり q=(K/K)gq = (K/K')\,gKK'gg の正規化定数)という比例関係が全点で成り立ちます。

最後に比例定数が 11 であることを言います。qqgg も確率密度なので q=g=1\int q = \int g = 1 です。q=(K/K)gq = (K/K')g の両辺を積分すると 1=(K/K)11 = (K/K')\cdot 1、すなわち K=KK = K' を得ます。よって q=gq = g、つまり qqN(μ,A1)=N(A1b,A1)\mathcal{N}(\boldsymbol{\mu},A^{-1}) = \mathcal{N}(A^{-1}\boldsymbol{b},\,A^{-1}) の密度です。

1 次元で確かめておきます。 M=1M=1A=a>0A=a>0bRb\in\mathbb{R} なら q(w)=Cexp(12aw2+bw)q(w) = C\exp(-\tfrac12 a w^{2}+bw) で、補題は「これは平均 b/ab/a、分散 1/a1/a の正規分布」と言っています。実際 12aw2+bw=a2(wba)2+b22a-\tfrac12 a w^2 + bw = -\tfrac{a}{2}\bigl(w - \tfrac{b}{a}\bigr)^{2} + \tfrac{b^{2}}{2a} なので、確かに平均 b/ab/a、分散 1/a1/a です。aa は分散の逆数、すなわち精度です。Theorem 6.1SN1=S01+σ2ΦTΦS_N^{-1} = S_0^{-1}+\sigma^{-2}\Phi^{\mathsf{T}}\Phi が「精度は足し算になる」と読めるのは、このためです。

Report an error in this article ・Operated by: Mugen Giken LLCPricingTermsLegal notice

© 2026 夢現技研合同会社 ・Feeding the text to an LLM is welcome. Code samples are MIT licensed.