# 最尤推定, MAP推定を用いた回帰
この記事は[古川研究室 Workout_calendar](https://qiita.com/flab5420/private/aee38a24be700d146817) 12日目の記事です。

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

# 概要

はじめに，今回のイベントで3本記事を出す予定なので簡単に内容を紹介します．第1回目（本記事）では，最尤推定とMAP推定を用いた回帰について，第2回目では，ベイズ推定を用いた回帰について，第3回目では，正規分布について勉強した記事になる予定です．第1,2回目は，確率モデルのパラメータ推定に関する内容なので，数式を追って理論的に理解を目指します．第3回目では，これら手法の中でも頻出する正規分布について，勉強したことをまとめたいと思います．

具体的に，第1回目（本記事）では，最尤推定やMAP推定を用いて回帰モデルのパラメータ$\boldsymbol{w},\beta$を推定し，新規データの予測分布を求めます．また，最尤推定やMAP推定で得られる結果が最小二乗誤差や正則化最小二乗誤差になっていることを確認します．

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

### 「記事内で行ったこと」
- [x] (理論)最尤推定を用いたモデルのパラメータと予測分布の導出
- [x] (理論)MAP推定を用いたモデルのパラメータと予測分布の導出
- [x] (実装)最尤推定，MAP推定を用いた回帰の実装

#### ここで，簡単に用語などの説明を行いたいと思います

## 最尤，MAP，ベイズ推定について
データセット$X$が与えられている状況で，モデルのパラメータ$\boldsymbol{w}$を推定する統計的な手法に，最尤推定・MAP推定・ベイズ推定があります．これらの手法は，最尤推定とMAP推定がモデルのパラメータ$\boldsymbol{w}$に関して，最も「良い」答えをただ一つ求める点推定（それ以外の他の可能性は全て切り捨てる）なのに対して，ベイズ推定はパラメータ$\boldsymbol{w}$に関して，「良い」だけでなく「まあまあ良い」答えなど，あらゆる答え（可能性）を含んだ分布を推定できることが大きな違いです．ここで，ベイズの定理は次のような式になっています．
\begin{eqnarray}
p(\boldsymbol{w}|\boldsymbol{x})=\frac{p(\boldsymbol{x}|\boldsymbol{w})p(\boldsymbol{w})}{p(\boldsymbol{x})}
\end{eqnarray}

以下は，この定理を使ったパラメータ推定方法です．

①最尤推定(事前知識なしで尤度$p(\boldsymbol{x}|\boldsymbol{w})$を最大にするパラメータを選ぶ)

$$\boldsymbol{w}_{ML}=\underset{w}{\arg \max \,}p(\boldsymbol{x}|\boldsymbol{w})$$

②MAP推定(事前知識ありで事後分布$p(\boldsymbol{w}|\boldsymbol{x})$を最大にするパラメータを選ぶ)

$$\boldsymbol{w}_{MAP}=\underset{w}{\arg \max \,}p(\boldsymbol{x}|\boldsymbol{w})p(\boldsymbol{w})$$

③ベイズ推定（パラメータ$\boldsymbol{w}$の分布を推定する）

$$p(\boldsymbol{w}|\boldsymbol{x})=\frac{p(\boldsymbol{x}|\boldsymbol{w})p(\boldsymbol{w})}{p(\boldsymbol{x})}$$

## 尤度とは
$p(A|B)$の形は，$B$が固定で$A$が変数のときは条件付き確率（事後分布），$A$が固定で$B$が変数のときは尤度と呼ばれます．条件付き確率は，先に$B$を固定した後，対象$A$（変数）がどのぐらいで起こりやすいかを確率で表しています．一方尤度は，対象$A$（固定）が分かった後，対象$A$から見て$B$（変数）だったらその対象が出現するかもしれない，という「妥当さ（尤もらしさ）」を表しています．

本記事内では，データセット$X$（固定），モデルのパラメータ$\boldsymbol{w}$（変数）の条件の下で，最も大きな尤度$p(X|\boldsymbol{w})$は$X$を一番上手く表している分布になっているに違いないと仮定して，尤度が最も高いときの$\boldsymbol{w}$をモデルのパラメータに選びます．この推定方法は最尤推定と呼ばれます．

と書かれています．本記事では，仮定したモデルのパラメータを推定するために最尤推定とMAP推定を使用します．

# 最尤推定を用いた回帰

## 課題設定
**与えられているもの**
入力$x\in\mathbb{R}$と出力$y\in\mathbb{R}$のデータセット$X=\{(x_i,y_i)\}_{i=1}^N$

**推定したいもの**
$y \simeq f(x,\boldsymbol{w},\beta)$となるパラメータ$\boldsymbol{w},\beta$

最尤推定で，対数尤度関数$L(\boldsymbol{w},\beta)$から$\boldsymbol{w}_{ML},\beta_{ML}$を推定し，新規入力$x^*$に対する写像$f(x^*,\boldsymbol{w}_{ML})$と予測分布$p(y^*|x^*,X,\boldsymbol{w}_{ML},\beta_{ML})$を求める．

## 仮定

・データセット$X=\{(x_i,y_i)\}_{i=1}^N$は独立同分布(i.i.d.)とします．

・出力$y$は写像$f$とガウスノイズ$\epsilon$の和で表されるとします．

$$y=f(x)+\epsilon$$

・写像$f$に対して，次のモデルを仮定します．

$$f(x)=\boldsymbol{w}^T\boldsymbol{\phi}(x)$$

・ここで，基底関数$\boldsymbol{\phi}(x)$とパラメータ$\boldsymbol{w}$は$D$次元で，次のようになっています．

$$
\boldsymbol{\phi}(x)=(\phi_1(x),...,\phi_{D}(x))^T\\
\boldsymbol{w}=(w_1,...,w_{D})^T
$$

## 解析解を導出する手順

① データセット$X$から尤度$\ln p(\boldsymbol{y}|\boldsymbol{x},\boldsymbol{w},\beta)$を求め，対数尤度関数$L(\boldsymbol{w},\beta)$とみなす．

② 最尤推定で，対数尤度関数$L(\boldsymbol{w},\beta)$から$\boldsymbol{w}_{ML}$を求める

③ 最尤推定で，対数尤度関数$L(\boldsymbol{w},\beta)$から$\beta_{ML}$を求める

④ $\boldsymbol{w}_{ML}$，$\beta_{ML}$によって，新規データ$x^*$の予測分布$p(y^{*}|x^{*},\boldsymbol{w}_{ML},\beta_{ML})$を求める．

## ① データセット$X$から対数尤度関数$L(\boldsymbol{w},\beta)$を求める

はじめに，観測データ$(x_i,y_i)$が与えられている状況下で，写像$f(x_i)=\boldsymbol{w}^T\boldsymbol{\phi}(x_i)$と仮定したとき，観測データ$y_i$は，$y_i=f(x_i)+\epsilon_i$となります（$\epsilon_i$は，ガウスノイズ$N(0,\beta^{-1})$を表しています）．

今，尤度$p\left(y_i|x_i,\boldsymbol{w},\beta\right)$を求めるとすると，$y_i$はガウス分布に従うため，

\begin{eqnarray}
p\left(y_i|x_i,\boldsymbol{w},\beta\right)
&=&N\left(y_i|f(x_i),\beta^{^-1}\right)\\
&=&\sqrt\frac{\beta}{{2\pi}}\exp\left(-\frac{\beta}{2}(y_i-f(x_i))^{2}\right)\\
&=&\sqrt\frac{\beta}{{2\pi}}\exp\left(-\frac{\beta}{2}\left(y_i-\boldsymbol{w}^{T}\boldsymbol{\phi}\left(x_i\right)\right)^{2}\right)
\end{eqnarray}

となります．ここで，計算を簡単にするために対数尤度にすると，

$$\ln p\left(y_i|x_i,\boldsymbol{w},\beta\right)=\frac{1}{2}\ln \frac{\beta}{2\pi}-\frac{\beta}{2}\left(y_i-\boldsymbol{w}^{T}\boldsymbol{\phi}(x_i)\right)^{2}$$

となります．よって，全体のデータセット$X$の対数尤度は，

\begin{eqnarray}
\sum_{i=1}^{N}\ln p(y_i|x_i,\boldsymbol{w},\beta)
&=&\sum_{i=1}^{N} \{ \frac{1}{2}\ln \frac{\beta}{2\pi}-\frac{\beta}{2}(y_i-\boldsymbol{w}^{T}\boldsymbol{\phi}(x_i))^{2} \} \\
&=&-\frac{\beta}{2}\sum_{i=1}^{N}(y_i-\boldsymbol{w}^{T}\boldsymbol{\phi}(x_i))^{2}+\frac{N}{2}\ln \beta-\frac{N}{2}\ln {2\pi}\\
&=&-\frac{\beta}{2}\|\boldsymbol{y}-\boldsymbol{\Phi}\boldsymbol{w}\|^{2}+\frac{N}{2}\ln \beta-\frac{N}{2}\ln {2\pi}
\end{eqnarray}

となります．ここで，$N,D$はデータ数と次元を表し，$\boldsymbol{y},\boldsymbol{w},\boldsymbol{\Phi}$はそれぞれ，
\begin{eqnarray}
\boldsymbol{y}&=&\left(y_1,...,y_N\right)^{T}\\
\boldsymbol{w}&=&(w_1,...,w_D)^{T}\\
\boldsymbol{\Phi} &=&
\left(
\begin{array}{cccc}
\phi_{1}(x_1) & \cdots & \phi_{D}(x_1)\\
\vdots & \ddots & \vdots\\
\phi_{1}(x_N) & \cdots & \phi_{D}(x_N)\\
\end{array}
\right)\\
\end{eqnarray}

となっています．ここで，求めた対数尤度は，パラメータ$\boldsymbol{w},\beta$による関数とみなせるので，対数尤度関数$L(\boldsymbol{w},\beta)$と置くと，

\begin{eqnarray}
L(\boldsymbol{w},\beta)&=&-\frac{\beta}{2}\|\boldsymbol{y}-\Phi\boldsymbol{w}\|^{2}+\frac{N}{2}\ln \beta-\frac{N}{2}\ln {2\pi}\\
&=&-\beta\frac{1}{2}\|\boldsymbol{y}-\Phi\boldsymbol{w}\|^{2}+\frac{N}{2}\ln \beta-\frac{N}{2}\ln {2\pi}\\
&=&-\beta E(\boldsymbol{w})+\frac{N}{2}\ln \beta-\frac{N}{2}\ln {2\pi}
\end{eqnarray}

となります．ここで，$E(\boldsymbol{w})=\frac{1}{2}\|\boldsymbol{y}-\Phi\boldsymbol{w}\|^{2}$です．

この式のパラメータ$\boldsymbol{w}$について注目すると，第二項と第三項は定数であり，第一項は$\beta$（スカラー）倍された目的変数とモデルの出力の二乗誤差$E(\boldsymbol{w})$になっていることが分かります．つまり，パラメータ$\boldsymbol{w}$に関して尤度関数を最大化するということは二乗誤差を最小にするという問題と等価であることを示しています．このことから，回帰における最小二乗法を解くという行為は最尤推定という統計的手法によって妥当性が裏付けられた解き方であるということが分かります．ただし，ガウスノイズの仮定の下でのみ成り立つことに注意しなければいけません．

## ② 最尤推定で，対数尤度関数$L(\boldsymbol{w},\beta)$から$\boldsymbol{w}_{ML}$を求める

次に，$\boldsymbol{w}_{ML}$を求めるため，対数尤度関数$L(\boldsymbol{w},\beta)$を$\boldsymbol{w}$で偏微分し，最大化します．

\begin{eqnarray}
\frac{\partial L(\boldsymbol{w},\beta)}{\partial \boldsymbol{w}}
&=&\frac{\partial}{\partial \boldsymbol{w}}\{-\frac{\beta}{2}\|\Phi\boldsymbol{w}-\boldsymbol{y}\|^{2}\}\\
&=&\frac{\partial}{\partial \boldsymbol{w}}\{-\frac{\beta}{2}(\Phi\boldsymbol{w}-\boldsymbol{y})^{T}(\Phi\boldsymbol{w}-\boldsymbol{y})\}\\
&=&-\frac{\beta}{2}\{\Phi^{T}(\Phi\boldsymbol{w}-\boldsymbol{y})+(\Phi\boldsymbol{w}-\boldsymbol{y})^{T}\Phi\}\\
&=&-\frac{\beta}{2}\{\Phi^{T}(\Phi\boldsymbol{w}-\boldsymbol{y})+\Phi^T(\Phi\boldsymbol{w}-\boldsymbol{y})\}\\
&=&-\beta \Phi^{T}(\Phi\boldsymbol{w}-\boldsymbol{y})
=0
\end{eqnarray}

$$途中式で，A^TB=B^TAを使用しました$$

$$※(AB)^T=B^{T}A^Tより，(A^TB)^T=B^TAを使用した．\\
但し，A^TBが対称行列であることを示さなければいけないが，\\今回の場面ではまだ示していない$$

従って，
\begin{eqnarray}
\Phi^{T}\Phi\boldsymbol{w}=\Phi^{T}\boldsymbol{y}
\end{eqnarray}

となり，最尤推定で求められた$\boldsymbol{w}_{ML}$は，

\begin{eqnarray}
\boldsymbol{w}_{ML}
&=&(\Phi^{T}\Phi)^{-1}\Phi^T\boldsymbol{y}\\
&=&\Phi^{\dagger}\boldsymbol{y}
\end{eqnarray}

となります．この式は，最小二乗法の正規方程式として知られています．また，$N×D$行列$\Phi^+$はムーア＝ペンローズの一般化逆行列もしくは疑似逆行列と呼ばれます．これは，通常の逆行列の概念への拡張とみなすことができ，もし$\Phi$が正則であれば，$(AB)^{-1}=B^{-1}A^{-1}$より，

\begin{eqnarray}
\Phi^{\dagger}
&=&(\Phi^{T}\Phi)^{-1}\Phi^T\\
&=&\Phi^{-1}(\Phi^T)^{-1}\Phi^T\\
&=&\Phi^{-1}
\end{eqnarray}

となって，$\Phi^{\dagger}=\Phi^{-1}$が成り立つことが分かります．

## ③ 最尤推定で，対数尤度関数$L(\boldsymbol{w},\beta)$から$\beta_{ML}$を求める

尤度関数$L(\boldsymbol{w},\beta)$は

\begin{eqnarray}
L(\boldsymbol{w},\beta)
=-\beta \frac{1}{2}\|\boldsymbol{y}-\boldsymbol{\Phi}\boldsymbol{w}\|^{2}+\frac{N}{2}\ln \beta-\frac{N}{2}\ln {2\pi}
\end{eqnarray}

より，$\beta$について最尤推定を行うと，

\begin{eqnarray}
\frac{\partial L\left(\boldsymbol{w},\beta\right)}{\partial \beta}
&=&-\frac{1}{2}\|\boldsymbol{y}-\boldsymbol{\Phi}\boldsymbol{w}\|^{2}+\frac{N}{2\beta}
&=&0
\end{eqnarray}

従って，$\beta_{ML}^{-1}$は，

\begin{eqnarray}
\frac{1}{\beta_{ML}}
=\frac{1}{N}\sum_{i=1}^{N}(y_i-\boldsymbol{w}^{T}\boldsymbol{\phi}(x_i))^{2}
\end{eqnarray}

となります，つまり，この式から分散（$\beta_{ML}$の逆数）は，平均二乗誤差と等価であることが分かりました．

## ④ $\boldsymbol{w}_{ML}$，$\beta_{ML}$によって，予測分布$p(y^{*}|x^{*},\boldsymbol{w}_{ML},\beta_{ML}^{-1})$を求める．

導出した$\boldsymbol{w}_{ML}$と$\beta_{ML}$を用いると，予測分布$p(y^{*}|x^{*},\boldsymbol{w}_{ML},\beta_{ML})$は，次の式になります．

\begin{eqnarray}
p(y^{*}|x^{*},\boldsymbol{w}_{ML},\beta_{ML})
&=&N(y^{*}|f(x^{*}),\beta_{ML}^{-1})\\
&=&N(y^{*}|\boldsymbol{w}_{ML}^{T}\boldsymbol{\phi}(x^{*}),\beta_{ML}^{-1})
\end{eqnarray}

# MAP推定を用いた回帰

## 課題設定
**与えられているもの**
観測データ$x\in\mathbb{R}$とラベル$y\in\mathbb{R}$のデータセット$X=\{(x_i,y_i)\}_{i=1}^N$

**推定したいもの**
$y \simeq f(x,\boldsymbol{w})$となるパラメータ$\boldsymbol{w}$

事後分布$p(\boldsymbol{w}|X,\alpha)$を最大化する$\boldsymbol{w}_{MAP}$を推定し，新規データ$x^*$に対する写像$f(x^*,\boldsymbol{w}_{MAP})$と予測分布$p(y^*|x^*,X,\boldsymbol{w}_{MAP},\alpha)$を求める．

## 仮定
・$\boldsymbol{w}$の事前$p(\boldsymbol{w}|\alpha)$は，滑らかであまり大きな$\boldsymbol{w}$でないと仮定する．ここで，$\alpha$は，ガウス分布の精度を表すハイパーパラメータです．

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

・ベイズの定理から，事後分布$p(\boldsymbol{w}|X,\alpha)$を尤度$p(X|\boldsymbol{w})$と事前分布$p(\boldsymbol{w}|\alpha)$で指定します．

・他の仮定は，最尤推定と同じです．

### 解析解を導出する手順
① 事前知識から事前分布$p(\boldsymbol{w|\alpha})$を仮定する

② 尤度$p(X|\boldsymbol{w})$と事前分布$p(\boldsymbol{w}|\alpha)$から事後分布$p(\boldsymbol{w}|X,\beta,\alpha)$を指定する

③ 事後分布$p(\boldsymbol{w}|X,\alpha)$を最大化する$\boldsymbol{w}_{MAP}$を求める

④ $\boldsymbol{w}_{MAP}$を用いて，新規データ$x^*$に対する$y^*$を予測する

### ① 事前知識から事前分布$p(\boldsymbol{w|\alpha})$を仮定する

最尤推定と同様に$y_i=\boldsymbol{w}^T\boldsymbol{\phi}(x_i)+\epsilon_i$（$\epsilon$はガウスノイズ）とすると，予測分布は，$p(y_i|x_i)=N(y_i|f(x_i),\beta^{-1})$と仮定できます．MAP推定が最尤推定と異なる点は，ベイズ的アプローチである事前知識を適用できることです．ここで，事前分布$p(\boldsymbol{w}|\alpha)$は滑らかで$\boldsymbol{w}$はあまり大きな値でないと仮定すると．次のガウス分布で定義できます．

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

ここで，$\alpha$はガウス分布の精度を表しており，ハイパーパラメータです．

### ② 尤度$p(X|\boldsymbol{w})$と事前分布$p(\boldsymbol{w}|\alpha)$から事後分布$p(\boldsymbol{w}|X,\alpha)$を指定する

ベイズの定理から，データセット$X$に対するパラメータ$\boldsymbol{w}$の事後分布$p(\boldsymbol{w}|X,\alpha)$は，

$$p(\boldsymbol{w}|X,\alpha)=\frac{p(X|\boldsymbol{w})p(\boldsymbol{w}|\alpha)}{p(X)}$$

となります．尤度$p(X,\boldsymbol{w})$と事前分布$p(\boldsymbol{w}|\alpha)$を用いて対数事後分布$p(\boldsymbol{w}|X,\alpha)$を次のように指定します．

\begin{eqnarray}
\ln p(\boldsymbol{w}|X) &=& \ln p(X|\boldsymbol{w})+\ln p(\boldsymbol{w}|\alpha) - \ln p(X)\\
&=&\ln p(X|\boldsymbol{w})+\ln p(\boldsymbol{w}|\alpha) + const
\end{eqnarray}

ここで，$p(X)$は，$p(X)=\int p(X|\boldsymbol{w})p(\boldsymbol{w})d\boldsymbol{w}$で，$\boldsymbol{w}$によらない定数項$const$として無視できる．更に，対数事後分布$\ln p(\boldsymbol{w}|X)$をパラメータ$\boldsymbol{w}$の目的関数$M(\boldsymbol{w})$とみなすと，

\begin{eqnarray}
M(\boldsymbol{w})
&=&-\frac{\beta}{2}\|\Phi\boldsymbol{w}-\boldsymbol{y}\|^{2}-\frac{\alpha}{2}\|\boldsymbol{w}\|^{2}+\frac{n}{2}\ln \beta-\frac{n}{2}\ln {2\pi} + const\\
&=&-\frac{\beta}{2}\|\Phi\boldsymbol{w}-\boldsymbol{y}\|^{2}-\frac{\alpha}{2}\|\boldsymbol{w}\|^{2} + const
\end{eqnarray}

となります．この式の，第3項と第4項はパラメータ$\boldsymbol{w}$に対して定数であるため，$const$にまとめています．また，$M(\boldsymbol{w})$を$\beta$で割って，$\lambda=\frac{\alpha}{\beta}$とするとRidge回帰の目的関数$E(\boldsymbol{w})$と一致することが分かります．

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

つまり，ノイズがガウス分布に従い，パラメータ$\boldsymbol{w}$の事前分布もガウス分布に従うとき，MAP推定とRidge回帰の結果は等価であることを示しています．

また，MAP推定は尤度だけでなく事前分布も考慮できるので，最尤推定に比べ情報量が増えていることが分かります．

### ③ 事後分布$p(\boldsymbol{w}|X,\alpha)$を最大化する$\boldsymbol{w}_{MAP}$を求める

事後分布を最大にするようなパラメータ$\boldsymbol{w}_{MAP}$は，$M(\boldsymbol{w})$を$\boldsymbol{w}$に関して偏微分することで求められます．

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

移項すると次の式のようになります．

$$(\Phi^{T}\Phi\boldsymbol+\lambda I)\boldsymbol{w}=\Phi^{T}\boldsymbol{y}$$

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

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

となります．

### ④ $\boldsymbol{w}_{MAP}$を用いて，新規データ$x^*$に対する$y^*$を予測する

推定した$\boldsymbol{w}_{MAP}$を用いた，新規データ$x^*$に対する予測$y^*$は次の式になります．

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

# スクラッチ実装
導出した式を使い，スクラッチで実装しました．

## 最尤推定を用いた回帰の実装

# 最尤推定
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(precision=3)

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] for x in X])
    return result

