ベイズ線形回帰ってナンだ？（Bayesian linear regression）

|データ数| ４ |１２|
|:-----------------:|:------------------:|:------------------:|
|**推定曲線**|![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78212/368d2875-8de2-13eb-c53d-6fa7e07e0dd9.png)| ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78212/b2bfdbc6-41c3-e0c9-6e44-c76242dcfa06.png)|

　　推定曲線（赤）と信頼度（データ数=4,12）

# はじめに

NumPy 上で、ベイズ線形回帰 (Bayesian linear regression) の**逐次学習 (sequential learning、オンライン学習)**を動かしてみました。

データ数を１から２４まで増やして、逐次学習した際の動画をつぎに示します。
（ YouTube では、 https://youtu.be/ruD9246nYP0 ）
動画と数式の表示は、「Google Chrome ウェブブラウザ」を推奨します。　
![animated_1.0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78212/ddc5b949-dab1-44c1-35d9-0f74b22b102f.png)


今回の投稿では、実装に必要な数式をまとめました。
数式は「パターン認識と機械学習　C.M. ビショップ　丸善出版」（「**PRML**」と略記します）より引用しました。

また、線形回帰の基本であるリッジ回帰も説明します。
リッジ回帰とベイズ回帰の特徴をつぎに示します。

### リッジ回帰の特徴

- 重み$\boldsymbol w$を推定するために、 二乗和誤差と重み$\boldsymbol w$の正規化項の和を最小にする（**過学習を防ぐ**）
- シンプル、計算量が少なく高速
- **逐次学習** に使用できない

### ベイズ回帰の特徴

- 重み$\boldsymbol w$を推定するために、 $\boldsymbol w$の**正規分布**を計算する。このとき訓練データの確率を最大にする（最尤推定）。
- 推定の**信頼度**が得られる
- **逐次学習**に使用できる

### scikit-learn

Pythonのオープンソース機械学習ライブラリ **scikit-learn** に、多様な線形回帰APIが用意されています。

