| データ数 | ４ | １２ |
|---|---|---|
| 推定曲線 |  |
 |
推定曲線（赤）と信頼度（データ数=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