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

# 正規線形回帰モデルのベイズ推定
## 正規線形回帰モデル
特徴量（説明変数）ベクトル$\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})$を

```math
\begin{align}
f(\boldsymbol{x};\boldsymbol{w}) = \boldsymbol{x}^T \boldsymbol{w}
\end{align}
```
とするモデルです。$\boldsymbol{w}$は係数ベクトルパラメータです。 正規線形回帰モデルはこの線形回帰に、正規分布に従うノイズ$\epsilon$が乗ったとするモデルです。

```math
\begin{align}
y_n = \boldsymbol{x}_n^T \boldsymbol{w} + \epsilon_n
\end{align}
```
確率分布で書くと

```math
\begin{align}
p(y|\boldsymbol{x},\boldsymbol{w}) &= \mathcal{N}(y; \boldsymbol{x}^T \boldsymbol{w}, \lambda^{-1})
\end{align}
```
です。 

このモデルについてベイズ推定してみます。

## 正規線形回帰モデルのベイズ推定
やることは正規分布の平均についてのベイズ推定のときと同じです。
対数尤度関数は

```math
\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}$の事前分布を

```math
\begin{align}
p(\boldsymbol{w}) &= \mathcal{N}(\boldsymbol{w};\boldsymbol{\mu}_0,\boldsymbol{\Sigma}_0)
\end{align}
```
とすれば、事後分布も正規分布：

```math
\begin{align}
p(\boldsymbol{w}|\mathcal{D}_N)
 &= \mathcal{N}(\boldsymbol{w};\boldsymbol{\mu}_N,\boldsymbol{\Sigma}_N)
\tag{1}
\end{align}
```
となります。事後分布の対数は

```math
\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}
```
です。あるいは計画行列：

```math
\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$を用いると

```math
\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)式の対数：

```math
\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}
```
と係数を対応させて

```math
\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}$の事後分布は

```math
\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}$とすると

```math
\begin{align}
\boldsymbol{\mu}_N \rightarrow (\boldsymbol{x}_N^T \boldsymbol{x}_N)^{-1} \boldsymbol{x}_N^T \boldsymbol{y}_N
\end{align}
```
となり、最小二乗解となります。

予測分布は

```math
\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}
```
です。

次回は[「ベイズ推定５：ロジスティック回帰のラプラス近似によるベイズ推定」](https://qiita.com/msekino/items/69887aa05abd05b4ef82)です。