N = 10
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])
y = np.array([[i] for i in y])

w_ML = np.dot(np.dot(np.linalg.inv(np.dot(Phi(x).T, Phi(x))), Phi(x).T), y.ravel())
w_ML = np.array([[w] for w in w_ML])

# テスト
predict = np.dot(w_ML.T, Phi(x_test).T).ravel()
variance_ML = np.sum(np.square(y-np.dot(Phi(x), w_ML))) / N
sigma = np.sqrt(variance_ML)

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].fill_between(x_test.ravel(), predict+sigma*3, predict-sigma*3, alpha=0.1, label='3sigma')
ax[1].fill_between(x_test.ravel(), predict+sigma, predict-sigma, alpha=0.1, label='sigma')
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_ML.ravel()):
    print('w{}: {:.2f}'.format(i,w),end=' ')

### 実行結果

真の関数がsinで，10点の観測（学習）データを3次多項式でフィッテングした時の結果です（左図）．観測データ点に対して，ある程度上手く曲線を引けています．テストデータ（右図）に対しても，$x$の定義域が-1.5から1.5程度まではある程度上手くできています．また，パラメータ$\boldsymbol{w}$を確認してみると，あまり大きな値になっていないことが分かります．


同じ条件で，9次多項式を使った結果です．10個の観測データ全てを通る曲線になっています（左図）．しかし，テストデータに対して全く上手く回帰できていないことが分かります（右図）．これは，学習データに対して過学習しているのが原因と考えられます．また，過学習が原因で，最尤推定で求めた$\beta$も非常に小さい値になっていることが右図から分かります．ここで，パラメータ$\boldsymbol{w}$を確認してみると極端に大きな値になっており，それが原因で傾きがとても急な曲線になっていることが分かります．

