ベイズ推定入門

#0. 概要
ベイジアンニューラルネットワーク等の分野が隆盛を見せている。この技術を手足のごとく使いこなせるようになるには、ベイズ推定を理解する必要があるそうだ。今回はベイズ推定の中でも最も簡単であろう線形回帰を実装したいと思う。線形回帰をベイズ化することで嬉しい点は、推論時の確率密度関数を計算することができ、推論誤差を分散として可視化できる点にある。言い換えると不確実性を可視化できるということである。

#1. 理論
理論について説明する。

## 1.1. 基底関数
基底関数とは、ある関数を表現する基底ベクトルのことである。すなわち、対象となる空間に属する全ての関数は、この基底関数の線型結合で$\phi$で表される。

線形回帰やニューラルネットワーク等の線形結合で言い換えると、<b>重みが掛けられる側(入力)のことである。</b>

基底関数は以下のようなものがある。

多項式基底

``` math
\phi _{i}\left( x\right) =x^{i} (i=0,\ldots ,M-1)
```

ガウス基底

``` math
\phi _{i}\left( x\right) =\exp \left\{ -\dfrac {\left( x-\mu_{i} \right) ^{2}}{2s^{2}}\right\} (i=0,\ldots ,M-1)
```

シグモイド基底

``` math
\phi _{i}\left( x\right) =\sigma \left( \dfrac {x-\mu_{i} }{s}\right) (i=0,\ldots ,M-1)
```

これら基底関数を$w_i$倍することで教師$y$との誤差を最小にするパラメータベクトル$w$を推定することが目標である。

## 1.2. 線形回帰

ここで最もシンプルな線形回帰の問題を考える。

まずは
・基底関数に線形回帰を採用し、
・$w$と線形結合させ、それを関数$f(x)$とする

``` math
f(x)=w^{T}\phi(x)=\sum ^{n-1}_{i=0}W_{i}\phi _{i}(x)
```

次に、損失関数を定義する。ここで教師データを$t$とする。$n$は与えられるデータの数である。

``` math
E(w)=\dfrac {1}{2}\sum ^{N}_{n=1}\left( f\left( x_n\right) -t_{n}\right) ^{2}
```

$E(w)$を最小化することで、与えられたデータの関数の形が求まる。

上記が一般的な線形回帰である。

## 1.2. ベイズ推定
ベイズ推定では線形回帰と異なり、<b>推定結果を確率分布で出力する。</b>
すなわち、予測値の分布をプロットすることができる。

予測分布は以下の数式で表記できる。

``` math
p(t_{estimate}|x_{new}, x_{observced}, t_{observed})
```

この$p$は確率分布を意味しており、$x_{observced}, t_{observced}$は既に観測されていて新たな入力$x_{new}$が得られた際の出力$t_{estimate}$の確率分布という意味である。

ここで周辺確率密度関数が以下のように定義される場合、

``` math
P(x) = \int_{y}P(x,y)dy
```

以下のように変形できる。

```math
p(t_{estimate}|x_{new}, x_{observced}, t_{observed}) = \int p(t, w | x_{new}, x_{ovserved}, t_{observed})dw
```

ここで、$t, w | x_{new}, x_{ovserved}, t_{observed}$は同時分布のため、以下のような変形則を用い、

``` math
\begin{eqnarray*}
P(X,Y|Z) &=& \frac{P(X,Y,Z)}{P(Z)} \\
&=& \frac{P(X|Y,Z)P(Y,Z)}{P(Z)} \\
&=& \frac{P(X|Y,Z)P(Y|Z)P(Z)}{P(Z)} \\
&=& P(X|Y,Z)P(Y|Z)  
\end{eqnarray*}
```

以下のように変形できる。

```math
p(t_{estimate}|x_{new}, x_{observced}, t_{observed}) = \int p(t_{estimate} | w, x_{new}, x_{ovserved}, t_{observed})p(w | x_{new}, x_{ovserved}, t_{observed})dw
```

