Help us understand the problem. What is going on with this article?

ベイズ線形回帰ってナンだ?(Bayesian linear regression)

データ数 12
推定曲線 image.png image.png

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

はじめに

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

データ数を1から24まで増やして、逐次学習した際の動画をつぎに示します。
( YouTube では、 https://youtu.be/ruD9246nYP0
動画と数式の表示は、「Google Chrome ウェブブラウザ」を推奨します。 
animated_1.0.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

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

リッジ回帰は、過学習を避ける手法である。
二乗和誤差と重み$\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$の正規分布は、訓練データの確率を最大にする(最尤推定)ように計算する。

 正規分布

Gaussian_Distribution (7).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)と呼ぶ。xが平均±$\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)で、1回の訓練データを示す。
$\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%の領域をつぎに示す。

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

BayesianRegression(alpha=10.0, fit_intercept=False)4.png

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

BayesianRegression(alpha=10.0, fit_intercept=False)12.png

プログラム

リッジ回帰(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

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

tshimizu8
電子工作、機械学習、脳、人工知能、哲学、ロボットなどをテーマに研究、Makeしています。
https://groups.google.com/forum/?hl=ja#!forum/robot-android-group-japan-akb
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした