正規線形回帰モデルのベイズ推定

正規線形回帰モデル

特徴量（説明変数）ベクトル$\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}

です。

次回は「ベイズ推定５：ロジスティック回帰のラプラス近似によるベイズ推定」です。