## 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次多項式を使いMAP推定を行った結果です．$\alpha$は0.03，$\beta$は0.3（既知と仮定）で行っています．実験結果は，学習データ，テストデータともに$x$の定義域が-1から1の範囲内で，データに対し，上手く曲線がフィッテングできていることが分かります．更に，事前分布（正則化項）によって，パラメータ$\boldsymbol{w}$はあまり大きな値でないことが分かります．


同じ条件下で，$\alpha$を0にした結果です．事前分布による$\boldsymbol{w}$の制約がないためこのようになります．この結果から，最尤推定の9次多項式の結果とMAP推定で$\alpha$が0の結果は，全く同じになることが確認できます．


同じ条件下で，$\alpha$を100にした結果です．事前分布の制約が強すぎる（ガウス分布の分散が小さい）のが原因で，上手く曲線が引けていないことが分かります．結果，パラメータ$\boldsymbol{w}$は，ほとんど0に近い値となってます．

# 感想
数式の導出とスクラッチで実装したことによって，最尤推定やMAP推定を使った回帰のイメージがつかめたと思います．また，MAP推定について，モデルの設計で付きまとうバイアス・バリアンス問題をハイパーパラメータの選択問題に置き換えられることを，実際に数値を色々変えてみて確認しました．

## 参考文献
[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)