MAP推定を用いた回帰

この記事は[古川研究室 Workout_calendar](https://qiita.com/flab5420/private/aee38a24be700d146817) 17日目の記事です．

本記事は古川研究室の学生が学習の一環として書いたものです．内容が曖昧であったり表現が多少異なったりする場合があります．

# 概要

[第1回目](https://qiita.com/tsnry7913/items/e56e08aa23d141175a04)は，最尤推定を用いた回帰の記事でした．第2回目である本記事では，MAP推定を用いた回帰についてまとめたいと思います．

この記事で行ったことは以下になります．

- (理論)MAP推定を用いたモデルのパラメータ$\mathbf{w}$の導出
- (実装)MAP推定を用いた回帰のスクラッチ実装

## MAP推定による回帰のタスク設定
**与えられているもの**

入力$\mathbf x=(x_1,...,x_N), x\in\mathbf R$と出力$\mathbf y=(y_1,...,y_N), y\in\mathbf R$のデータセット

**仮定**

・$\mathbf x$と$\mathbf y$のデータセットは，独立同分布(i.i.d.)に従うとします．

・$y_i,f,\epsilon$を次のように仮定します（$\beta$はガウスノイズの精度でスカラーです）．

```math
\begin{eqnarray}
y_i&=&f(x_i;\mathbf w)+\epsilon_i
\end{eqnarray}
```

```math
\begin{eqnarray}
f(x_i;\mathbf w)=\mathbf{w}^T\boldsymbol{\phi}(x_i)
\end{eqnarray}
```

```math
\begin{eqnarray}
\epsilon_i&\sim& N(0,\beta^{-1})
\end{eqnarray}
```

$\;\;\;$ここで，基底関数$\boldsymbol{\phi}(x)$とパラメータ$\mathbf{w}$は$D$個あるので次のようになります．

```math
\begin{eqnarray}
\boldsymbol{\phi}(x_i)&=&(\phi_1(x_i),...,\phi_{D}(x_i))^T
\end{eqnarray}
```

```math
\begin{eqnarray}
\mathbf{w}&=&(w_1,...,w_{D})^T
\end{eqnarray}
```

・事前分布は，あまり大きな$\mathbf w$を取らないと仮定します（ハイパーパラメータ$\alpha$はガウス分布の精度でスカラーです）．

$$p(\mathbf　w|\alpha)=N(\mathbf　w|\mathbf 0,\alpha^{-1}\mathbf I)$$

**MAP推定のタスク**

事後分布が最大となるパラメータ$\mathbf w_{MAP}$を推定します．

```math
\begin{eqnarray}
\mathbf {w}_{MAP}
&=&\underset{\mathbf{w}}{\arg \max\;}p(\mathbf y|\mathbf x,\mathbf　w,\beta)p(\mathbf w|\alpha)\\\
\end{eqnarray}
```

推定した$\mathbf w_{MAP}$から，新規入力$x^*$に対する出力$y^\*$を予測できます．

# 解析解の導出
1. 事後分布を最大化するパラメータ$\mathbf{w}_{MAP}$を求める

    1-1 尤度と事前分布から変形して目的関数を求める
    1-2 目的関数を$\mathbf{w}$で偏微分して，極値を求める

2. $\mathbf{w}_{MAP}$を用いて，新規入力$x^*$に対する出力$y^*$を予測する

## 1. 事後分布を最大化するパラメータを求める
### 1-1 尤度と事前分布から変形して目的関数を求める

MAP推定が最尤推定と異なる点は，ベイズ的アプローチである事前分布を使用できることです．ここで，事前分布$p(\mathbf{w}|\alpha)$はあまり大きな$\mathbf{w}$を取らないガウス分布と仮定すると，

```math
\begin{eqnarray}
p(\mathbf{w}|\alpha)
&=&N(\mathbf{w}|\mathbf 0,\alpha^{-1}\mathbf I)\\
&=&\left(\frac{\alpha}{2\pi}\right)^{\frac{D}{2}}\exp\left(-\frac{\alpha}{2}\mathbf w^T \mathbf w\right)
\end{eqnarray}
```

となります．ハイパーパラメータ$\alpha$はガウス分布の精度を表すスカラーです．また，尤度は出力$y_i=f(x_i;\mathbf w)+\epsilon_i$から次のようになります．

```math
\begin{eqnarray}
p(\mathbf y|\mathbf x,\mathbf w,\beta)
&=&\prod_{i=1}^{N} \sqrt \frac{\beta}{2\pi}\exp\left(-\frac{\beta}{2}\left(y_i-\mathbf w^T\boldsymbol \phi(x_i)\right)^2\right)
\end{eqnarray}
```

これら尤度と事前分布から，データセットに対する事後分布はベイズの定理より

```math
$$p(\mathbf w|\mathbf x,\mathbf y,\alpha,\beta)=\frac{p(\mathbf y|\mathbf x,\mathbf w,\beta)p(\mathbf w|\alpha)}{p(\mathbf y|\mathbf x, \alpha,\beta)}$$
```

となって，更に事後分布に対数を取ると次のようになります．

```math
\begin{eqnarray}
\ln p\left(\mathbf w|\mathbf x,\mathbf y, \alpha, \beta\right)
&=& \ln p\left(\mathbf y|\mathbf x, \mathbf w, \beta\right)+\ln p\left(\mathbf w|\alpha\right) - \ln p\left(\mathbf y|\mathbf x,\alpha,\beta\right)\\
&=&\ln \left(\prod_{i=1}^{N} N(y_i|\mathbf w^T\boldsymbol \phi(x_i),\beta^{-1})\right)+\ln N(\mathbf{w}|\mathbf 0,\alpha^{-1}\mathbf I) + const\\
&=&\ln \left(\prod_{i=1}^{N} \sqrt \frac{\beta}{2\pi}\exp\left(-\frac{\beta}{2}\left(y_i-\mathbf w^T\boldsymbol \phi(x_i)\right)^2\right)\right)+\ln \left(\left(\frac{\alpha}{2\pi}\right)^{\frac{D}{2}}\exp\left(-\frac{\alpha}{2}\mathbf w^T \mathbf w\right)\right) + const\\
&=&\sum_{i=1}^N\left(\frac{1}{2}\ln \frac{\beta}{2\pi}-\frac{\beta}{2}\left(y_i-\mathbf w^T\boldsymbol \phi(x_i)\right)^2\right)+\left(\frac{D}{2}\ln\frac{\alpha}{2\pi}-\frac{\alpha}{2}\mathbf w^T\mathbf w\right)+const\\
&=&-\frac{\beta}{2}\sum_{i=1}^{N}\left(y_i-\mathbf w^T\boldsymbol \phi(x_i)\right)^2-\frac{\alpha}{2}\mathbf w^T\mathbf w+\frac{N}{2}\ln \frac{\beta}{2\pi}+\frac{D}{2}\ln \frac{\alpha}{2\pi}+const\\
&=&-\frac{\beta}{2}\sum_{i=1}^{N}\left(y_i-\mathbf w^T\boldsymbol \phi(x_i)\right)^2-\frac{\alpha}{2}\mathbf w^T\mathbf w+const\\
&=&-\frac{\beta}{2}\|\mathbf{y}-\Phi\mathbf{w}\|^{2}-\frac{\alpha}{2}\|\mathbf{w}\|^{2} + const
\end{eqnarray}
```

ここで，1行目の右辺の周辺尤度$p(\mathbf y|\mathbf x,\alpha,\beta)$と5行目の第3項と第4項は，パラメータ$\mathbf{w}$に対して定数であるため$const$にまとめています．また，対数事後分布$\ln p(\mathbf{w}|\mathbf x,\mathbf y, \alpha, \beta)$をパラメータ$\mathbf{w}$の目的関数$M(\mathbf{w})$とみなすと次のようになります．

```math
\begin{eqnarray}
M(\mathbf{w})
&=&-\frac{\beta}{2}\|\mathbf{y}-\Phi\mathbf{w}\|^{2}-\frac{\alpha}{2}\|\mathbf{w}\|^{2} + const
\end{eqnarray}
```

ここで，Ridge回帰（制約付き二乗誤差）の目的関数$E(\mathbf w)$は

```math
\begin{eqnarray}
E(\mathbf{w})=\frac{1}{2}\|\mathbf{y}-\Phi\mathbf{w}\|^{2}+\frac{\lambda}{2}\|\mathbf{w}\|^{2}
\end{eqnarray}
```

となり，$M(\mathbf{w})$を$\beta$で割って$\lambda=\frac{\alpha}{\beta}$とすると，$M(\mathbf{w})$の最大化と$E(\mathbf{w})$の最小化は同じになります．従って，ノイズ$\epsilon$と事前分布$p(\mathbf w|\alpha)$がガウス分布のとき，MAP推定とRidge回帰（制約付き二乗誤差）の結果は等価になることが分かります．

### 1-2 目的関数を偏微分して，極値を求める

事後分布が最大になるパラメータ$\mathbf{w}_{MAP}$は，目的関数$M(\mathbf{w})$が2次形式で上に凸形状なので，$\mathbf{w}$で偏微分して0になる極値から求められます．

```math
\begin{eqnarray}
\frac{\partial M(\mathbf{w})}{\partial \mathbf{w}}
&=&\frac{\partial}{\partial \mathbf{w}}\{-\frac{\beta}{2}(\mathbf{y}-\Phi\mathbf w)^{T}(\mathbf{y}-\Phi\mathbf{w})-\frac{\alpha}{2}\mathbf w^T\mathbf w+const\}\\
&=&\frac{\partial}{\partial \mathbf{w}}\{-\frac{\beta}{2}(\mathbf{y}^T\mathbf{y}-\mathbf{y}^T\Phi\mathbf{w}-\mathbf{w}^T\Phi^T\mathbf{y}+\mathbf{w}^T\Phi^T\Phi\mathbf{w})-\frac{\alpha}{2}\mathbf w^T\mathbf w+const\}\\
&=&\frac{\partial}{\partial \mathbf{w}}\{-\frac{\beta}{2}(\mathbf{y}^T\mathbf{y}-2\mathbf{w}^T\Phi^T\mathbf{y}+\mathbf{w}^T\Phi^T\Phi\mathbf{w})-\frac{\alpha}{2}\mathbf w^T\mathbf w+const\}\\
&=&-\frac{\beta}{2}(0-2\Phi^T\mathbf{y}+2\Phi^T\Phi\mathbf{w})-\alpha\mathbf{w}+0\}\\
&=&\beta\Phi^T\mathbf{y}-\beta\Phi^T\Phi\mathbf{w}-\alpha\mathbf{w}\\
&=&0
\end{eqnarray}
```

次に，両辺を$\beta$で割って式変形をすると次の式になります（$\lambda=\frac{\alpha}{\beta}$としています）．

```math
\begin{eqnarray}
(\Phi^{T}\Phi\mathbf+\lambda I)\mathbf{w}&=&\Phi^{T}\mathbf{y}\\
\end{eqnarray}
```

従って，MAP推定で求められた$\mathbf{w}_{MAP}$は，

```math
\begin{eqnarray}
\mathbf{w}_{MAP}
=(\Phi^{T}\Phi+\lambda I)^{-1}\Phi^T\mathbf{y}
\end{eqnarray}
```

となります．

## 2. 推定したパラメータから，新規入力に対する出力を予測する

推定した$\mathbf{w}_{MAP}$を用いると新規入力$x^\*$に対する出力$y^\*$は次の式になります．

```math
\begin{eqnarray}
y^{*}\simeq f(x^{*})=\mathbf{w}_{MAP}^{T}\boldsymbol{\phi}(x^*)
\end{eqnarray}
```

# MAP推定を用いた回帰の実装
導出した式を使い，スクラッチで実装しました．

## MAP推定を用いた回帰の実装
```
# MAP推定
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

def Phi(X):
    # X      : [[x0],[x1],...]
    # phi    : e.g. [1, x, x**2]
    # return : [[phi_1(X),...,phi_n(X)]]
    result = np.array([[1, x[0], x[0]**2, x[0]**3, x[0]**4, x[0]**5, x[0]**6, x[0]**7, x[0]**8, x[0]**9] for x in X])
    return result

N = 10
alpha = 0.003 # ハイパーパラメータ
beta = 0.3   # 既知と仮定しています
np.random.seed(1111)
x = np.linspace(-1,1,N)
t = np.sin(np.pi*np.linspace(-1.5,1.5,100))
y = np.sin(np.pi*x) + np.random.normal(0,beta,N)
x_test = np.linspace(-1.5,1.5,100)
y_test = np.sin(np.pi*x_test) + np.random.normal(0,beta,len(x_test))

# 縦ベクトルに変換
x = np.array([[i] for i in x])
x_test = np.array([[i] for i in x_test])

# MAP推定はここが変化するだけ
w_MAP = np.dot(np.dot(np.linalg.inv(np.dot(Phi(x).T, Phi(x))+(alpha/beta)*np.eye(10)), Phi(x).T), y)
w_MAP = np.array([[w] for w in w_MAP])

predict = np.dot(w_MAP.T, Phi(x_test).T).ravel()

fig, ax = plt.subplots(1,2,figsize=(12,4))
ax[0].plot(np.linspace(-1.5,1.5,len(t)),t,label='sin',color='b')
ax[0].plot(x_test.ravel(), predict.ravel(),label='predict curve',color='r')
ax[0].scatter(x.ravel(), y, label='train_data',color='b')
ax[1].plot(x_test.ravel(), predict.ravel(),label='predict curve',color='r')
ax[1].scatter(x_test.ravel()[::5],y_test[::5], label='test_data',color='0.1')
ax[0].set_title('train')
ax[1].set_title('test')
ax[0].set_xlim(-1.5,1.5)
ax[0].set_ylim(-1.5,1.5)
ax[1].set_xlim(-1.5,1.5)
ax[1].set_ylim(-1.5,1.5)
ax[0].legend(loc='upper right',fontsize=9)
ax[1].legend(loc='upper right',fontsize=9)
plt.show()

for i, w in enumerate(w_MAP.ravel()):
    print('w{}: {:.3f}'.format(i,w),end=' ')
```

### 実験結果

**① 基底関数が9次多項式で$\alpha=0$のとき**
![](https://i.imgur.com/DMhDa87.png)

最尤推定（[第1回目](https://qiita.com/tsnry7913/items/e56e08aa23d141175a04#%E5%AE%9F%E9%A8%93%E7%B5%90%E6%9E%9C3%E6%AC%A1%E5%A4%9A%E9%A0%85%E5%BC%8F%E3%81%AE%E3%81%A8%E3%81%8D)）の実験で使用したデータセットと同じものを使って，$\alpha=0$（分散無限大），$\beta=0.3$（既知）で9次多項式のパラメータ$\mathbf{w}$をMAP推定したのがこの結果です．$\mathbf{w}$への制約が働いていないため過学習していることが，曲線が学習データ全てを通っていることや$\mathbf w$が極端に大きな値になっていることから分かります．また，この結果は最尤推定の9次多項式の[結果](https://qiita.com/tsnry7913/items/e56e08aa23d141175a04#%E5%AE%9F%E9%A8%93%E7%B5%90%E6%9E%9C9%E6%AC%A1%E5%A4%9A%E9%A0%85%E5%BC%8F%E3%81%AE%E3%81%A8%E3%81%8D)と全く同じになっていることが分かります．

**② 基底関数が9次多項式で$\alpha=0.03$のとき**
![](https://i.imgur.com/VyCnPJ2.png)

同じ条件で，$\alpha=0.03$に変更したのがこの結果です．学習データ，テストデータともに$x$が-1から1の範囲内でデータに上手く曲線フィッテングできることが分かります（一方で，$x$が-1から1の範囲外については学習していないので上手くできていません）．更に，事前分布によって表現力が制約され，パラメータ$\mathbf w$があまり大きな値を取れていないことが分かります．

**③ 基底関数が9次多項式で$\alpha=100$のとき**
![](https://i.imgur.com/dfKnvw6.png)

同じ条件で，$\alpha=100$（分散小さめ）にしたのがこの結果です．事前分布の制約が強すぎるため，上手く学習データを通るような曲線が引けていないことが分かります．また，パラメータ$\mathbf{w}$の値を見ると，全て0付近の値になっていることが分かります．

# まとめ
数式の導出から，MAP推定を使った回帰がRidge回帰と等価になることを示しました．更に，実験から事前分布（正則化）のハイパーパラメータ$\alpha$を調整することでモデルの表現能力を制約し，モデルを変えずに適切な回帰ができることを確認しました．

## 参考文献
[1] [パターン認識と機械学習（上）　C.M. ビショップ](https://www.amazon.co.jp/%E3%83%91%E3%82%BF%E3%83%BC%E3%83%B3%E8%AA%8D%E8%AD%98%E3%81%A8%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92-%E4%B8%8A-C-M-%E3%83%93%E3%82%B7%E3%83%A7%E3%83%83%E3%83%97/dp/4621061224)