コンテンツにスキップ

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

対象読者: 学部レベルの確率(条件付き確率・ガウス分布)と線形代数・微分積分を学び、本テーマ第 1〜6 章の内容(最尤推定・線形回帰・ロジスティック回帰)を読んだ人。ベイズ統計は仮定しません。

前提となる記事:主成分分析:分散最大化はなぜ固有値問題になるのか

生 Markdown

コインを 3 回投げて 3 回とも表が出たとき、「このコインは絶対に裏が出ない」と結論してよいでしょうか。直感では「まだ決めつけられない」と感じるはずです。この記事は、その「決めつけない」感覚を、確率の言葉で正確に書く方法を扱います。同じ考え方で、機械学習でおなじみの「重みを小さく抑える工夫」(正則化)の正体も分かり、さらに予測が「どのくらい当てになるか」まで数字で出せるようになります。

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

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

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

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

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

  1. なぜ残差の 2 乗なのですか。絶対値でも 4 乗でもよいはずです。
  2. なぜ罰則が ∥w∥2\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) が定まる仕組みのことです。

定義 2.1(尤度・事前分布・事後分布)↗

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

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

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

L(θ)=∏i=1Np(zi∣θ),ℓ(θ)=∑i=1Nlog⁡p(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 を信じたとき、手元のデータはどれくらい出やすかったか」を測る物差しです。

定理 2.2(ベイズの定理)✓

θ\theta と D\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 事前分布」である。

証明(定理 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 の等号は、定義 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 を含まない因子をすべて捨てて、残った形が何という分布かを見抜くことに尽きます。

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

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

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

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

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

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

注意 2.4(頻度論とベイズの分かれ目)

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

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

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

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

定義 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⁡θ∈Θ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)

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

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

例 3.2(コイン投げの最尤推定)

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

L(θ)=θk(1−θ)n−k,ℓ(θ)=klog⁡θ+(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).

微分して

ℓ′(θ)=kθ−n−k1−θ=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)}.

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

ℓ′′(θ)=−kθ2−n−k(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 乗和」が、実はガウスノイズの仮定と同じものであることを示します。

命題 3.3(ガウス雑音の下で最尤推定は最小二乗法)✓

入力 x1,…,xN∈Rd\boldsymbol{x}_1,\ldots,\boldsymbol{x}_N \in \mathbb{R}^{d} は固定された既知の値とし、特徴写像 ϕ:Rd→RM\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 は独立で ε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=(y1,…,yN)T\boldsymbol{y} = (y_1,\ldots,y_N)^{\mathsf{T}} とおくと、w∈RM\boldsymbol{w} \in \mathbb{R}^{M} の最尤推定量は

w^ML=arg min⁡w∈RM∥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}

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

証明(命題 3.3)

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

