正規線形回帰モデルのベイズ推定
正規線形回帰モデル
特徴量（説明変数）ベクトル$\boldsymbol{x}_n$に対する目的変数値$y_n \in \mathbb{R}$の組が$N$個得られているとします。これを学習データ$\mathcal{D}_N \equiv \{(\boldsymbol{x}_n,y_n)\}_{n=1}^N$として、任意の$\boldsymbol{x}$に対する予測$\hat{f}(\boldsymbol{x})$を作ることを考えます。
線形回帰モデルは関数$f(\boldsymbol{x})$を
\begin{align}
f(\boldsymbol{x};\boldsymbol{w}) = \boldsymbol{x}^T \boldsymbol{w}
\end{align}
とするモデルです。$\boldsymbol{w}$は係数ベクトルパラメータです。 正規線形回帰モデルはこの線形回帰に、正規分布に従うノイズ$\epsilon$が乗ったとするモデルです。
\begin{align}
y_n = \boldsymbol{x}_n^T \boldsymbol{w} + \epsilon_n
\end{align}
確率分布で書くと
\begin{align}
p(y|\boldsymbol{x},\boldsymbol{w}) &= \mathcal{N}(y; \boldsymbol{x}^T \boldsymbol{w}, \lambda^{-1})
\end{align}
です。
このモデルについてベイズ推定してみます。
正規線形回帰モデルのベイズ推定
やることは正規分布の平均についてのベイズ推定のときと同じです。
対数尤度関数は
\begin{align}
\ln \prod_{n=1}^N p(y_n|\boldsymbol{x}_n,\boldsymbol{w})
&= \sum_{n=1}^N \ln \mathcal{N}(y_n;\boldsymbol{x}_n^T \boldsymbol{w}, \lambda^{-1}) \\
&= -\frac{\lambda}{2} \sum_{n=1}^N (y_n - \boldsymbol{x}_n^T \boldsymbol{w})^2 + {\rm const.}
\end{align}
で$\boldsymbol{w}$についての二次形式なので、共役事前分布は多変量正規分布です。よって、$\boldsymbol{w}$の事前分布を
\begin{align}
p(\boldsymbol{w}) &= \mathcal{N}(\boldsymbol{w};\boldsymbol{\mu}_0,\boldsymbol{\Sigma}_0)
\end{align}
とすれば、事後分布も正規分布：
\begin{align}
p(\boldsymbol{w}|\mathcal{D}_N)
&= \mathcal{N}(\boldsymbol{w};\boldsymbol{\mu}_N,\boldsymbol{\Sigma}_N)
\tag{1}
\end{align}
となります。事後分布の対数は
\begin{align}
&\ln p(\boldsymbol{w}|\mathcal{D}_N) \\
&= \ln \prod_{n=1}^N p(y_n|\boldsymbol{x}_n,\boldsymbol{w}) + \ln p(\boldsymbol{w}) + {\rm const.} \\
&= \sum_{n=1}^N \ln \mathcal{N}(y_n;\boldsymbol{x}_n^T \boldsymbol{w},\lambda^{-1})
+ \ln \mathcal{N}(\boldsymbol{w};\boldsymbol{\mu}_0,\boldsymbol{\Sigma}_0)
+ {\rm const.} \\
&= - \frac{\lambda}{2} \sum_{n=1}^N (y_n - \boldsymbol{x}_n^T \boldsymbol{w})^2
- \frac{1}{2} (\boldsymbol{w} - \boldsymbol{\mu}_0)^T \boldsymbol{\Sigma}_0^{-1} (\boldsymbol{w} - \boldsymbol{\mu}_0)
+ {\rm const.} \\
&= -\frac{1}{2} \left\{
\boldsymbol{w}^T \left(
\lambda \sum_{n=1}^N \boldsymbol{x}_n \boldsymbol{x}_n^T + \boldsymbol{\Sigma}_0^{-1}
\right) \boldsymbol{w}
- 2 \left(
\lambda \sum_{n=1}^N y_n \boldsymbol{x}_n + \boldsymbol{\Sigma}_0^{-1} \boldsymbol{\mu}_0
\right)^T \boldsymbol{w}
\right\} + {\rm const.}
\end{align}
です。あるいは計画行列：
\begin{align}
\boldsymbol{x}_N &\equiv \left(
\begin{array}{l}
\boldsymbol{x}_1^T \\
\vdots \\
\boldsymbol{x}_N^T
\end{array}
\right)
\end{align}
と$\boldsymbol{y}_N = (y_1,\cdots,y_N)^T$を用いると
\begin{align}
&-\frac{\lambda}{2} \sum_{n=1}^N (y_n - \boldsymbol{x}_n^T \boldsymbol{w})^2
- \frac{1}{2} (\boldsymbol{w} - \boldsymbol{\mu}_0)^T \boldsymbol{\Sigma}_0^{-1} (\boldsymbol{w} - \boldsymbol{\mu}_0)
+ {\rm const.} \\
&= -\frac{\lambda}{2} \parallel \boldsymbol{y}_N - \boldsymbol{x}_N \boldsymbol{w} \parallel^2
- \frac{1}{2} (\boldsymbol{w} - \boldsymbol{\mu}_0)^T \boldsymbol{\Sigma}_0^{-1} (\boldsymbol{w} - \boldsymbol{\mu}_0)
+ {\rm const.} \\
&= -\frac{1}{2} \left\{
\boldsymbol{w}^T \left(
\lambda \boldsymbol{x}_N^T \boldsymbol{x}_N + \boldsymbol{\Sigma}_0^{-1}
\right) \boldsymbol{w}
- 2 \left(
\lambda \boldsymbol{x}_N^T \boldsymbol{y}_N + \boldsymbol{\Sigma}_0^{-1} \boldsymbol{\mu}_0
\right)^T \boldsymbol{w}
\right\} + {\rm const.}
\end{align}
と書けます。(1)式の対数：
\begin{align}
-\frac{1}{2} \left\{
\boldsymbol{w}^T \boldsymbol{\Sigma}_N^{-1} \boldsymbol{w} - 2 \boldsymbol{\mu}_N^T \boldsymbol{\Sigma}_N^{-1} \boldsymbol{w}
\right\} + {\rm const.}
\end{align}
と係数を対応させて
\begin{align}
\boldsymbol{\Sigma}_N
&= \left(
\lambda \boldsymbol{x}_N^T \boldsymbol{x}_N + \boldsymbol{\Sigma}_0^{-1}
\right)^{-1} \\
\boldsymbol{\mu}_N
&= \boldsymbol{\Sigma}_N \left(
\lambda \boldsymbol{x}_N \boldsymbol{y}_N + \boldsymbol{\Sigma}_0^{-1} \boldsymbol{\mu}_0
\right) \\
&=\left(
\lambda \boldsymbol{x}_N^T \boldsymbol{x}_N + \boldsymbol{\Sigma}_0^{-1}
\right)^{-1}
\left(
\lambda \boldsymbol{x}_N^T \boldsymbol{y}_N + \boldsymbol{\Sigma}_0^{-1} \boldsymbol{\mu}_0
\right)
\end{align}
を得ます。
以上より$\boldsymbol{w}$の事後分布は
\begin{align}
p(\boldsymbol{w}|\mathcal{D}_N)
&= \mathcal{N}(\boldsymbol{w};\boldsymbol{\mu}_N,\boldsymbol{\Sigma}_N) \\
\boldsymbol{\mu}_N
&=\left(
\lambda \boldsymbol{x}_N^T \boldsymbol{x}_N + \boldsymbol{\Sigma}_0^{-1}
\right)^{-1}
\left(
\lambda \boldsymbol{x}_N^T \boldsymbol{y}_N + \boldsymbol{\Sigma}_0^{-1} \boldsymbol{\mu}_0
\right) \\
\boldsymbol{\Sigma}_N
&= \left(
\lambda \boldsymbol{x}_N^T \boldsymbol{x}_N + \boldsymbol{\Sigma}_0^{-1}
\right)^{-1}
\end{align}
と求まりました。$\boldsymbol{\Sigma}_0$を対角行列とし、${\rm diag}(\boldsymbol{\Sigma}_0) \rightarrow {\bf \infty}$とすると
\begin{align}
\boldsymbol{\mu}_N \rightarrow (\boldsymbol{x}_N^T \boldsymbol{x}_N)^{-1} \boldsymbol{x}_N^T \boldsymbol{y}_N
\end{align}
となり、最小二乗解となります。
予測分布は
\begin{align}
p(Y|\boldsymbol{x},\mathcal{D}_N)
&= \int p(Y|\boldsymbol{x},\boldsymbol{w}) p(\boldsymbol{w}|\mathcal{D}_N) d\boldsymbol{w} \\
&= \int \mathcal{N}(Y;\boldsymbol{x}^T \boldsymbol{w},\lambda^{-1}) \mathcal{N}(\boldsymbol{w};\boldsymbol{\mu}_N,\boldsymbol{\Sigma}_N) d\boldsymbol{w} \\
&= \mathcal{N}(Y;\boldsymbol{x}^T \boldsymbol{\mu}_N,\boldsymbol{x}^T \boldsymbol{\Sigma}_N \boldsymbol{x} + \lambda^{-1})
\end{align}
です。