はじめに
『パターン認識と機械学習』の独学時のまとめです。一連の記事は「数式の行間埋め」または「R・Pythonでのスクラッチ実装」からアルゴリズムの理解を補助することを目的としています。本とあわせて読んでください。
この記事は、4.5.1項の内容です。ロジスティック回帰に事前分布を導入したベイズロジスティック回帰のMAP推定と近似事後分布を導出します。
【実装編】
www.anarchive-beta.com
【前節の内容】
www.anarchive-beta.com
【他の節一覧】
www.anarchive-beta.com
【この節の内容】
4.5.1 ラプラス近似
ベイズロジスティック回帰のパラメータの事後分布は解析的に求められません。そこで、事後分布の近似分布を求めます。
まずは、パラメータのMAP解(最大事後確率解)を求めます。そして、MAP解における事後分布に対して、ラプラス近似(ガウス分布による近似)による近似分布を導出します。ロジスティック回帰については4.3.2-3項、ラプラス近似については4.4.0項を参照してください。
・事後分布の導出
ベイズロジスティック回帰のパラメータの事後分布を求めます。ロジスティック回帰については「4.3.2-3：ロジスティック回帰の導出【PRMLのノート】 - からっぽのしょこ」を参照してください。
4.3.2項と同様に、基底関数$\boldsymbol{\phi}(\mathbf{x}) = (\phi_0(\mathbf{x}), \phi_1(\mathbf{x}), \cdots, \phi_{M-1}(\mathbf{x}))$によって変換した入力を$\boldsymbol{\phi}_n = \boldsymbol{\phi}(\mathbf{x}_n)$で表し、重みパラメータ$\mathbf{w} = (w_0, w_1, \cdots, w_{M-1})^{\top}$との重み付き和をシグモイド関数$\sigma(\cdot)$で変換します。
$$
y_n = \sigma(\mathbf{w}^{\top} \boldsymbol{\phi}_n)
$$
尤度関数を、$\mathbf{y} = (y_1, \cdots, y_N)^{\top}$をパラメータとするベルヌーイ分布の積とします。
$$
p(\mathbf{t} | \mathbf{w})
= \prod_{n=1}^N
\mathrm{Bern}(t_n | \mathbf{w})
= \prod_{n=1}^N
y_n^{t_n} (1 - y_n)^{1-t_n}
\tag{4.89}
$$
$\mathbf{w}$の事前分布を、平均$\mathbf{m}_0$・分散共分散行列$\mathbf{S}_0$の$M$次元ガウス分布とします。
$$
p(\mathbf{w})
= \mathcal{N}(\mathbf{w} | \mathbf{m}_0, \mathbf{S}_0)
= \frac{1}{(2 \pi)^{\frac{M}{2}}}
\frac{1}{|\mathbf{S}_0|^{\frac{1}{2}}}
\exp \left\{
- \frac{1}{2}
(\mathbf{w} - \mathbf{m}_0)^{\top}
\mathbf{S}_0^{-1}
(\mathbf{w} - \mathbf{m}_0)
\right\}
\tag{4.140}
$$
$\mathbf{w}$の事後分布は、ベイズの定理より
$$
\begin{align}
p(\mathbf{w} | \mathbf{t})
&= \frac{p(\mathbf{t} | \mathbf{w}) p(\mathbf{w})}{p(\mathbf{t})}
\\
&\propto
p(\mathbf{t} | \mathbf{w}) p(\mathbf{w})
\tag{4.141}
\end{align}
$$
で求められます。周辺分布(分母)は$\mathbf{w}$に影響しないため省略して、$\mathbf{w}$に関する比例関係に注目します。
式を分かりやすくするために事後分布の対数をとり、また具体的な式を代入します。
$$
\begin{aligned}
\ln p(\mathbf{w} | \mathbf{t})
&= \ln p(\mathbf{t} | \mathbf{w})
+ \ln p(\mathbf{w})
+ \mathrm{const.}
\\
&= \sum_{n=1}^N
\ln \mathrm{Bern}(t_n | \mathbf{w})
+ \ln \mathcal{N}(\mathbf{w} | \mathbf{m}_0, \mathbf{S}_0)
+ \mathrm{const.}
\\
&= \sum_{n=1}^N \Bigl\{
t_n \ln y_n
+ (1 - t_n) \ln (1 - y_n)
\Bigr\} \\
&\qquad
- \frac{M}{2} \ln (2 \pi)
- \frac{1}{2} \ln |\mathbf{S}_0|
- \frac{1}{2}
(\mathbf{w} - \mathbf{m}_0)^{\top}
\mathbf{S}_0^{-1}
(\mathbf{w} - \mathbf{m}_0)
+ \mathrm{const.}
\end{aligned}
$$
$\mathbf{w}$と無関係な周辺分布を$\mathrm{const.}$とおきました。後で微分する際に、$\mathrm{const.}$は0になりに消えます。
括弧を展開して、$\mathbf{w}$と無関係な項を$\mathrm{const.}$にまとめて式を整理します。
$$
\begin{align}
\ln p(\mathbf{w} | \mathbf{t})
&= \sum_{n=1}^N \Bigl\{
t_n \ln y_n
+ (1 - t_n) \ln (1 - y_n)
\Bigr\}
- \frac{1}{2}
(\mathbf{w} - \mathbf{m}_0)^{\top}
\mathbf{S}_0^{-1}
(\mathbf{w} - \mathbf{m}_0)
+ \mathrm{const.}
\\
&= \sum_{n=1}^N \Bigl\{
t_n \ln y_n
+ (1 - t_n) \ln (1 - y_n)
\Bigr\} \\
&\qquad
- \frac{1}{2} \Bigl\{
\mathbf{w}^{\top} \mathbf{S}_0^{-1} \mathbf{w}
- \mathbf{w}^{\top} \mathbf{S}_0^{-1} \mathbf{m}_0
- \mathbf{m}_0^{\top} \mathbf{S}_0^{-1} \mathbf{w}
+ \mathbf{m}_0^{\top} \mathbf{S}_0^{-1} \mathbf{m}_0
\Bigr\}
+ \mathrm{const.}
\\
&= \sum_{n=1}^N \Bigl\{
t_n \ln y_n
+ (1 - t_n) \ln (1 - y_n)
\Bigr\}
- \frac{1}{2}
\mathbf{w}^{\top} \mathbf{S}_0^{-1} \mathbf{w}
+ \mathbf{w}^{\top} \mathbf{S}_0^{-1} \mathbf{m}_0
+ \mathrm{const.}
\end{align}
$$
$\mathbf{w}^{\top} \mathbf{S}_0^{-1} \mathbf{m}_0$はスカラなので、$\mathbf{w}^{\top} \mathbf{S}_0^{-1} \mathbf{m}_0 = (\mathbf{w}^{\top} \mathbf{S}_0^{-1} \mathbf{m}_0)^{\top} = \mathbf{m}_0^{\top} (\mathbf{S}_0^{-1})^{\top} \mathbf{w}$と転置できます。また、精度行列$\mathbf{S}_0^{-1}$は対称行列なので、$\mathbf{S}_0^{-1} = (\mathbf{S}_0^{-1})^{\top}$です。
以上で、対数をとった$\mathbf{w}$の事後分布の式が得られました。
・勾配の導出
対数をとった事後分布の勾配(1階微分)を求めます。
$\mathbf{w}$に関して$\mathbf{w}$の対数事後分布を微分します。
$$
\begin{aligned}
\frac{\partial \ln p(\mathbf{w} | \mathbf{t})}{\partial \mathbf{w}}
&= \frac{\partial}{\partial \mathbf{w}} \left[
\sum_{n=1}^N \Bigl\{
t_n \ln y_n
+ (1 - t_n) \ln (1 - y_n)
\Bigr\}
- \frac{1}{2}
\mathbf{w}^{\top} \mathbf{S}_0^{-1} \mathbf{w}
+ \mathbf{w}^{\top} \mathbf{S}_0^{-1} \mathbf{m}_0
+ \mathrm{const.}
\right]
\\
&= \sum_{n=1}^N
\frac{\partial}{\partial \mathbf{w}} \Bigl\{
t_n \ln y_n
+ (1 - t_n) \ln (1 - y_n)
\Bigr\}
- \frac{1}{2}
\frac{\partial \mathbf{w}^{\top} \mathbf{S}_0^{-1} \mathbf{w}}{\partial \mathbf{w}}
+ \frac{\partial \mathbf{w}^{\top} \mathbf{S}_0^{-1} \mathbf{m}_0}{\partial \mathbf{w}}
\end{aligned}
$$
1つ目の項は、ロジスティック回帰の対数尤度関数の微分です。4.3.2項で求めた負の対数尤度関数の微分(4.91)の符号を反転させた(-1を掛けた)
$$
\frac{\partial \ln p(t_n | \sigma(\mathbf{w}^{\top} \boldsymbol{\phi}_n))}{\partial \mathbf{w}}
= - (y_n - t_n) \boldsymbol{\phi}_n
\tag{4.91'}
$$
となります。
2つ目の項に関して、ベクトル$\mathbf{x}$と行列$\mathbf{A}$の積の微分$\frac{\partial \mathbf{x}^{\top} \mathbf{A} \mathbf{x}}{\partial \mathbf{x}} = (\mathbf{A} + \mathbf{A}^{\top}) \mathbf{x}$より
$$
\frac{\partial \mathbf{w}^{\top} \mathbf{S}_0^{-1} \mathbf{w}}{\partial \mathbf{w}}
= \Bigl\{
\mathbf{S}_0^{-1} + (\mathbf{S}_0^{-1})^{\top}
\Bigr\}
\mathbf{w}
= 2 \mathbf{S}_0^{-1} \mathbf{w}
$$
となります。
3つ目の項は、$\mathbf{S}_0^{-1} \mathbf{m}_0$を1つの行列とみなして、ベクトルと行列の積の微分$\frac{\partial \mathbf{x}^{\top} \mathbf{A}}{\partial \mathbf{x}} = \mathbf{A}$より
$$
\frac{\partial \mathbf{w}^{\top} \mathbf{S}_0^{-1} \mathbf{m}_0}{\partial \mathbf{w}}
= \mathbf{S}_0^{-1} \mathbf{m}_0
$$
となります。
それぞれ代入します。
$$
\begin{aligned}
\frac{\partial \ln p(\mathbf{w} | \mathbf{t})}{\partial \mathbf{w}}
&= - \sum_{n=1}^N
(y_n - t_n) \boldsymbol{\phi}_n
- \mathbf{S}_0^{-1} \mathbf{w}
+ \mathbf{S}_0^{-1} \mathbf{m}_0
\\
&= - \sum_{n=1}^N
(y_n - t_n) \boldsymbol{\phi}_n
- \mathbf{S}_0^{-1} (
\mathbf{w} - \mathbf{m}_0
)
= - \boldsymbol{\Phi}^{\top} (\mathbf{y} - \mathbf{t})
- \mathbf{S}_0^{-1} (
\mathbf{w} - \mathbf{m}_0
)
\equiv
\nabla \ln p(\mathbf{w} | \mathbf{t})
\end{aligned}
$$
最後の変形は4.3.3項で確認しました。
以上で、対数をとった$\mathbf{w}$の事後分布の勾配が得られました。
・ヘッセ行列の導出
対数をとった事後分布のヘッセ行列(2階微分)を求めます。
$\mathbf{w}$に関して$\mathbf{w}$の対数事後分布の勾配を更に微分します。
$$
\begin{aligned}
\frac{\partial \ln p(\mathbf{w} | \mathbf{t})}{\partial \mathbf{w} \partial \mathbf{w}^{\top}}
&= \frac{\partial}{\partial \mathbf{w}}
\frac{\partial \ln p(\mathbf{w} | \mathbf{t})}{\partial \mathbf{w}}
\\
&= \frac{\partial}{\partial \mathbf{w}} \left\{
- \sum_{n=1}^N
(y_n - t_n) \boldsymbol{\phi}_n
- \mathbf{S}_0^{-1} (
\mathbf{w} - \mathbf{m}_0
)
\right\}
\\
&= - \sum_{n=1}^N
\frac{\partial}{\partial \mathbf{w}}
(y_n - t_n) \boldsymbol{\phi}_n
- \frac{\partial \mathbf{S}_0^{-1} \mathbf{w}}{\partial \mathbf{w}}
+ \frac{\partial \mathbf{S}_0^{-1} \mathbf{m}_0}{\partial \mathbf{w}}
\end{aligned}
$$
1つ目の項は、ロジスティック回帰の対数尤度関数の2階微分です。4.3.3項で求めた負の対数尤度関数の2階微分(4.97)の符号を反転させた(-1を掛けた)
$$
\frac{
\partial^2 \ln p(t_n | \sigma(\mathbf{w}^{\top} \boldsymbol{\phi}_n))
}{
\partial \mathbf{w} \partial \mathbf{w}^{\top}
}
= - \sum_{n=1}^N
y_n (1 - y_n)
\boldsymbol{\phi}_n \boldsymbol{\phi}_n^{\top}
= - \boldsymbol{\Phi}^{\top} \mathbf{R} \boldsymbol{\Phi}
\tag{4.97'}
$$
となります。$\mathbf{R}$は、対角要素が$R_{n,n} = y_n (1 - y_n)$の対角行列です。
2つ目の項は、行列とベクトルの積の微分$\frac{\partial \mathbf{A} \mathbf{x}}{\partial \mathbf{x}} = \mathbf{A}^{\top}$より
$$
\frac{\partial \mathbf{S}_0^{-1} \mathbf{w}}{\partial \mathbf{w}}
= (\mathbf{S}_0^{-1})^{\top}
= \mathbf{S}_0^{-1}
$$
となります。
3つ目の項は、$\mathbf{w}$を含まないため0行列となり消えます。
それぞれ代入します。
$$
\frac{\partial \ln p(\mathbf{w} | \mathbf{t})}{\partial \mathbf{w} \partial \mathbf{w}^{\top}}
= - \sum_{n=1}^N
y_n (1 - y_n)
\boldsymbol{\phi}_n \boldsymbol{\phi}_n^{\top}
- \mathbf{S}_0^{-1}
= - \boldsymbol{\Phi}^{\top} \mathbf{R} \boldsymbol{\Phi}
- \mathbf{S}_0^{-1}
\equiv
\nabla \nabla \ln p(\mathbf{w} | \mathbf{t})
$$
ベイズロジスティック回帰の負の(符号を反転させた)対数事後分布のヘッセ行列を$\mathbf{S}_N^{-1}$とおきます。
$$
\mathbf{S}_N^{-1}
\equiv
- \nabla \nabla \ln p(\mathbf{w} | \mathbf{t})
= \mathbf{S}_0^{-1}
+ \sum_{n=1}^N
y_n (1 - y_n)
\boldsymbol{\phi}_n \boldsymbol{\phi}_n^{\top}
= \mathbf{S}_0^{-1}
+ \boldsymbol{\Phi}^{\top} \mathbf{R} \boldsymbol{\Phi}
\tag{4.143}
$$
以上で、対数をとった$\mathbf{w}$の事後分布のヘッセ行列が得られました。
・パラメータの更新式の導出
ニュートン-ラフソン法によって、ベイズロジスティック回帰のパラメータのMAP解(最大事後確率解)を求めます。ニュートン-ラフソン法については「4.3.2-3：ロジスティック回帰の導出【PRMLのノート】 - からっぽのしょこ」を参照してください。
ニュートン-ラフソン法の更新式
$$
\mathbf{w}^{(\mathrm{new})}
= \mathbf{w}^{(\mathrm{old})}
- \mathbf{H}^{-1}
\nabla E(\mathbf{w})
\tag{4.92}
$$
に、負の対数事後分布$- \ln p(\mathbf{w} | \mathbf{t})$のヘッセ行列$\mathbf{H} = \mathbf{S}_N^{-1} = - \nabla \nabla \ln p(\mathbf{w} | \mathbf{t})$、勾配$\nabla E(\mathbf{w}) = - \nabla \ln p(\mathbf{w} | \mathbf{t})$を代入します。
$$
\begin{aligned}
\mathbf{w}^{(\mathrm{new})}
&= \mathbf{w}^{(\mathrm{old})}
- (\mathbf{S}_N^{-1})^{-1}
(- \nabla \ln p(\mathbf{w}^{(\mathrm{old})} | \mathbf{t}))
\\
&= \mathbf{w}^{(\mathrm{old})}
- \Bigl(
\mathbf{S}_0^{-1}
+ \boldsymbol{\Phi}^{\top} \mathbf{R} \boldsymbol{\Phi}
\Bigr)^{-1}
\Bigl\{
\boldsymbol{\Phi}^{\top} (\mathbf{y} - \mathbf{t})
+ \mathbf{S}_0^{-1} (
\mathbf{w}^{(\mathrm{old})} - \mathbf{m}_0
)
\Bigr\}
\end{aligned}
$$
以上で、ベイズロジスティック回帰の重みパラメータの更新式が得られました。(4.3.3項のときのように変形した方がいいのかもしれないけどよく分からない。というかなぜあんな変形をしたんだ？)
ロジスティック回帰の更新式(4.99)(の変形前)と比較すると、事前分布のパラメータ$\mathbf{m}_0, \mathbf{S}_0$によって調整されているのが分かります。
繰り返しパラメータを更新して、収束した値(負の対数事後確率密度が最小となる値)をMAP解$\mathbf{w}_{\mathrm{MAP}}$とします。
・近似事後分布の導出
ラプラス近似によって、事後分布の近似分布を求めます。ラプラス近似については「4.4.0：ラプラス近似の導出【PRMLのノート】 - からっぽのしょこ」を参照してください。
$\mathbf{w}$の対数事後分布$\ln p(\mathbf{w} | \mathbf{t})$をMAP解$\mathbf{w}_{\mathrm{MAP}}$の周りで2次までテイラー展開します。
$$
\begin{aligned}
\ln p(\mathbf{w} | \mathbf{t})
&\simeq
\ln p(\mathbf{w}_{\mathrm{MAP}} | \mathbf{t})
+ (\mathbf{w} - \mathbf{w}_{\mathrm{MAP}})^{\top}
\nabla \ln p(\mathbf{w} | \mathbf{t})|_{\mathbf{w}=\mathbf{w}_{\mathrm{MAP}}} \\
&\qquad
+ \frac{1}{2}
(\mathbf{w} - \mathbf{w}_{\mathrm{MAP}})^{\top}
\nabla \nabla \ln p(\mathbf{w} | \mathbf{t})|_{\mathbf{w}=\mathbf{w}_{\mathrm{MAP}}}
(\mathbf{w} - \mathbf{w}_{\mathrm{MAP}})
\end{aligned}
$$
2番目の項に関して、$\mathbf{w}_{\mathrm{MAP}}$は負の対数事後確率密度が最小(事後確率密度が最大)となる値なので
$$
\nabla \ln p(\mathbf{w} | \mathbf{t})|_{\mathbf{w}=\mathbf{w}_{\mathrm{MAP}}}
= \left.
\frac{\partial \ln p(\mathbf{w} | \mathbf{t})}{\partial \mathbf{w}}
\right|_{\mathbf{w}=\mathbf{w}_{\mathrm{MAP}}}
= \mathbf{0}
$$
となります(MAP解は極小値(極大値)なので微分すると0ベクトルになります)。
また、3番目の項に関して
$$
\mathbf{S}_N^{-1}
= - \nabla \nabla \ln p(\mathbf{w} | \mathbf{t})
\tag{4.143}
$$
で置き換えると、$\ln p(\mathbf{w} | \mathbf{t})$は
$$
\ln p(\mathbf{w} | \mathbf{t})
\simeq
\ln p(\mathbf{w}_{\mathrm{MAP}} | \mathbf{t})
- \frac{1}{2}
(\mathbf{w} - \mathbf{w}_{\mathrm{MAP}})^{\top}
\mathbf{S}_N^{-1}
(\mathbf{w} - \mathbf{w}_{\mathrm{MAP}})
$$
で近似できます。
近似式の両辺の指数をとり、$\ln$を外します。
$$
p(\mathbf{w} | \mathbf{t})
\simeq
p(\mathbf{w}_{\mathrm{MAP}} | \mathbf{t})
\exp \left\{
- \frac{1}{2}
(\mathbf{w} - \mathbf{w}_{\mathrm{MAP}})^{\top}
\mathbf{S}_N^{-1}
(\mathbf{w} - \mathbf{w}_{\mathrm{MAP}})
\right\}
$$
この式は、平均$\mathbf{w}_{\mathrm{MAP}}$・分散共分散行列$\mathbf{S}_N$の$M$次元ガウス分布の形をしています。よって、多変量ガウス分布の式(2.43)と対応させると正規化係数が分かります。
$$
q(\mathbf{w} | \mathbf{t})
\equiv
\frac{1}{(2 \pi)^{\frac{M}{2}}}
\frac{1}{|\mathbf{S}_N|^{\frac{1}{2}}}
\exp \left\{
- \frac{1}{2}
(\mathbf{w} - \mathbf{w}_{\mathrm{MAP}})^{\top}
\mathbf{S}_N^{-1}
(\mathbf{w} - \mathbf{w}_{\mathrm{MAP}})
\right\}
= \mathcal{N}(\mathbf{w} | \mathbf{w}_{\mathrm{MAP}}, \mathbf{S}_N)
$$
以上で、MAP解における$\mathbf{w}$の近似事後分布$q(\mathbf{w} | \mathbf{t})$が得られました。
平均$\mathbf{m}_0$・分散共分散行列$\mathbf{S}_0$のガウス分布を事前分布として、近似事後分布を求めました。このことは、「事前分布$p(\mathbf{w})$が尤度関数$p(\mathbf{t} | \mathbf{w})$によって(観測データ$\mathbf{t}$によって)近似事後分布$q(\mathbf{w} | \mathbf{t})$に更新された」と言えます。
この近似事後分布を事前分布として用いて、新たなMAP解と近似事後分布を求めます。これを収束するまで繰り返すことで、解析的に求められない真の事後分布に近付けることができます(たぶん)。
参考文献
- C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上下』,丸善出版,2012年.
おわりに
なんで最後まで本に載ってないのー。
【次節の内容】
続かずにたぶん戻る