p(yi∣w)=12πσ2exp⁡ ⁣(−(yi−wTϕ(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} を与えたとき)独立で、定義 2.1 より尤度は積になります。対数を取ると

ℓ(w)=∑i=1Nlog⁡p(yi∣w)=−N2log⁡(2πσ2)  −  12σ2∑i=1N(yi−wTϕ(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(yi−wTϕ(xi))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} は、計画行列の定義から Φw\Phi\boldsymbol{w} の第 ii 成分が ϕ(xi)Tw\boldsymbol{\phi}(\boldsymbol{x}_i)^{\mathsf{T}}\boldsymbol{w} であることによります。したがって

ℓ(w)=const−12σ2∥y−Φw∥2.\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−Φw∥2\lVert \boldsymbol{y} - \Phi\boldsymbol{w}\rVert^{2} を最小化することは同値です。

後半は正規方程式です。J(w)=∥y−Φw∥2J(\boldsymbol{w}) = \lVert \boldsymbol{y}-\Phi\boldsymbol{w}\rVert^{2} を展開すると J=yTy−2yTΦ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 は正定値(v≠0\boldsymbol{v} \ne \boldsymbol{0} に対し vTΦTΦv=∥Φv∥2>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. 最大事後確率推定:事前分布を足す”

定義 4.1(MAP 推定量)↗

事前分布 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 推定量)といいます。

命題 4.2(MAP の目的関数)✓

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

θ^MAP=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\}

が成り立つ。とくに事前分布が Θ\Theta 上で定数(一様)ならば θ^MAP=θ^ML\hat{\theta}_{\mathrm{MAP}} = \hat{\theta}_{\mathrm{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(\mathcal{D}) は θ\theta に依存しない正の定数なので、p(θ∣D)p(\theta\mid\mathcal{D}) を最大にする θ\theta と p(D∣θ)p(θ)p(\mathcal{D}\mid\theta)p(\theta) を最大にする θ\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) を最大にする θ\theta とも同じです。

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

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

定義 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} に属することをいいます。

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

次の例では、コインの表の確率 θ\theta の事前分布としてベータ分布を使います。ベータ分布は区間 (0,1)(0,1) 上の分布で、密度が θa−1(1−θ)b−1\theta^{a-1}(1-\theta)^{b-1} に比例するものです(a,ba,b で山の位置や幅が変わります)。密度は全体の積分が 11 になるよう定数を決めれば形から一意に定まるので、「θ\theta に依存する部分の形が同じなら同じ分布」と読み取れます。この例のねらいは、そう読み取ることで積分を一度も計算せずに事後分布が求まる、と確かめることです。

例 4.4(ベータ事前分布とコイン投げ)

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

p(θ)=1B(a,b) θa−1(1−θ)b−1,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 に依存しない正規化定数です。定理 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}.

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

具体的な数値を最後まで追います。n=3n = 3、k=3k = 3(3 回投げて 3 回とも表)、事前分布は Beta(2,2)\mathrm{Beta}(2,2)(「0.50.5 のあたりが怪しい」という穏やかな信念、表裏 2 回ずつを見たのと同じ重み)とします。事後分布は 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−θ=4−5θθ(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/θ2−1/(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/7≈0.7145/7 \approx 0.714 です。

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

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

という 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 = 1 で 00 に落ちるので、積は 0.80.8 付近で山を作ります。データが増えれば尤度の山が鋭くなり、事前分布の影響は相対的に小さくなっていきます。

次の注意のねらいは、MAP が「パラメータの名前の付け方」で変わってしまうことを、同じ事後分布を使った具体的な計算で確かめることです。θ\theta を θ2\theta^{2} に取り替えるだけで、山の頂上の位置が動きます。

注意 4.5(MAP は座標変換で不変ではない)

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

例 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/2−2η)=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. 正則化の正体:罰則項は事前分布である”

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

次の定理の考え方は一言でいえば、「対数をとると、ガウス分布の指数部分がそのまま 2 乗の罰則になる」ということです。ガウス分布の密度は exp⁡(−(2 乗))\exp(-(\text{2 乗})) の形なので、その対数は 2 乗の項になり、MAP の目的関数「対数尤度 ++ 対数事前分布」に足すと、損失に罰則を加えたものに見えます。

定理 5.1(L2 正則化はガウス事前分布の MAP 推定)✓

命題 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))に加えて、重みの事前分布を

w∼N(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 min⁡w∈RM{∥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}}

で与えられ、これはリッジ回帰(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}

である。

証明(定理 5.1)

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

log⁡p(y∣w)=−N2log⁡(2πσ2)−12σ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}

です(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⁡(−∥w∥2/(2τ2))p(\boldsymbol{w}) = (2\pi\tau^{2})^{-M/2}\exp\bigl(-\lVert\boldsymbol{w}\rVert^{2}/(2\tau^{2})\bigr) なので

log⁡p(w)=−M2log⁡(2πτ2)−12τ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}.

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

log⁡p(y∣w)+log⁡p(w)=const−12σ2∥y−Φw∥2−12τ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}

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

w^MAP=arg min⁡w{12σ2∥y−Φw∥2+12τ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\}

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

w^MAP=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\}

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

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

J(w)=yTy−2yTΦ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 は正定値です。実際、v≠0\boldsymbol{v}\ne\boldsymbol{0} に対し

