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

データ数	４	１２
推定曲線

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

はじめに

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

データ数を１から２４まで増やして、逐次学習した際の動画をつぎに示します。
（ YouTube では、 https://youtu.be/ruD9246nYP0 ）
動画と数式の表示は、「Google Chrome ウェブブラウザ」を推奨します。　

今回の投稿では、実装に必要な数式をまとめました。
数式は「パターン認識と機械学習　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）

入力ベクトル $\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）

リッジ回帰は、過学習を避ける手法である。
二乗和誤差と重み$\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 。

\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$の正規分布は、訓練データの確率を最大にする（最尤推定）ように計算する。

　正規分布

$\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$はノイズの精度パラメータである。

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

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

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

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

　

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

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

　ベイズ線形回帰の例

推定曲線と信頼度６８％の領域をつぎに示す。

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

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

プログラム

リッジ回帰（RidgeRegression）

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$ 　

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

テストベッド

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実行

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