search
LoginSignup
10

More than 3 years have passed since last update.

posted at

updated at

Organization

MAP推定を用いた回帰

この記事は古川研究室 Workout_calendar 17日目の記事です.

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

概要

第1回目は,最尤推定を用いた回帰の記事でした.第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$はガウスノイズの精度でスカラーです).

\begin{eqnarray}
y_i&=&f(x_i;\mathbf w)+\epsilon_i
\end{eqnarray}
\begin{eqnarray}
f(x_i;\mathbf w)=\mathbf{w}^T\boldsymbol{\phi}(x_i)
\end{eqnarray}
\begin{eqnarray}
\epsilon_i&\sim& N(0,\beta^{-1})
\end{eqnarray}

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

\begin{eqnarray}
\boldsymbol{\phi}(x_i)&=&(\phi_1(x_i),...,\phi_{D}(x_i))^T
\end{eqnarray}
\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}$を推定します.

\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}$を取らないガウス分布と仮定すると,

\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$から次のようになります.

\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}

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

$$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)}$$

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

\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})$とみなすと次のようになります.

\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)$は

\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になる極値から求められます.

\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}$としています).

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

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

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

となります.

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

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

\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$のとき

最尤推定(第1回目)の実験で使用したデータセットと同じものを使って,$\alpha=0$(分散無限大),$\beta=0.3$(既知)で9次多項式のパラメータ$\mathbf{w}$をMAP推定したのがこの結果です.$\mathbf{w}$への制約が働いていないため過学習していることが,曲線が学習データ全てを通っていることや$\mathbf w$が極端に大きな値になっていることから分かります.また,この結果は最尤推定の9次多項式の結果と全く同じになっていることが分かります.

② 基底関数が9次多項式で$\alpha=0.03$のとき

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

③ 基底関数が9次多項式で$\alpha=100$のとき

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

まとめ

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

参考文献

[1] パターン認識と機械学習(上) C.M. ビショップ

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
What you can do with signing up
10