vTAv=∥Φ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

となります(λ=σ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=A−1Φ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 で最尤推定に戻り、これは 命題 4.2 の後半(一様事前分布なら MLE)とも整合します。

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

ΦTΦ\Phi^{\mathsf{T}}\Phi は対称半正定値なので、スペクトル定理(系 4.3[スペクトル定理])により直交行列 UU と対角行列 Λ=diag⁡(d1,…,dM)\Lambda = \operatorname{diag}(d_1,\ldots,d_M)(dj≥0d_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}} なので、定理 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 = 1、d1=100d_1 = 100、d2=0.01d_2 = 0.01 なら、縮小率は第 1 方向で 100/101≈0.990100/101 \approx 0.990、第 2 方向で 0.01/1.01≈0.00990.01/1.01 \approx 0.0099 です。データの分散が大きい方向はほとんど手つかず、分散がほとんどない方向は 1/1001/100 に潰されます。

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

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

命題 5.3(L1 正則化はラプラス事前分布の MAP 推定)✓

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

p(w)=∏j=1M12bexp⁡ ⁣(−∣wj∣b)=1(2b)Mexp⁡ ⁣(−∥w∥1b)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)

に従うものとする(∥w∥1=∑j∣wj∣\lVert\boldsymbol{w}\rVert_1 = \sum_j |w_j|)。このとき

w^MAP=arg min⁡w∈RM{∥y−Φw∥2+λ∥w∥1},λ=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 正則化つき最小二乗法)の解に一致する。

証明(命題 5.3)

対数事前分布は

log⁡p(w)=−Mlog⁡(2b)−1b∥w∥1\log p(\boldsymbol{w}) = -M\log(2b) - \frac{1}{b}\lVert\boldsymbol{w}\rVert_{1}

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

const−12σ2∥y−Φw∥2−1b∥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}

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

w^MAP=arg min⁡w{∥y−Φw∥2+2σ2b∥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\}

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

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

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

罰則項対応する事前分布事前密度(θ\theta 依存部分)正則化係数解の特徴
なし一様(非正則な事前分布)∝1\propto 1λ=0\lambda = 0最尤推定に一致
λ∥w∥22\lambda \lVert \boldsymbol{w}\rVert_{2}^{2}ガウス N(0,τ2I)\mathcal{N}(\boldsymbol{0},\tau^{2}I)exp⁡(−∥w∥22/(2τ2))\exp\bigl(-\lVert\boldsymbol{w}\rVert_2^{2}/(2\tau^{2})\bigr)λ=σ2/τ2\lambda = \sigma^{2}/\tau^{2}全成分を滑らかに縮小
λ∥w∥1\lambda \lVert \boldsymbol{w}\rVert_{1}ラプラス(尺度 bb)exp⁡(−∥w∥1/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 行目は、事前学習済みモデルからの微調整で「元の重みから離れすぎない」ようにする正則化が、そのまま「元の重みを平均とするガウス事前分布」であることを示しています。正則化を発明するとは、事前分布を設計することです。

-303wラプラス(尖った山)ガウス(滑らかな山)
分散を 1 にそろえたガウス分布(実線)とラプラス分布(尺度 b = 1/√2、破線)の密度。ラプラス分布は原点で尖り、裾は重い

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

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

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

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

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

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

次の定理の証明の方針は 3 段です。(1) ベイズの定理の比例形で「尤度 ×\times 事前分布」を書く。(2) 指数の中を w\boldsymbol{w} の 2 次式に整理する(長い行列の展開ですが、w\boldsymbol{w} について 2 次の項と 1 次の項を集めるだけの作業です)。(3) 付録の補題で、その形の密度がどのガウス分布かを読み取る。本質は (3) で、2 次式の形をしていればガウス分布だ、という点です。

定理 6.1(ベイズ線形回帰の事後分布)✓

命題 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_0 は MM 次の対称正定値行列、m0∈RM\boldsymbol{m}_0\in\mathbb{R}^M)とし、w\boldsymbol{w} と雑音は独立とする。このとき事後分布は再びガウス分布

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

