ロジスティック回帰のラプラス近似によるベイズ推定
ロジスティック回帰
特徴量（説明変数）ベクトル$\boldsymbol{x}_n$に対するバイナリ値を取る目的変数値$y_n \in \{0,1\}$の組が$N$個得られているとします。
ロジスティック回帰は$(\boldsymbol{x}_n,y_n)$を
\begin{align}
y_n &\sim {\rm Bern}(y_n; \sigma(\boldsymbol{x}_n^T \boldsymbol{w})) \\
\sigma(x) &= \frac{1}{1+\exp(-x)}
\end{align}
とモデル化します。$\boldsymbol{w}$は係数ベクトルパラメータ、$\sigma(x)$はシグモイド関数で、${\rm Bern}(y;\pi)$は確率$\pi$で$y=1$、確率$1-\pi$で$y=0$となるベルヌーイ分布：
\begin{align}
{\rm Bern}(y;\pi) &= \pi^y (1-\pi)^{1-y}
\end{align}
です。
事後分布のラプラス近似
学習データを$\mathcal{D}_N \equiv {(\boldsymbol{x}_n,y_n)}_{n=1}^N$とすると対数尤度は
\begin{align}
\ln p(\mathcal{D}_N|\boldsymbol{w})
&= \ln \prod_{n=1}^N {\rm Bern}(y_n; \sigma(\boldsymbol{x}_n^T \boldsymbol{w})) \\
&= \sum_{n=1}^N \left\{
y_n \ln \sigma(\boldsymbol{x}_n^T \boldsymbol{w}) + (1-y_n) \ln (1 - \sigma(\boldsymbol{x}_n^T \boldsymbol{w}))
\right\}
\end{align}
です。対数を取ったときに$\boldsymbol{w}$についてこのような関数形になる（名前の付いた）分布はありませんね。。よって予測分布が解析的に求まるとか、既存のサンプラーを使うことのできるような共役事前分布はありません。
そこでここでは事後分布を正規分布で近似（ラプラス近似）することにします。事前分布も簡単のため正規分布とします。
\begin{align}
p(\boldsymbol{w}) &= \mathcal{N}(\boldsymbol{w};\boldsymbol{\mu}_0,\boldsymbol{\Sigma}_0)
\end{align}
事後分布の対数は
\begin{align}
\ln p(\boldsymbol{w}|\mathcal{D}_N)
&= \ln p(\mathcal{D}_N|\boldsymbol{w}) + \ln p(\boldsymbol{w}) + {\rm const.} \\
&= \ln \prod_{n=1}^N {\rm Bern}(y_n; \sigma(\boldsymbol{x}_n^T \boldsymbol{w})) + \ln \mathcal{N}(\boldsymbol{w};\boldsymbol{\mu}_0,\boldsymbol{\Sigma}_0) +{\rm const.} \\
&= \sum_{n=1}^N \left\{
y_n \ln \sigma(\boldsymbol{x}_n^T \boldsymbol{w}) + (1-y_n) \ln (1 - \sigma(\boldsymbol{x}_n^T \boldsymbol{w}))
\right\} \\
&\hspace{4mm}
- \frac{1}{2} (\boldsymbol{w} - \boldsymbol{\mu}_0)^T \boldsymbol{\Sigma}_0^{-1} (\boldsymbol{w} - \boldsymbol{\mu}_0) + {\rm const.}
\end{align}
です。分布の確率密度関数の値が一番大きくなる値をモード（最頻値）と言います。事後分布のモードは事後確率を最大化する値であり、MAP（Maximum A Posteriori）推定値と言います。事後分布のモードを$\hat{\boldsymbol{w}}$として、事後分布の対数を$\hat{\boldsymbol{w}}$の周りで二次近似すると、モードでは一階微分は$\boldsymbol{0}$なので
\begin{align}
\ln p(\boldsymbol{w}|\mathcal{D}_N)
&\approx \ln p(\hat{\boldsymbol{w}}|\mathcal{D}_N) + \frac{1}{2} (\boldsymbol{w}-\hat{\boldsymbol{w}})^T \left(\frac{\partial^2}{\partial \boldsymbol{w} \boldsymbol{w}^T} \ln p(\boldsymbol{w}|\mathcal{D}_N) \bigg|_{\boldsymbol{w}=\hat{\boldsymbol{w}}} \right) (\boldsymbol{w}-\hat{\boldsymbol{w}}) \\
&= \ln p(\hat{\boldsymbol{w}}|\mathcal{D}_N) + \frac{1}{2} (\boldsymbol{w}-\hat{\boldsymbol{w}})^T \boldsymbol{H}(\hat{\boldsymbol{w}}) (\boldsymbol{w}-\hat{\boldsymbol{w}}) \\
\boldsymbol{H}(\hat{\boldsymbol{w}})
&\equiv \frac{\partial^2}{\partial \boldsymbol{w} \boldsymbol{w}^T} \ln p(\boldsymbol{w}|\mathcal{D}_N) \bigg|_{\boldsymbol{w}=\hat{\boldsymbol{w}}}
\end{align}
となります。これは$\boldsymbol{w}$についての二次形式なので、事後分布を正規分布で近似していることになります。
\begin{align}
p(\boldsymbol{w}|\mathcal{D}_N)
&\approx \mathcal{N}(\boldsymbol{w};\boldsymbol{\mu}_N, \boldsymbol{\Sigma}_N) \\
\boldsymbol{\mu}_N
&\equiv \hat{\boldsymbol{w}} \\
\boldsymbol{\Sigma}_N
&\equiv -\boldsymbol{H}(\hat{\boldsymbol{w}})^{-1}
\end{align}
このように事後分布を、モード周りで正規分布で近似することをラプラス近似と呼びます。$\boldsymbol{H}(\boldsymbol{w})$のような二階導関数行列をヘッセ行列と呼びます。$\boldsymbol{H}(\hat{\boldsymbol{w}})$は「事後分布の対数」のヘッセ行列の、モード$\hat{\boldsymbol{w}}$での値です。
ニュートン法による学習
次にモード$\hat{\boldsymbol{w}}$の求め方です。勾配法では
\begin{align}
\hat{\boldsymbol{w}}^{(t+1)} &\leftarrow \hat{\boldsymbol{w}}^{(t)} + \alpha \boldsymbol{g}(\hat{\boldsymbol{w}}^{(t)}) \\
\boldsymbol{g}(\boldsymbol{w}) &\equiv \frac{\partial}{\partial \boldsymbol{w}} \ln p(\boldsymbol{w}|\mathcal{D}_N)
\end{align}
と更新します（ここでは最大化のため$\alpha \boldsymbol{g}(\boldsymbol{w})$の符号は$+$になっています）。このとき$\alpha$は学習率で、勾配法ではパラメータとなります。
勾配法は勾配が$\boldsymbol{0}$になるところで収束します。関数値が$0$となるところを求める高速なアルゴリズムにニュートン法があります。簡単のため一変数関数で考えると、関数$f(x)$の$x^{(t)}$での一次近似は
\begin{align}
f(x) = f(x^{(t)}) + f^\prime(x^{(t)}) (x - x^{(t)})
\end{align}
であり、$f(x)=0$を解くと
\begin{align}
x = x^{(t)} - \frac{f(x^{(t)})}{f^\prime(x^{(t)})}
\end{align}
です。この更新を反復して、$f(x)=0$の解を求めるのがニュートン法です。
今回の場合は勾配$\boldsymbol{g}(\boldsymbol{w})$が$\boldsymbol{0}$となるところを求めるので、勾配の微分、すなわち目的関数（事後分布の対数）の二階微分$\boldsymbol{H}(\boldsymbol{w})$が出てきて、ニュートン法の更新式は
\begin{align}
\hat{\boldsymbol{w}}^{(t+1)} &\leftarrow \hat{\boldsymbol{w}}^{(t)} - \boldsymbol{H}(\hat{\boldsymbol{w}}^{(t)})^{-1} \boldsymbol{g}(\hat{\boldsymbol{w}}^{(t)})
\end{align}
となります。ニュートン法が収束したとき、ニュートン法の更新で計算しているヘッセ行列が、そのまま事後分布の共分散行列を求めるのに利用できます。
ニュートン法とラプラス近似は相性の良い組み合わせですね。
次回は「ベイズ推定６：周辺尤度」です。