- sklearn.linear_model.Ridge
- sklearn.linear_model.BayesianRidge　 (逐次学習用APIである partial_fit が見当たらない）

### プログラム環境

- numpy                    1.17.5
- matplotlib               3.1.3   
- scikit-learn             0.22.1                   
- Python 3.6.9

  
# 数学的定式化（Mathematical formulation）

### 回帰（Regression）
![RidgeRegression(alpha=0.0, fit_intercept=False).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78212/0454f3ca-65fe-0761-67e9-47facad4de83.png)




入力ベクトル $\boldsymbol x$と、訓練ベクトル $\boldsymbol t$ の組が複数与えられた時、
　　$\boldsymbol t=y(\boldsymbol x)$ 
となる、連続関数$y$を推定することを**回帰**と呼ぶ。$y$を**回帰関数**と呼ぶ。


### 線形回帰（Linear regression）（PRML 3.1）
重みを$\boldsymbol w$ としたとき、$\boldsymbol w$ の線形結合で与えらた$y$を、**線形回帰関数**と呼ぶ。
**線形回帰関数を推定することは、重み$\boldsymbol w$を推定することになる。**

$$
y(\boldsymbol x, \boldsymbol w) =  \sum_{j = 0}^{M} w_{j} x_{j}= w_0 x_0 + w_1 x_1  + w_2 x_2 + \cdots + w_M x_M  = \boldsymbol w^{T} \boldsymbol x 　\tag{3.1}
$$

ただし、　$x_{0} =1$ で、$\boldsymbol w_{0}$ はオフセットを指定し、バイアスと呼ばれる。

$M = 1$ のとき単回帰、$M \geqq 2$ のとき重回帰と呼ぶ。

**scikit-learn** ライブラリでは属性 **intercept_**（切片） と **coef**（係数）_ で、$w_0$ と $w_1,..., w_M$ にアクセスできる。

### 多項式回帰（Polynomial regression）（PRML 3.1）
つぎの形をした線形回帰関数$y$を、**多項式回帰関数**と呼ぶ。

$$
y(\boldsymbol x, \boldsymbol w) = w_0 + w_1 x + w_2 x^2 + \cdots + w_M x^M \tag{1.1}
$$

基底関数$ \phi_j(x) = x^j$ とすると

$$
y(\boldsymbol x, \boldsymbol w) = \sum_{j=0}^{M} w_j \phi_j(x) = \boldsymbol w^T \phi (\boldsymbol x) \tag{3.3}
$$

### リッジ回帰（Ridge Regression）（PRML 3.1.4）
![RidgeRegression(alpha=10.0, fit_intercept=False).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78212/fac7e9b9-32af-2264-30c0-2eafffef26e8.png)


リッジ回帰は、過学習を避ける手法である。
二乗和誤差と重み$\boldsymbol w$の正規化項の和である、誤差関数(Error function)$\boldsymbol E$ を最小にする。 $\lambda$ は正規化係数である。

$$
E(\boldsymbol  w) = \frac{1}{2} \sum_{n=1}^N \left( t_n - \boldsymbol w^T \phi (x_n) \right)^2 + \frac{\lambda}{2} \boldsymbol w^T  \boldsymbol w \tag{3.27}
$$

 



誤差関数を最小にする $\boldsymbol w$ は、次の式で計算できる。ここで$\boldsymbol \Phi と \boldsymbol t $ は、それぞれ入力データと訓練データ、$I$ は単位行列。

$$
\boldsymbol w = (\lambda I + \boldsymbol \Phi ^T \boldsymbol \Phi)^{-1} \boldsymbol \Phi ^T \boldsymbol t \tag{3.28}
$$

$ (\boldsymbol \Phi ^T \boldsymbol \Phi)^{-1} \boldsymbol \Phi ^T$ は $\boldsymbol \Phi$ のムーア - ペンローズの擬似逆行列と呼ばれる。

$\lambda$ をゼロにすると、過学習ペナルティがなくなり、単純な二乗和誤差の線形回帰になる。

$\lambda$ の scikit-learn でのパラメータ名は alpha 。

```math
\boldsymbol E = || \boldsymbol y - \boldsymbol t||_2^2 + \alpha ||\boldsymbol w||_2^2
```



### ベイズ線形回帰（Bayesian Linear Regression）（PRML 3.3）

ベイズ線形回帰は重み$\boldsymbol w$を推定するために、 $\boldsymbol w$の**正規分布**（**平均**、**分散**）を計算する方法である。
$\boldsymbol w$の正規分布は、訓練データの確率を最大にする（最尤推定）ように計算する。

#### 　正規分布

![Gaussian_Distribution (7).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78212/4dc6f675-14bf-4ec8-686f-c77bb8ff83b9.png)

$\boldsymbol x$ の確率分布が、正規分布(normal distribution)、またはガウス分布( Gaussian distribution)である場合、つぎのように表記する。

$$
\mathcal{N} \bigl( \boldsymbol x \mid \mu, \sigma^{2} \bigr)　\tag{1.46}
$$
ここで、$\mu, \sigma^{2} $は、分布の平均(mean)と分散(variance)である。 $\sigma$は標準偏差(SD:standard deviation)と呼ばれる。
また、分散の逆数　$\beta = 1/\sigma^{2}$ を精度(precision)と呼ぶ。ｘが平均±$\sigma$の信頼度は、68%である。

標準正規分布 (Standard Normal Distribution) は、平均が 0 で分散が 1 の正規分布である。

#### 　ベイズ定理と最尤推定


ベイズ定理をつぎに示す。 $p(\boldsymbol w )$ を確率、$\mathcal{D}$を観測した事後に$\boldsymbol w$に関する不確実性を事後確率分布 $p(\boldsymbol{w} | \mathcal{D})$ とする。

$$
\begin{align}
p(\boldsymbol{w} | \mathcal{D})=\frac{p(\mathcal{D} | \boldsymbol{w})p(\boldsymbol{w})}{p(\mathcal{D}) } \tag{1.43}
\end{align}
$$

右辺の $p(\mathcal{D} | \boldsymbol{w})$ は、データ集合 $\mathcal{D}$ の評価であって、$\boldsymbol{w}$の関数とみなせる。これを尤度関数（Likelihood function）と呼ぶ。


最尤推定（Maximum likelihood）は、$p(\mathcal{D} | \boldsymbol{w})$ を最大にする値であり、観測されたデータ集合の確率を最大にする $\boldsymbol{w}$ の値を選ぶことに相当する。
尤度の最大化は、 $\boldsymbol{w}$ を決めるという観点からは、二乗和誤差の最小化と等価であり、ベイズ線形回帰とリッジ回帰の結果は一致する。

#### 　ベイズ線形回帰

重みを$\boldsymbol w$、確率 $p(\boldsymbol w )$ が正規分布、$p(\boldsymbol w \mid \boldsymbol t)$ を事後確率分布とし、最尤推定とするとつぎの漸化式が得られる。
ここで　$\boldsymbol m_{0}$と$\boldsymbol S_{0}$は事前分布、$\boldsymbol m_{N}$と$\boldsymbol S_{N}$は事後分布である。$\boldsymbol \Phi$ は計画行列 (design matrix)で、１回の訓練データを示す。
 $\beta$はノイズの精度パラメータである。


```math
\begin{align}
p(\boldsymbol w ) = \mathcal{N} \bigl( \boldsymbol w \mid \boldsymbol m_{0}, \boldsymbol S_{0} \bigr) \tag{3.48} \\
p(\boldsymbol w \mid \boldsymbol t) = \mathcal{N} \bigl( \boldsymbol w \mid \boldsymbol m_{N}, \boldsymbol S_{N} \bigr) \tag{3.49} \\
\boldsymbol m_{N} = \boldsymbol S_{N} \bigl( \boldsymbol S_{0}^{-1} \boldsymbol m_{0} + \beta \boldsymbol \Phi^{T} \boldsymbol t \bigr) \tag{3.50} \\
\boldsymbol S_{N}^{-1} = \boldsymbol S_{0}^{-1} + \beta \boldsymbol \Phi^{T} \boldsymbol \Phi \tag{3.51} \\
\end{align}
```

漸化式の初期値をつぎに示す。

```math
p(\boldsymbol　w | \alpha ) = N(w | 0, \alpha^{-1} I ) \tag{3.52} 　\\

\lambda = \frac{\alpha}{\beta } 
```
 　

<!-- 追加開始
```math
\boldsymbol m_{N} = \beta \boldsymbol S_{N} \boldsymbol \Phi^{T}\boldsymbol t \tag{3.53} 
```
```math
\boldsymbol S_{N}^{-1} = \boldsymbol \alpha  \boldsymbol I + \beta \boldsymbol \Phi^{T} \boldsymbol \Phi \tag{3.54} 
```
追加終了-->

回帰で推定したい重みは　$\boldsymbol w ＝ \boldsymbol m_{N}$、信頼度はつぎで与えられる。

```math
\begin{align}
\sigma^{2}_{N}(x) = \frac{1}{\beta } + \phi(x)^{T}\boldsymbol S_{N}\phi(x) \tag{3.59}
\end{align}
```

ハイパーパラメータ $\alpha$ と $\beta$ は、「エビデンス近似」で設定できる。　(3.92)(3.95)

#### 　ベイズ線形回帰の例
推定曲線と信頼度６８％の領域をつぎに示す。

「訓練データ数＝４」
データ数が少なく推定曲線と目的曲線が大きく異なっている。
データ点から離れると信頼度６８％の領域が絞れずに大きく広がる。

![BayesianRegression(alpha=10.0, fit_intercept=False)4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78212/368d2875-8de2-13eb-c53d-6fa7e07e0dd9.png)


「訓練データ数＝１２」
データ点が追加され情報が増えたため、全体的には推定曲線が目的曲線に近づく。信頼度68％の領域も絞りこまれている。
データ点が少ない右端では、曲線の相違が大きく、信頼度68％の領域も広がったままである。

![BayesianRegression(alpha=10.0, fit_intercept=False)12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78212/b2bfdbc6-41c3-e0c9-6e44-c76242dcfa06.png)




# プログラム
### リッジ回帰（RidgeRegression）
```python:
class RidgeRegression(): 
    def __init__(self, alpha=1.0, fit_intercept=True):
        self.lambda_, self.fit_intercept = alpha, fit_intercept
        self.w, self.intercept_, self.coef_ = None, None, None
    def fit(self, X, T):
        if (self.fit_intercept): X = np.insert(X, 0, 1., X.ndim-1)
        Phi = X
        PhiT, lambdaI = Phi.T, self.lambda_ * np.identity(Phi.shape[1]) 
        PhiT_Phi = PhiT@Phi
        self.w = (np.linalg.inv(lambdaI + PhiT_Phi)@PhiT)@T #(3.28)
        if (self.fit_intercept): self.intercept_, self.coef_ = self.w[0], self.w[1:]         
        else: self.intercept_, self.coef_ = 0, self.w   
    def predict(self, X):  
        if self.fit_intercept: X =  np.insert(X, 0, 1, X.ndim-1) 
        return X @ self.w 
```

### ベイズ線形回帰（Bayesian Linear Regression）
パラメータ sigma は、ノイズの標準偏差 $\sigma$  　 

```python:
class BayesianRegression(): 
    def __init__(self, alpha=1.0, sigma=0.2, fit_intercept=True):
        self.lambda_ = alpha
        self.sigma = sigma # noise standard deviation
        self.fit_intercept = fit_intercept
        self.M = 0 # len(w) 
        self.w, self.intercept_, self.coef_,  = None, None, None
        self.m, self.S = None, None # mean, variance 
        self.beta = 1.0 / (self.sigma ** 2) # precision
        self.alpha = self.lambda_ * self.beta 
        if self.alpha == 0.: print("set alpha=0.001 (Prevent division by zero)")
    def fit(self, X, T):
        self.partial_fit(X, T)
    def partial_fit(self, X, T): 
        if (self.fit_intercept): X = np.insert(X, 0, 1., X.ndim-1)
        if self.M==0 and len(X)>0: # Initialize 
            self.M = len(X[0])
            self.m = np.zeros(self.M)
            self.S = np.identity(self.M) / self.alpha # (3.52)
        for i in range(len(X)):
            Phi = X
            phi, t = Phi[i], T[i]
            phi_T = phi.reshape(-1, 1)
            S_inv = np.linalg.inv(self.S)
            self.S = np.linalg.inv(S_inv + self.beta * phi_T * phi) # (3.51)
            self.m = self.S @ (S_inv@self.m + self.beta *  phi * t) # (3.50)
        self.w = self.m
        if (self.fit_intercept): self.intercept_, self.coef_ = self.w[0], self.w[1:]         
        else: self.intercept_, self.coef_ = 0, self.w   
    def predict(self, X, return_std=False):
        if self.fit_intercept: X =  np.insert(X, 0, 1, X.ndim-1) 
        if not return_std: return X @ self.w
        else:
            faiT, fai = X, X.T
            sigma = np.sqrt(1.0/self.beta + np.diag(faiT @ self.S @ fai)) # (3.59)
            return X @ self.w , sigma
```
### テストベッド
```python:
import numpy as np
import matplotlib.pyplot as plt
from sklearn import linear_model
from sklearn import gaussian_process

class test():
    def __init__(self, constructor, x_n=12, M=10, std=False, partial=False, graph=[]):
        def t_fn(x): return np.sin(x) # target function
        def b_fn(x,i): return x**i  # polynomial basic function
        def polynomialize(x, M): return np.array([b_fn(x,i) for i in range(M)]).T
        def graph_out(std, x, y, model, id=''):
            g_n = 100
            g_x = np.linspace(x_min, x_max, g_n)
            g_X =polynomialize(g_x, M)
            if std: g_y, ys = model.predict(g_X, return_std=True)
            else:   g_y = model.predict(g_X)
            plt.ylim(y_min, y_max)
            plt.scatter(x, y, color='steelblue', label='Training data #'+ str(len(x)))
            plt.plot(g_x, t_fn(g_x), color='steelblue', label='target')
            plt.plot(g_x, g_y, color='r', label=constructor)
            if std: plt.fill_between(g_x, g_y-ys, g_y+ys, color="pink", alpha=0.5, label="std(68%)")
            plt.xlabel('x'); plt.ylabel('y'); plt.legend(loc="lower left"); plt.grid()
            plt.savefig(constructor + id + '.png'); plt.show()
        x_min, x_max = 0, 2*np.pi
        y_min, y_max = -2.5, 2.0
        SN = 5 # S/N (signal-noise ratio)
        x = np.concatenate([np.array([ x_min, x_max]), np.random.uniform(x_min, x_max, x_n - 2)])
        X = polynomialize(x, M)
        T = t_fn(x) + np.random.randn(x_n)/SN
        model = self.model = eval(constructor)
        if not partial:
            model.fit(X, T)
            graph_out(std, x, T, model)
        else:
            for i in range(len(T)):
                model.partial_fit(X[i:i+1],T[i:i+1])
                if i+1 in graph: graph_out(std, x[:i+1], T[:i+1], model, id=str(i+1))
```

### Test実行
```python:
def case(seed):
    def init(): np.random.seed(seed)
    init(); test('RidgeRegression(alpha=0.0, fit_intercept=False)')
    init(); test('RidgeRegression(alpha=10.0, fit_intercept=False)')
    init(); test('BayesianRegression(alpha=10.0, fit_intercept=False)', partial=True, graph=[4,12], std=True)  
    ''' 
    init(); test('BayesianRegression(alpha=10.0, fit_intercept=True)', partial=True, graph=[4,12], std=True)  
    init(); test('BayesianRegression(alpha=10.0, fit_intercept=False)', partial=True, graph=range(1,13), std=True)  
    # sklearn library
    init(); test('linear_model.LinearRegression()') 
    init(); test('linear_model.Ridge(alpha=0.0)') 
    init(); print("RidgeCV.alpha =", test('linear_model.RidgeCV()').model.alpha_)
    init(); test('linear_model.BayesianRidge()', std=True)
    init(); test('linear_model.SGDRegressor(alpha=1000.)', x_n=12, M = 3) # Stochastic Gradient Descent
    init(); test('linear_model.ARDRegression(fit_intercept=False)')
    init(); test('gaussian_process.GaussianProcessRegressor(alpha=10.0)', std=True)     

    ''' 
case(1)
! ls -lt
```
# 参考資料

「Pattern Recognition and Machine Learning, Christopher Bishop, 2011/4」、「パターン認識と機械学習, 2012/4」
https://www.microsoft.com/en-us/research/uploads/prod/2006/01/Bishop-Pattern-Recognition-and-Machine-Learning-2006.pdf

sklearn.linear_model.BayesianRidge
https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.BayesianRidge.html#sklearn.linear_model.BayesianRidge

Curve Fitting with Bayesian Ridge Regression
https://scikit-learn.org/stable/auto_examples/linear_model/plot_bayesian_ridge_curvefit.html#sphx-glr-auto-examples-linear-model-plot-bayesian-ridge-curvefit-py

Scikit-learn 0.22.1 (stable) documentation (PDF 48.5 MB)
https://scikit-learn.org/stable//_downloads/scikit-learn-docs.pdf

初心に帰ってベイズ線形多項式回帰(2):エビデンス近似編::8/15　
（ハイパーパラメータ $\alpha、\beta$ のエビデンス近似）
https://fgjiutx.hatenablog.com/entry/2018/08/16/012401

https://ja.wikipedia.org/wiki/数学記号の表#代数学の記号