であり、

SN=(S0−1+1σ2ΦTΦ)−1,mN=SN(S0−1m0+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)

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

証明(定理 6.1)

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

vTAv=vTS0−1v+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

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

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

p(w∣y)∝p(y∣w) p(w)∝exp⁡(−12σ2∥y−Φw∥2−12(w−m0)TS0−1(w−m0))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−Φw∥2=yTy−2yTΦ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},(w−m0)TS0−1(w−m0)=wTS0−1w−2m0TS0−1w+m0TS0−1m0(\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

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

−12wT(S0−1+1σ2ΦTΦ)⏟=Aw  +  (S0−1m0+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(w∣y)=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 段より対称正定値なので、補題 7.5 が適用でき、この密度は N(A−1b, A−1)\mathcal{N}(A^{-1}\boldsymbol{b},\, A^{-1}) の密度に一致します。A−1=SNA^{-1}=S_N、A−1b=SN(S0−1m0+σ−2ΦTy)=mNA^{-1}\boldsymbol{b} = S_N(S_0^{-1}\boldsymbol{m}_0 + \sigma^{-2}\Phi^{\mathsf{T}}\boldsymbol{y}) = \boldsymbol{m}_N なので主張が従います。

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

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

w^MAP=E[w∣y]=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}}

が成り立つ。

証明(系 6.2)

ガウス分布の密度 ∝exp⁡(−12(w−mN)TSN−1(w−mN))\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 を取り、SN−1S_N^{-1} が正定値なので他の点では真に負です。よって最頻値は平均 mN\boldsymbol{m}_N に一致します。

具体形は代入するだけです。S0−1=τ−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}.

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

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

定義 6.3(事後予測分布)↗

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

p(y∗∣x∗,D)=∫p(y∗∣x∗,w) p(w∣D) 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(y∗∣x∗,w^)p(y_{*}\mid\boldsymbol{x}_{*},\hat{\boldsymbol{w}}) を使うプラグイン予測と対比されます。

次の定理は、「新しい入力での予測のばらつき」を、観測ノイズの分だけの部分とパラメータが分からない分の部分に分けて数式で書くものです。証明では、ガウス分布に従う量の線形結合や和がまたガウス分布になることを使います。特性関数(確率分布を E[eitX]\mathbb{E}[e^{itX}] という関数で表したもの)は、その確認に使う道具です。

定理 6.4(ベイズ線形回帰の事後予測分布)✓

定理 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(y∗∣x∗,D)=N(y∗ ∣ mNTϕ∗,  σ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{パラメータの不確実性(認識的不確実性)}}

と分解される。

証明(定理 6.4)

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

E[eitaTw]=exp⁡ ⁣(it aTmN−t22aTSNa)\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)

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

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

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

M=1M=1、ϕ(x)=x\boldsymbol{\phi}(x)=x(原点を通る直線 y=wxy = wx)とし、σ2=1\sigma^{2}=1、事前分布は w∼N(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=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.

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

SN−1=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)=16⋅8=43≈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.

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

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

x∗x_{*}予測平均 43x∗\tfrac{4}{3}x_{*}認識的分散 x∗2/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 を使ったプラグイン予測をしていたら、どの x∗x_{*} でも分散は σ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 の中にしか書かれていません。

0510x予測平均プラグイン予測の 95% 区間事後予測分布の 95% 区間データ
上の数値例(データ 2 点)の予測平均 4x/3(実線)と 95% 予測区間(帯)。破線は点推定を使ったプラグイン予測の区間で、幅は一定です。データのない遠方ほど帯が扇状に広がります。
尤度・事前分布から、最尤推定・MAP 推定・事後予測分布へ至る道筋

注意 6.6(現実のモデルではどうするか)

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

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

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

演習 7.1易

