時刻	行動
10:00	起床
13:00	外出
17:00	読書など

ベイズ線形回帰

おとといの記事日記(8/2):初心に帰って回帰分析 - 日記帳では，線形回帰モデルのパラメータを最尤推定というアプローチで決めた．

しかしながら，このモデルでは過学習を抑えるために正則化項というパラメータを設ける必要があった．そこで，ベイズ的な手法を導入することで，過学習を抑えることができ，なおかつモデルのパラメータ数がデータ数より多くなっても問題ないようだ．

解説の前にガウス分布は以下の法則を満たすことを紹介する．

[A]

${\bf x}$の周辺分布と${\bf x}$が与えられた際の${\bf y}$の条件付き分布が

$$ p({\bf x})=\mathcal N({\bf x}|\mu, \Lambda^{-1}) $$ $$ p({\bf y}|{\bf x})=\mathcal N({\bf y}|A{\bf x}+b,L^{-1}) $$

として与えられるとき，${\bf y}$の周辺分布と${\bf y}$が与えられた際の${\bf x}$の条件付き分布は

$$ p({\bf y})=\mathcal N({\bf y}|A\mu +b,L^{-1}+A\Lambda^{-1}A^{-1}) $$ $$ p({\bf x}|{\bf y})=\mathcal N({\bf x}|\Sigma(A^{T}L({\bf y}-b)+\Lambda \mu),\Sigma) $$ となる．ただし，$\Sigma=(\Lambda+A^{T}LA)^{-1}$である．

回帰モデルの尤度関数は入力$X$,目標${\bf t}$，基底関数$\phi(\cdot)$,精度$\beta=\sqrt{\sigma}^{-1}$に対し， $$ p({\bf t}|X,{\bf w},\beta)=\prod_{n=1}^{N}\mathcal N(t_n|{\bf w}^{T}\phi({\bf x}),\beta^{-1}) $$ と与えられることから，共役事前分布*1は

$$ p({\bf w})=\mathcal N({\bf w}|m_0,S_0) $$

のように書くことができる．

上で紹介した式変形[A]を使えば，事後分布$p({\bf w}|{\bf t})$は $$ p({\bf w}|{\bf t})=\mathcal N({\bf w}|m_N,S_N) $$ $$ m_N=S_N(S_0^{-1}m_0+\beta\Phi^{T}{\bf t}),S_N^{-1}=S_0^{-1}+\beta\Phi^{T}\Phi $$ $$ \Phi= \left[ \begin{array}{cccc} \phi_0(x_{1})&\phi_{1}(x_{1})&\cdots & \phi_{M -1}(x_{1})\\ \phi_0(x_{2})&\phi_1(x_{2})&\cdots&\phi_{M -1}(x_{2})\\ \vdots&\vdots&\ddots&\vdots\\ \phi_0(x_{N})&\phi_1(x_{N})&\cdots&\phi_{M -1}(x_{N})\\ \end{array} \right] $$ となる．また，ガウス分布は最頻値と期待値が一致するので，事後分布を最大化させる重み${\bf w}_{MAP}$は単に${\bf w}_{MAP}=m_N$である．

ここで，事前分布を簡単のために $$ p({\bf w}|\alpha)=\mathcal N({\bf w}|0,\alpha^{-1}I) $$ として与える．このとき，$m_N,S_N$は $$ m_N=\beta S_N\Phi^{T}{\bf t},S_N^{-1}=\alpha I+\beta\Phi^{T}\Phi $$ となる．

ここで，事後分布の対数を取ると $$ \ln p({\bf w}|{\bf t})=-\frac{\beta}{2}\sum^{N}_{n=1}(t_n-{\bf w}^{T}\phi({\bf x})_n)^{2}-\frac{\alpha}{2}{\bf w}^{T}{\bf w}+\text{Const} $$ となり，この関数の${\bf w}$に関する最大化は二乗和誤差と正則化項の和の最小化に一致することに等価であることが分かる．

このとき，入力${\bf x}$に対する$t$の予測確率分布は $$ p(t|{\bf t},\alpha, \beta)=\int p(t|{\bf w},\beta)p({\bf w}|{\bf t},\alpha,\beta)\text{d}{\bf w} $$ を評価すればよく，条件付き分布が $$ p(t|{\bf x},{\bf w},\beta)=\mathcal N(t|m_N\phi({\bf x}),\beta^{-1}) $$ で，事後分布が上のように与えられていることから式変形[A]から，予測分布は $$ p(t|{\bf t},\alpha, \beta)=\mathcal N\left(t\bigg|m_N\phi({\bf x}),\frac{1}{\beta}+\phi({\bf x})^{T}S_N\phi({\bf x})\right) $$ として与えられる．