ここで各変数の関係を考える。まず$w, x$は新たな入力のため独立事象である。予測値$t_{estimate}$は$w$が与えられた条件下において$w$と条件付き独立である。

$w$は訓練データから決まることを念頭に置くと、$w$が決定すれば訓練データ$x, t$が何であっても予測値に影響を与えない。

そこで、不要な条件を取り除くと以下のようになる。

```math
p(t_{estimate}|x_{new}, x_{observced}, t_{observed}) = \int p(t_{estimate} | w, x_{new})p(w|x_{ovserved}, t_{observed})dw
```

$p(t_{estimate} | w, x_{new})$は予測値は入力と重みから決まるということである。
$p(w|x_{ovserved}, t_{observed})$は観測された学習データ$x, t$から$w$が求められるということ。

このように容易に理解できる形へと変形できた。式全体としては周辺化の計算といえる。

では、過去に学習したデータセットに近い距離における推定値の予測分布の分散が小さくなるのはなぜか？
それは、全ての$w$について事後確率の重み付けをしながら分布$p(t_{estimate} | w, x_{new})$を積分している。よって$w$の事後確率が高い(分布の分散が小さい)ところでは$p(t_{estimate} | w, x_{new})$も高い値を取るためである。



## 1.3. ガウス過程回帰

まず初期分布を生成する。


ガウス基底は以下である。

``` math
\phi _{i}\left( x\right) =\exp \left\{ -\dfrac {\left( x-\mu_{i} \right) ^{2}}{2s^{2}}\right\} (i=0,\ldots ,M-1)
```

この時、$\mu_i$は基底関数の位置を示す。均等に配置した場合は以下のようなグラフができる。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/054f82f1-d547-84d1-af0d-7ab8df411c2a.png)

多項式基底は以下である。

``` math
\phi _{i}\left( x\right) =x^{i} (i=0,\ldots ,M-1)
```

こちらも均等に12個生成するとこのような結果が得られる。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/b061a91d-130c-38a2-3296-6fad005a0322.png)

<b>
結局はこの一本一本を学習させていき、全ての関数で平均値を取ると推定値になる。
また、各々の分散を取るとそれが信頼度になる。
</b>

多項式の例だと、入力値の1乗の関数、2乗の関数、3乗の関数…といった具合である。

``` python
def phi(x): 
    s = 0.1
    return np.append(1, np.exp(-(x - np.arange(0, 1 + s, s)) ** 2 / (2 * s * s)))

def poly_phi(x): 
    m = 12
    return x ** np.arange(0, m)

X = np.linspace(0, 1)
Y = np.array([phi(x) for x in X])
ax.append(fig.add_subplot((max_plot_num), 1, ax_count))
ax[-1].set_xlabel('x')
ax[-1].set_ylabel('y')
ax[-1].set_title('dist')
for i in range(1, 12):
    ax[-1].plot(X, Y.T[i])
plt.legend()  
ax_count += 1

X = np.linspace(0, 1)
Y = np.array([poly_phi(x) for x in X])
ax.append(fig.add_subplot((max_plot_num), 1, ax_count))
ax[-1].set_xlabel('x')
ax[-1].set_ylabel('y')
ax[-1].set_title('dist')
for i in range(1, 12):
    ax[-1].plot(X, Y.T[i])
plt.legend()  
ax_count += 1
```

次に、このグラフを行列式にするため、計画行列を作成する。

``` math
\Phi = \left(\begin{array}{cccc}
\phi_1(x_1) & \cdots  & \phi_{M-1}(x_1) \\
\vdots & \ddots  & \vdots \\
\phi_1(x_N) & \cdots& \phi_{M-1}(x_N)
\end{array}
\right)
```

以下で各入力に対する各次数の計算結果が求められる。

``` python
Phi = np.array([poly_phi(x) for x in X])
```