例 4.4 の設定(nn 回中 kk 回表、事前分布 Beta(a,b)\mathrm{Beta}(a,b)、事後分布 Beta(k+a, n−k+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}

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

解答

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

E[θ∣D]=k+a(k+a)+(n−k+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+b⋅kn+a+bn+a+b⋅aa+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−γ=1−nn+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 なので、これは確かに凸結合です。

n→∞n\to\infty のとき γ=nn+a+b→1\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) そのものです。

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

演習 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(λ)∝λa−1e−bλp(\lambda)\propto \lambda^{a-1}e^{-b\lambda}(a>0a>0、b>0b>0、λ>0\lambda>0)を取ったときの事後分布を求め、ガンマ分布がポアソン尤度の共役事前分布であることを示してください。
  3. MAP 推定量を求め、n→∞n\to\infty での挙動を述べてください。
解答

1. 尤度と対数尤度は

L(λ)=∏i=1ne−λλxixi!=e−nλλS∏i1xi!,ℓ(λ)=−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!).

S≥1S \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 と書けます。例 3.2 の端の場合と同じ病理です。

2. 定理 2.2 の比例形より

p(λ∣D)∝(e−nλλS)⋅(λa−1e−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}.

右辺は形状 S+aS+a、率 n+bn+b のガンマ分布 Gamma(S+a, n+b)\mathrm{Gamma}(S+a,\, n+b) の密度の λ\lambda 依存部分です。S+a>0S+a>0、n+b>0n+b>0 なのでこれは正しく正規化可能な密度であり、事後分布は Gamma(S+a, n+b)\mathrm{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λ−(n+b)\frac{S+a-1}{\lambda} - (n+b) です。S+a>1S+a>1 のとき、これが 00 になる λ\lambda は

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

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

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

λ^MAP=nxˉ+a−1n+b=xˉ⋅nn+b+a−1n+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) で消え、演習 7.1 と同じ結論になります。

演習 7.3標準

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

J(w)=12(w−z)2+λ∣w∣(w∈R)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(w−z)2+λ2w2J_2(w) = \frac{1}{2}(w-z)^{2} + \frac{\lambda}{2}w^{2} の最小点を求め、1 と比較して「なぜ L1 だけが厳密な 00 を生むのか」を説明してください。
解答

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

∂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}

です。場合分けします。

  • w>0w>0:条件は w−z+λ=0w-z+\lambda=0、すなわち w=z−λw = z-\lambda。これが w>0w>0 を満たすのは z>λz>\lambda のときだけです。
  • w<0w<0:条件は w−z−λ=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−λ≤0≤−z+λ-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-\lambda、z<−λz<-\lambda のとき −(∣z∣−λ)=−(−z−λ)=z+λ-(|z|-\lambda) = -(-z-\lambda) = z+\lambda、∣z∣≤λ|z|\le\lambda のとき max⁡(∣z∣−λ,0)=0\max(|z|-\lambda,0)=0 なので、これは Sλ(z)S_\lambda(z) に一致します。

2. J2J_2 は滑らかで狭義凸なので J2′(w)=(w−z)+λ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 w は w=0w=0 で 00 になり、罰則が 00 の近傍で ww を押し戻す力を失うため、どんなに小さくても z≠0z\ne0 なら w⋆≠0w^\star\ne0 になります。

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

演習 7.4難

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

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

1. yi∈{0,1}y_i\in\{0,1\} なので p(yi∣xi,w)=p^i yi(1−p^i)1−yip(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_i、yi=0y_i=0 なら 1−p^i1-\hat p_i になります)。各データが独立なので対数尤度は和になり、

log⁡p(D∣w)=∑i[yilog⁡p^i+(1−yi)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}).