データから最適な$\alpha,\beta$を決める方法は，後日に紹介，実装するとして，まずはこのベイズ線形回帰モデルをPythonで実装する．

実装

import numpy as np
import matplotlib.pyplot as plt

class BayesLinearPolynomialRegression(object):
    def __init__(self, M=3, alpha=1e-2, beta=1):
        self.M = M
        self.alpha = alpha
        self.beta = beta
        
    def polynomialize(self, x):
        phi = np.array([x**i for i in range(self.M + 1)]).T
        return phi
    
    def fit(self, X, t):
        phi = self.polynomialize(X)
        self.w_s = np.linalg.inv(self.alpha * np.eye(phi.shape[1])
                                 + self.beta * phi.T.dot(phi))
        self.w_m = self.beta * self.w_s.dot(phi.T.dot(t))
        
    def predict(self, x):
        phi = self.polynomialize(x)
        t_m = phi.dot(self.w_m)
        t_v = 1/self.beta + np.sum(phi.dot(self.w_s) * phi, axis=1)
        t_s = np.sqrt(t_v)
        return t_m, t_v

polynomializeメソッドはサンプル点集合${\bf x}={x_n}$に基底関数を施し，以下のように定義される$\Phi$を生成するメソッドである*2

今回は多項式基底関数を用いたため，最右辺のような行列を生成する． $$ \Phi= \left[ \begin{array}{cccc} \phi_0(x_{1})&\phi_{1}(x_{1})&\cdots & \phi_{M -1}(x_{1})\\ \phi_0(x_{2})&\phi_1(x_{2})&\cdots&\phi_{M -1}(x_{2})\\ \vdots&\vdots&\ddots&\vdots\\ \phi_0(x_{N})&\phi_1(x_{N})&\cdots&\phi_{M -1}(x_{N})\\ \end{array} \right]= \left[ \begin{array}{cccc} x_{1}^{0}&x_{1}^{1}&\cdots & x_{1}^{M -1}\\ x_{2}^{0}&x_{2}^{1}&\cdots & x_{2}^{M -1}\\ \vdots&\vdots&\ddots&\vdots\\ x_{N}^{0}&x_{N}^{1}&\cdots & x_{N}^{M -1}\\ \end{array} \right] $$

fit及びpredictは，上記の数式を愚直に実装すればよいのだが，predictにおいて複数の点を予測するため，$\phi$に対してではなく$\Phi$に対する予測をさせるため若干の式変形を行う．その際，内積が要素積の総和で表されることを利用する．

def sin(x):
    return np.sin(2*np.pi*x)

def generate_data(func, noise_scale, N):
    x = np.random.uniform(size=N)
    t = func(x) + np.random.normal(scale=noise_scale, size=x.shape)
    return x, t

def display_predict(truefunc, X_train, t_train, X_test, t_test_m, t_test_s,text=None,name=None):
    fig = plt.figure(figsize=(8,6))
    plt.fill_between(X_test, t_test_m-t_test_s, t_test_m+t_test_s,color='r',alpha=0.25)
    plt.plot(X_test, truefunc(X_test), label='Target Curve', c='blue')
    plt.scatter(X_train, t_train, facecolors='none', edgecolors='black', label='Observated Points', alpha=0.6)
    plt.plot(X_test, t_test_m, c='r', label='Predict Curve')
    if text is not None:
        plt.title(text)
    plt.xlabel('$x$')
    plt.ylabel('$t$')
    plt.legend()
    plt.ylim((-1.5,1.5))
    if name is not None:
        plt.savefig('{}.png'.format(name))
        plt.close()
    else:
        plt.show()

補助関数を定義する．上から，目的曲線関数，データ生成関数，グラフ表示関数である．display_predictの引数はそれぞれ実曲線，訓練データX軸，訓練データt軸，テストデータX軸，テストデータt軸の平均，標準偏差．．．である．この関数を実行すると下記のような，$\pm \sigma$の範囲が赤く塗られた図が生成される．

f:id:Fgjiutx:20180805213332p:plain

X_train, t_train = generate_data(sin, 0.15,60)

model = BayesLinearPolynomialRegression(M=20, alpha=1e-3, beta=5)

model.fit(X_train, t_train)

X_test = np.linspace(0,1,200)
t_m, t_s = model.predict(X_test)

display_predict(sin, X_train, t_train, X_test, t_m, t_s)

データ数，$\alpha,\beta$を変えた際の予測の変化の様子は以下のようになる．