Pythonで実装 PRML 第3章 ベイズ線形回帰

パターン認識と機械学習、第3章からベイズ線形回帰が予測分布を収束させていく様子を再現します。
やっていることは、1章でやったベイズ推定とほぼ一緒で、基底関数がガウス関数(3.4)になっているだけです。
ですので、コードもほぼほぼ流用してます。

1章と比較すると、3章では数式の表現がベクトルを扱うように書きなおされてますが、この辺何か意図あるのでしょうか。
(3.16)の特徴ベクトルの概念をここで紹介したからこの書き方なのかと思ってますが、どうなんでしょう。


##実装の大まかな流れ
①求めたいのは予測分布(3.57)。

```math
 p(t|{\bf x}, {\bf t}, \alpha, \beta) = N (t| {\bf m}^T_N \phi({\bf x}),  \sigma^2_N ({\bf x})) (3.58)
```

②基底関数は(3.4)の$\phi_i(x) = \exp\{- \frac{(x -\mu_j)^2}{2s^2}\}$。
基底関数は9個と指定されているがで、配置の仕方に指示がないので均等に分布するように配置。


③予測分布の平均、分散を求める式は下記の3つなので、これをそれぞれ実装して予測分布を求める。

```math
{\bf m}_N = \beta{\bf S}_N\Phi(x_n)^{\bf T}{\bf t}(3.53)
```

```math
\sigma^2_N(x) = \beta^{-1} + \phi({\bf x})^{\bf T} {\bf S}_N \phi({\bf x}). (3.59)
```

```math
{\bf S}^{-1}_N = \alpha {\bf I} + \beta\Phi^{\bf T}\Phi(3.54)
```


#コード
図を描くにあたって、サンプルデータが増えるに従って収束する様子を@naoya_tさんの[こちらの記事(ベイズ線形回帰（PRML§3.3）の図版再現)](http://qiita.com/naoya_t/items/80ea108cebc694f5cd63)以上に美しく描くアイデアを出す自信が無かったので、ほぼ写経な感じで参考にさせて頂きました。ありがとうございました。

```python
import numpy as np
from numpy.linalg import inv
import pandas as pd
from pylab import *
import matplotlib.pyplot as plt

def func(x_train, y_train):
    # (3.4) Gaussian basis function
    def phi(s, mu_i, M, x):
        return np.array([exp(-(x - mu)**2 / (2 * s**2)) for mu in mu_i]).reshape((M, 1))

    #(3.53)' ((1.70)) Mean of predictive distribution
    def m(x, x_train, y_train, S):
        sum = np.array(zeros((M, 1)))
        for n in xrange(len(x_train)):
            sum += np.dot(phi(s, mu_i, M, x_train[n]), y_train[n])
        return Beta * phi(s, mu_i, M, x).T.dot(S).dot(sum)
    
    #(3.59)'((1.71)) Variance of predictive distribution   
    def s2(x, S):
        return 1.0/Beta + phi(s, mu_i, M, x).T.dot(S).dot(phi(s, mu_i, M, x))

    #(3.53)' ((1.72))
    def S(x_train, y_train):
        I = np.identity(M)
        Sigma = np.zeros((M, M))
        for n in range(len(x_train)):
            Sigma += np.dot(phi(s, mu_i, M, x_train[n]), phi(s, mu_i, M, x_train[n]).T)
        S_inv = alpha*I + Beta*Sigma
        S = inv(S_inv)
        return S
    
    #params for prior probability
    alpha = 0.1
    Beta = 9
    s = 0.1

    #use 9 gaussian basis functions
    M = 9

    # Assign basis functions
    mu_i = np.linspace(0, 1, M)
    
    S = S(x_train, y_train)

    #Sine curve
    x_real = np.arange(0, 1, 0.01)
    y_real = np.sin(2*np.pi*x_real)
    
    #Seek predictive distribution corespponding to entire x
    mean = [m(x, x_train, y_train, S)[0,0] for x in x_real]
    variance = [s2(x, S)[0,0] for x in x_real]
    SD = np.sqrt(variance)
    upper = mean + SD
    lower = mean - SD
    
    plot(x_train, y_train, 'bo')
    plot(x_real, y_real, 'g-')
    plot(x_real, mean, 'r-')
    fill_between(x_real, upper, lower, color='pink')
    xlim(0.0, 1.0)
    ylim(-2, 2)
    title("Figure 3.8")
    show()

if __name__ == "__main__":
    # Maximum observed data points (underlining function plus some noise)
    x_train = np.linspace(0, 1, 26)

    #Set "small level of random noise having a Gaussian distribution"
    loc = 0
    scale = 0.3
    y_train =  np.sin(2*np.pi*x_train) + np.random.normal(loc,scale,26)


    #Sample data pick up
    def randidx(n, k):
        r = range(n)
        shuffle(r)
        return sort(r[0:k])

    for k in (1, 2, 5, 26):
        indices = randidx(size(x_train), k)
        func(x_train[indices], y_train[indices])
```

##結果
![Screen Shot 2015-09-21 at 00.20.33.png](https://qiita-image-store.s3.amazonaws.com/0/62779/0ef3b525-bef2-42e2-dee2-8c941e8495f5.png)
![Screen Shot 2015-09-21 at 00.21.00.png](https://qiita-image-store.s3.amazonaws.com/0/62779/47f2ce6a-e0f1-92c4-c2c0-e31744353cc2.png)
![Screen Shot 2015-09-21 at 00.21.21.png](https://qiita-image-store.s3.amazonaws.com/0/62779/53ebc4bf-22d4-93b4-8020-223b4dca58a1.png)
![Screen Shot 2015-09-21 at 00.21.32.png](https://qiita-image-store.s3.amazonaws.com/0/62779/a82e5c8f-79c6-dac1-0809-80f82e470d40.png)

##参考
[ベイズ線形回帰（PRML§3.3）の図版再現](http://qiita.com/naoya_t/items/80ea108cebc694f5cd63)