次に 1−σ(t)=1−11+e−t=e−t1+e−t=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 なら 1−p^i=σ(−wTxi)=σ(ti)1-\hat p_i = \sigma(-\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i) = \sigma(t_i)(2yi−1=−12y_i-1=-1 だから)なので、どちらの場合も ii 番目の項は −log⁡σ(ti)-\log\sigma(t_i) です。−log⁡σ(t)=log⁡(1+e−t)-\log\sigma(t) = \log(1+e^{-t}) なので E(w)=∑ilog⁡(1+e−ti)E(\boldsymbol{w}) = \sum_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}) の最大点です。1 より第 1 項は −E(w)-E(\boldsymbol{w})、第 2 項は 定理 5.1 の証明第 1 段と同じ計算で const−12τ2∥w∥2\mathrm{const} - \frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2} です。符号を反転して定数を落とせば、F(w)=E(w)+12τ2∥w∥2F(\boldsymbol{w}) = E(\boldsymbol{w}) + \frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2} の最小化になります。これは本文と同じ E(w)+λ∥w∥2E(\boldsymbol{w}) + \lambda\lVert\boldsymbol{w}\rVert^{2} という書き方をすれば、正則化係数 λ=1/(2τ2)\lambda = 1/(2\tau^{2}) の L2 正則化つきロジスティック回帰にほかなりません。

3. 最尤推定量の非存在。 γ=min⁡i(2yi−1)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}_0(c>0c>0)と取ると ti=c(2yi−1)w0Txi≥cγt_i = c(2y_i-1)\boldsymbol{w}_0^{\mathsf{T}}\boldsymbol{x}_i \ge c\gamma なので、log⁡(1+e−t)\log(1+e^{-t}) が tt について減少することから

E(cw0)=∑i=1nlog⁡(1+e−ti)≤nlog⁡(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.

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

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

F(w)=E(w)+12τ2∥w∥2≥12τ2∥w∥2F(\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)=nlog⁡2F(\boldsymbol{0}) = E(\boldsymbol{0}) = n\log 2 なので、劣位集合 C={w:F(w)≤nlog⁡2}C = \{\boldsymbol{w} : F(\boldsymbol{w})\le n\log 2\} は空でなく、その上では 12τ2∥w∥2≤nlog⁡2\frac{1}{2\tau^{2}}\lVert\boldsymbol{w}\rVert^{2}\le n\log 2 すなわち ∥w∥≤τ2nlog⁡2\lVert\boldsymbol{w}\rVert\le\tau\sqrt{2n\log 2} となるので有界です。FF は連続なので CC は閉、よって CC はコンパクトです。ワイエルシュトラスの定理より FF は CC 上で最小値を取り、CC の外では F>nlog⁡2≥min⁡CFF>n\log 2\ge \min_C F なので、それは RM\mathbb{R}^{M} 全体での最小値です。

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

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

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

補題 7.5(2 次形式で書かれた密度の同定)✓

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

q(w)=Cexp⁡ ⁣(−12wTAw+bTw)(w∈RM)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(A−1b, A−1)\mathcal{N}(A^{-1}\boldsymbol{b},\, A^{-1}) の密度である。

証明(補題 7.5)

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

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

−12(w−μ)TA(w−μ)=−12(wTAw−wTAμ−μ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μ=AA−1b=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μ=bTA−1b\boldsymbol{\mu}^{\mathsf{T}}A\boldsymbol{\mu} = \boldsymbol{b}^{\mathsf{T}}A^{-1}\boldsymbol{b} です。したがって

−12(w−μ)TA(w−μ)=−12wTAw+bTw−12bTA−1b.-\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⁡(12bTA−1b)>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).

一方、Σ=A−1\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')\,g(K′K' は gg の正規化定数)という比例関係が全点で成り立ちます。

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

1 次元で確かめておきます。 M=1M=1、A=a>0A=a>0、b∈Rb\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(w−ba)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 は分散の逆数、すなわち精度です。定理 6.1 の SN−1=S0−1+σ−2ΦTΦS_N^{-1} = S_0^{-1}+\sigma^{-2}\Phi^{\mathsf{T}}\Phi が「精度は足し算になる」と読めるのは、このためです。

このテーマの 100 本ノック(50 問)

執筆: Claude Code(Sonnet 5.5)・監査 3 回(Claude Code)・最終 2026-10-12

この記事の誤りを報告する ・運営: 夢現技研合同会社 ・料金 ・利用条件 ・特定商取引法に基づく表記

© 2026 夢現技研合同会社 ・コード例は MIT ライセンスです。