PRML第3章のベイジアン線形回帰をPythonで実装

この記事では、PRML第3章で述べられている、ベイジアン線型回帰とエビデンス近似を実装します。

対応するjupyter notebookは筆者のgithubリポジトリにあります。PRMLの他のアルゴリズムについては、連載(?)のまとめページから飛べるようになる予定です。

1 理論のおさらい

まず、記号を定義し、PRMLに書かれている内容(数式)をざっくり復習していきましょう。
PRMLに書かれている内容を理解している方は飛ばして2節に行ってしまって大丈夫です。

1.1 設定

• $N \in \mathbb{N}$を訓練データの個数、
• $d \in \mathbb{N}$を入力の次元、
• 訓練データのうち入力を$x_0, x_1, \dots , x_{N-1} \in \mathbb{R}^d$、目的値を$t_0, t_1, \dots, t_{N-1} \in \mathbb{R}$、また$\boldsymbol{t} = {}^t (t_0, t_1 \dots, t_{N-1}) \in \mathbb{R}^N$とします。

1.2 モデルと基底関数

1.2.1 モデル

1つの$(x,t)$の組に対して, 次の形のモデルを仮定します:
$$
p(t|x,w,\beta) = \mathcal{N}(t|w^{T}\phi(x),\beta^{-1})
= \sqrt{\frac{\beta}{2\pi}} \exp\left[-\frac{\beta}{2} \left[ t - w^{T}\phi(x) \right]^2 \right],
$$
ただし、

• $\beta > 0$はGaussianの精度(hyper parameterの1つ),
• $w \in \mathbb{R}^d$はパラメーター,
• $\phi : \mathbb{R}^d \rightarrow \mathbb{R}^M $は基底関数で、具体的な形は後で定めます。$\phi(x) := (\phi_0(x), \phi_1(x), \dots, \phi_{M-1}(x))^T$と表し、$M$は基底関数の次元とします。

また、入力$x_0, x_1, \dots, x_{N-1}$が与えられた下で、出力$t_0, t_1 \dots, t_{N-1}$はそれぞれ上記の分布に従って独立に生成されると仮定します(いわゆる独立同分布(i.i.d)の仮定ですね。)。

後のために、計画行列$\Phi$ ($N \times M$行列)を定義しておきます：
$$
\Phi = (\Phi_{i,j}), \ \ \Phi_{i,j} = \phi_j(x_i).
$$

1.2.2 計画行列

この記事では、基底関数としてガウシアン基底関数を用います。具体的には、次のように定めます(PRMLの(3.4)式)。
$$
\begin{align}
&{} \phi_0 (x) = 1 \
&{} \phi_j(x) = \exp\left[ -\frac{(x-\mu_j)^2}{2 s^2} \right] \ (j = 1, \dots, M-1),
\end{align}
$$
ここで、$\mu_j$は各ガウシアンの中心、$s$はガウシアンの幅を表します。当面の間は、これらの値は事前に適当に与えられているとします。

1.3 ベイジアン線型回帰

1.3.1 パラメーターの事後分布

前節のモデルより、尤度関数は次で与えられます(PRMLの(3.10)式)：
$$
\begin{align}
p(t_1, \dots, t_N | w, \beta) = \prod_{n=0}^{N-1} \mathcal{N}(t_n|w^T\phi(x_n), \beta^{-1}).
\end{align}
$$

パラメーター$w$の事前分布としては、次のガウス分布を用います(PRMLの(3.48)式)。
$$
\begin{align}
p(w|\alpha) = \mathcal{N}(w|m_0,S_0),
\end{align}
$$
ただし、$m_0 \in \mathbb{R}^M$はpriorの期待値で、$S_0$は$(M,M)$正定値実対称行列とします。

この設定の下で、パラメーター$w$の事後分布は次の分布で与えられます(PRMLの(3.49)-(3.51)式):
$$
\begin{align}
&{} p(w|t_1, \dots , t_n) = \mathcal{N}(w|m_N,S_N) \
&{} m_N := S_N (S_{0}^{-1}m_0 + \beta \Phi^T t) \
&{} S_N := \left( S_{0}^{-1} + \beta \Phi^T \Phi \right)^{-1} \
&{} \Phi = (\Phi_{i,j}), \ \ \Phi_{i,j} = \phi_j(x_i)
\end{align}
$$
特に、$m_0 = 0$, $S_0 = \frac{1}{\alpha} I$ ($\alpha > 0$)とおくと、
$$
\begin{align}
&{} m_N := \beta S_N \Phi^T t \
&{} S_N := \left( \alpha I + \beta \Phi^T \Phi \right)^{-1}
\end{align}
$$
となります。この後は、PRMLに従い、この形を仮定して進めていきます。

1.3.2 予測分布

$w$の事後分布を用いて$w$を積分して消去すると、$t$の予測分布を得ることができます(PRMLの(3.57)-(3.59)式)
$$
\begin{align}
p(t| x, \boldsymbol{t}, \alpha, \beta)
& := \int dw \ p(t|x, w,\beta) p(w | \boldsymbol{t}, \alpha, \beta) \
& = \mathcal{N}(t| m_{N}^{T} \phi(x), \sigma_{N}^{2}(x)) ,
\end{align}
$$
ただし、分散$\sigma_{N}^{2}(x)$は次で与えられます
$$
\begin{align}
\sigma_{N}^{2}(x) = \frac{1}{\beta} + \phi(x)^T S_N \phi(x)
\end{align}
$$

1.4 エビデンス近似

ここまでは、ハイパーパラメーター$\alpha, \beta$の値は与えられているとして話を進めてきました。ですが、実際に手でハイパーパラメーターの値を選ぶのは大変です。不適切なハイパーパラメーターを選ぶと、回帰の結果が望ましくないものになってしまうこともあります(後で実例を見ます。)。

そこで、できれば自動でハイパーパラメーターの値を決めてもらえると嬉しいところ。これをデータから行う手法の1つがエビデンス近似(または経験ベイズ)です。

1.4.1 エビデンス(周辺尤度)

エビデンス近似では、エビデンス$p(\boldsymbol{t}|\alpha,\beta)$を最大化するような$\alpha, \beta$を選択します。この記事で扱っているモデルについては、エビデンスの具体形は次で与えられます(PRMLの(3.86)式):
$$
\begin{align}
\log p(\boldsymbol{t}|\alpha,\beta)
&= \frac{M}{2} \log \alpha + \frac{N}{2} \log \beta - E(m_N) - \frac{1}{2} \log\left|\alpha I_M + \beta \Phi^T \Phi\right| - \frac{N}{2} \log ( 2\pi) \
E(m_N) &:=
\frac{\beta}{2} | \boldsymbol{t} - \Phi m_N |^2 + \frac{\alpha}{2} m_{N}^{T} m_{N}
\end{align}
$$

ここで、$m_N$は$\alpha, \beta$に依存することに注意!

1.4.2 最適化

エビデンスの対数を微分することにより、以下の停留条件が求まります(PRMLの(3.91), (3.92), (3.95)式) :
$$
\begin{align}
\lambda_i &: \ \beta\Phi^T \Phi \mbox{ の固有値 } \
\gamma &:= \sum_{i=1}^{M} \frac{\lambda_i}{ \lambda_i + \alpha} \
\alpha &= \frac{\gamma}{m_{N}^{T} m_{N} } \
\beta &= \frac{N-\gamma}{| \boldsymbol{t} - \Phi m_N |^2}
\end{align}
$$
PRMLの本文でも述べられている通り、この方程式をiterativeに解いていきます。

おさらいは以上です!

2 数式からコードへ

さあ、ここからはいよいよ数式をコードに落とし込んでいきます。

この記事では、

• 計画行列を返す函数gen_desmat_gaussian
• 線型回帰とエビデンス近似を行うBayesianRidgeRegression class
  の2つを作ることにしましょう。

2.1 計画行列を返す函数

訓練データ入力$x_0, \dots, x_{N-1}$を表す配列Xとparameter(ここでは$s$と$\mu_j$)を表す辞書paramsを入力とし、計画行列Phiを返す函数とします。
ここで、XはX[n, i]が$x_n$の第$i$成分を表すとしましょう。絵(?)で描くと
$$
X = \begin{pmatrix}
x_{0}^{T} \
x_{1}^{T} \
\vdots \
x_{N-1}^{T}
\end{pmatrix}
$$
という感じですね ¹。Phiについても同様とします(というよりも、前節で書いた$\Phi$の定義と同じ。)。

また、パラメーター$s$, $\mu_j$は直接渡しても良いのですが、他のパラメーターを渡すこともあるかなと思い、辞書型のparamsとしておきます。
params[mus]は(M-1, d) arrayで、params[mus][j, i]は$\mu_j$の第$i$成分を表すとします。

さて、ここから実際の計算です。「pythonではforやwhileのloopを直接書くと遅いので、numpy arrayを利用して計算をしましょう」という注意事項は、様々なところで口を酸っぱくして言われます(たとえばこちらの記事など)。　ここでの計算もこの路線で行きますが、ややトリッキーなので、詳しく見てみましょう。

$\Phi$の$(n, 0)$成分は全て1なので, $(n, m+1)$成分($n=0,\dots, N-1$, $m=0,\dots, M-2$)を考えることにします。
$\Phi_{n, m+1}$は以下のように表せます(ここで$X$は'X'、$\mathcal{M}$はparams['mus']を表します)

    \Phi_{n,m+1} = 
    \exp\left[ -\frac{1}{2 s^2} \sum_{i=0}^{d-1} \left(X_{n,i} - \mathcal{M}_{m,i} \right)^2 \right]

指数函数の部分はnp.expを成分ごとに作用させれば良いので、指数函数の中の部分を考えれば十分です。これを$A_{n,m}$と書くことにすると、

    A_{n,m} 
    := \sum_{i=0}^{d-1} \left( X_{n,i} - \mathcal{M}_{m,i} \right)^2 \\
    = \sum_{i=0}^{d-1} X_{n,i}^{2} + \sum_{i=0}^{d-1} \mathcal{M}_{m,i}^{2} 
 - 2 \sum_{i=0}^{d-1} X_{n,i} \mathcal{M}_{m,i}

となります。
これを、numpyのbroadcastingを利用して次のように計算します(下のコードブロックも参照ください。):

• 1: まず、-2*(X@(mus.T))で、上式の第3項に対応する(N,M-1) arrayを作ります。
• 2: 次に、np.reshape(np.sum(X**2, axis=1), (len(X), -1)) で第1項に対応する(N,1) arrayをつくり、先ほどのarrayと和をとります。このとき、numpy broadcastingによってarrayの形が調整されます。
• 3: 最後に、np.sum(mus**2, axis=1)で第2項に対応する(M-1,) arrayをつくり、和を取ります。

def gen_desmat_gaussian(X, params):    
    """
    This function generates a design matrix with basis function being Gaussian.
    
    Parameters
    ----------
    X :  2-D numpy array
        (N,d) numpy array, with X[n, i] = i-th element of x_n
    
    params : dictionary
        Dictionary of parameters of Gaussian basis function with
        params['mus'] : (M-1,d) numpy array. params['mus'][j] being the center of the jth Gaussian
        params['s'] : double. positive real number, which stands for the width of Gaussians
    
    Returns
    ----------
    Phi: 2-D numpy array
        (N,M) array, with Phi[n, m] = $\phi_m(x_n)$
    
    """

    s = params['s']
    mus = params['mus']
    Phi = np.zeros((len(X),len(mus)+1))
    Phi[:,0] = np.ones(len(X)) # the 0-the basis is a constant
    A = ( -2*(X@(mus.T)) + np.reshape(np.sum(X**2, axis=1), (len(X), -1)) ) + np.sum(mus**2, axis=1)
    Phi[:,1:] = np.exp(-A/(2*s*s))
    return Phi

計画行列の部分はこれで完了です。

2.2 線型回帰を行うクラス

次に、線型回帰を行うためのBayesianRidgeRegressionclassを作りましょう。

このクラスのオブジェクトには、次の属性を持たせることにします:

• ハイパーパラメーター$\alpha$(self.alpha)、
• ハイパーパラメーター$\beta$(self.beta)、
• 事後分布の期待値$m_N$ (self.m)
• 事後分布の共分散行列$S_N$ (self.S)

メソッドとしては

• コンストラクタ(__init__)
• 事後分布の期待値と共分散行列の計算(calc_posterior_params)
• 予測分布の期待値と標準偏差の計算(predict)
• エビデンスの値の計算(calc_evidence)
• ハイパーパラメーター$\alpha$, $\beta$の最適化を行うメソッド('empirical_bayes')

を持たせます。

以下では、各メソッドを1つずつ書き下していきます。

2.2.1 コンストラクタ

属性は上でも述べたとおりです。alphaとbetaを指定できるようにしておきます。

    def __init__(self, alpha=1.0, beta=1.0):
        self.alpha = alpha
        self.beta = beta
        self.m = None # posterior mean
        self.S = None # posterior covariancematrix

2.2.2 事後分布のパラメーターを計算

訓練データの入力値を表すarray Phiと出力値を表すarray tを与えられたときに、self.mとself.Sを計算します。
self.mとself.Sは次の式から計算できますので、これをそのままコードに翻訳します。
$$
\begin{align}
&{} m_N := \beta S_N \Phi^T t \
&{} S_N := \left( \alpha I + \beta \Phi^T \Phi \right)^{-1}
\end{align}
$$
Python 3.5から@が行列の積の演算子として使えるようになったので、これを活用します。

 
    def calc_posterior_params(self, Phi, t):
        """
        This method calculates posterior mean and covariance matrix from the training data Phi and t.
        
        Parameters
        ----------
        Phi : 2-D numpy array
            (N,M) array, representing design matrix
        t : 1-D numpy arra
            (N,) array, representing target values
        """

        self.S = np.linalg.inv(self.alpha*np.identity(len(Phi[0])) + self.beta*(Phi.T)@Phi )
        self.m = self.beta * ( self.S @ (Phi.T) @ t)

2.2.3 予測分布の計算

与えられた入力$\xi$に対して²、予測分布は期待値$m_{N}^{T} \phi(\xi)$, 標準偏差$\sigma_{N}(\xi)$のガウス分布となるので, この2つを計算してやればOKです。ただし、
$$
\begin{align}
\sigma_{N}^{2}(\xi) = \frac{1}{\beta} + \phi(\xi)^T S_N \phi(\xi)
\end{align}
$$
です.

通常、「1つの$\xi$についてだけ予測値(分布)がほしい」というケースは少なく、複数の入力、$\xi_0, \xi_1, \dots, \xi_{N_{test}-1}$について同時に予測分布を知りたいということが多々あります。
そこで、この実装では入力PhiはPhi[n, j] = $\phi_j(\xi_{n})$を表すarrayとして、これらの入力に対して予測分布(のパラメーター)をまとめて返します

pred_mean[n] = $m_{N}^{T} \phi(\xi_n)$
= $\sum_{i}$ Phi[n, i] * self.m[i] = (Phi @ self.m )[n]

pred_std[n] = $\sigma_{N}(\xi_n)$ = $\frac{1}{\beta} + \phi(\xi_n)^T S_N \phi(\xi_n)$
= 1/self.beta + $\sum_{i,j}$ Phi[n, i] * self.S[i, j] * Phi[n, j]
= 1/self.beta + (Phi @ self.S @ (Phi.T))[n, n]

return_stdがTrueのときのみpred_stdを返すようにしておきます(scikit-learnのBayesianRidgeをまねています。)。

    def predict(self, Phi, return_std=False):
        """
        This method makes prediction on the input Phi, and returns predictive mean (and standard deviation)
        
        Parameters
        ----------
        Phi : 2-D numpy array
            (N_test, M) numpy array. M must be equal to "M" (the length in the second dimension) of the training data.
        return_std : boolean, default False
            If True, the method also returns predictive standard deviation
            
        Returns 
        ----------
        pred_mean : 1-D numpy array
            (N_test,) numpy array representing predictive mean
        pred_std : 1-D numpy array
            (N_test,) numpy array representing predictive mean
        """

        pred_mean = Phi @ self.m
        if not(return_std):
            return pred_mean
        else:
            pred_std = np.sqrt(1.0/self.beta + np.diag(Phi @ self.S @ (Phi.T) ) )
            return pred_mean, pred_std

2.2.4 エビデンス近似

まず、エビデンスの対数値を計算します。これは数式を愚直にコードに焼直せば大丈夫です。

$$
\begin{align}
\log p(\boldsymbol{t}|\alpha,\beta)
&= \frac{M}{2} \log \alpha + \frac{N}{2} \log \beta - E(m_N) - \frac{1}{2} \log\left|\alpha I_M + \beta \Phi^T \Phi\right| - \frac{N}{2} \log ( 2\pi) \
E(m_N) &:=
\frac{\beta}{2} | \boldsymbol{t} - \Phi m_N |^2 + \frac{\alpha}{2} m_{N}^{T} m_{N}
\end{align}
$$

    def calc_evidence(self, Phi, t):
        """
        This method calculates the evidence with respect to the data Phi and t
        
        Parameters
        ----------
        Phi : 2-D numpy array
            (N,M) array, representing design matrix
        t : 1-D numpy arra
            (N,) array, representing target values
            
        Returns
        ----------
        evidence : float
        
        """

        N, M = np.shape(Phi)
        evidence = 0.5*M*np.log(self.alpha) + 0.5*N*np.log(self.beta) \
            - 0.5*self.beta*np.linalg.norm( t - Phi @ self.m )**2 - 0.5*self.alpha*(self.m@self.m) \
            - 0.5*np.log( np.linalg.det( self.alpha*np.identity(M) + self.beta*(Phi.T)@Phi ) ) \
            - 0.5*N*np.log(2*np.pi)
        return evidence

この値そのものは$\alpha$と$\beta$の最適化には利用しないものの、他のハイパーパラメーターの最適化に利用することができます(次の実験の節で具体例を見ます。)。

次に、$\alpha$と$\beta$についての最適化を行うメソッドを書きます。
具体的には、以下の方程式の解を、代入を繰り返すことで求めます。
$$
\begin{align}
\lambda_i &: \ \beta\Phi^T \Phi \mbox{　の固有値 } \
\gamma &:= \sum_{i=1}^{M} \frac{\lambda_i}{ \lambda_i + \alpha} \
\alpha &= \frac{\gamma}{m_{N}^{T} m_{N} } \
\beta &= \frac{N-\gamma}{| \boldsymbol{t} - \Phi m_N |^2}
\end{align}
$$

$\lambda_i$は、最初に$\Phi^T \Phi$の固有値を求めておけば$\beta$をかけるだけでよいので、最初に$\Phi^T \Phi$を対角化しておきます。

    def empirical_bayes(self, Phi, t, tol, maxiter, show_message=True):
        """
        This method performs empirical bayes (or evidence approximation), 
        where hyper parameters alpha and beta are chosen in such a way that they maximize the evidence. 
        
        Parameters
        ----------
        Phi : 2-D numpy array
            (N,M) array, representing design matrix
        t : 1-D numpy arra
            (N,) array, representing target values
        tol : float
            The tolerance. 
            If the changes of alpha and beta are smaller than the value, the iteration is judged as converged.
        maxiter : int
            The maximum number of iteration
        show_message : boolean, default True
            If True, the message indicating whether the optimization terminated successfully is shown.
        """

        tmp_lambdas = eigh((Phi.T)@Phi)[0]
        cnt = 0
        while cnt < maxiter:
            lambdas = self.beta * tmp_lambdas
            self.calc_posterior_params(Phi, t)
            
            alpha_old = self.alpha
            beta_old = self.beta
            
            gamma = np.sum( lambdas/ (self.alpha + lambdas) )
            self.alpha = gamma/np.dot(self.m, self.m)
            self.beta = (len(t) - gamma) / ( np.linalg.norm(t -  Phi @ self.m )**2    )
            if (abs(self.alpha - alpha_old) < tol) and ( abs(self.beta - beta_old) < tol ):
                break
            cnt += 1
        if show_message:
            if cnt <= maxiter:
                print(f"Optimization terminated succesfully. The number of iteration : {cnt}")
            else:
                print("Maximum number of iteration exceeded.")

2.2.5 fit

最後に、scikit-learnのestimatorのようにfitメソッドを書いておきます(といっても、今までのメソッドを利用するだけですが。)。optimize_hyperparamsがTrueのときには、事前にエビデンス近似を行って$\alpha$と$\beta$の値を選んでくれます。

    def fit(self, Phi, t, tol=1e-4, maxiter=100, show_message=True, optimize_hyperparams=False):
        """
        This method performs fitting. 
        The user can choose whether or not to perform empirical Bayes.
        
        Parameters
        ----------
        Phi : 2-D numpy array
            (N,M) array, representing design matrix
        t : 1-D numpy arra
            (N,) array, representing target values
        tol : float
            The tolerance. 
            If the changes of alpha and beta are smaller than the value, the iteration is judged as converged.
        maxiter : int
            The maximum number of iteration
        show_message : boolean, default True
            If True, the message indicating whether the optimization terminated successfully is shown.
        optimize_hyperparams : boolean, default False
            If True, the hyper parameters alpha and beta are optimized by empirical Bayes.
        """

        if optimize_hyperparams:
            self.empirical_bayes(Phi, t, tol, maxiter, show_message)
        self.calc_posterior_params(Phi, t)            

2.2.6 BayesianRidgeRegression

まとめると、BayesianRidgeRegressionのコードは次の通りです。

class BayesianRidgeRegression:
    
    def __init__(self, alpha=1.0, beta=1.0):
        self.alpha = alpha
        self.beta = beta
        self.m = None # posterior mean
        self.S = None # posterior covariancematrix
        
    
    def calc_posterior_params(self, Phi, t):
        """
        This method calculates posterior mean and covariance matrix from the training data Phi and t.
        
        Parameters
        ----------
        Phi : 2-D numpy array
            (N,M) array, representing design matrix
        t : 1-D numpy arra
            (N,) array, representing target values
        """
        self.S = np.linalg.inv(self.alpha*np.identity(len(Phi[0])) + self.beta*(Phi.T)@Phi )
        self.m = self.beta * ( self.S @ (Phi.T) @ t)


    def predict(self, Phi, return_std=False):
        """
        This method makes prediction on the input Phi, and returns predictive mean (and standard deviation)
        
        Parameters
        ----------
        Phi : 2-D numpy array
            (N_test, M) numpy array. M must be equal to "M" (the length in the second dimension) of the training data.
        return_std : boolean, default False
            If True, the method also returns predictive standard deviation
            
        Returns 
        ----------
        pred_mean : 1-D numpy array
            (N_test,) numpy array representing predictive mean
        pred_std : 1-D numpy array
            (N_test,) numpy array representing predictive mean
        """
        pred_mean = Phi @ self.m
        if not(return_std):
            return pred_mean
        else:
            pred_std = np.sqrt(1.0/self.beta + np.diag(Phi @ self.S @ (Phi.T) ) )
            return pred_mean, pred_std
        
        
    def calc_evidence(self, Phi, t):
        """
        This method calculates the evidence with respect to the data Phi and t
        
        Parameters
        ----------
        Phi : 2-D numpy array
            (N,M) array, representing design matrix
        t : 1-D numpy arra
            (N,) array, representing target values
            
        Returns
        ----------
        evidence : float
        
        """
        N, M = np.shape(Phi)
        evidence = 0.5*M*np.log(self.alpha) + 0.5*N*np.log(self.beta) \
            - 0.5*self.beta*np.linalg.norm( t - Phi @ self.m )**2 - 0.5*self.alpha*(self.m@self.m) \
            - 0.5*np.log( np.linalg.det( self.alpha*np.identity(M) + self.beta*(Phi.T)@Phi ) ) \
            - 0.5*N*np.log(2*np.pi)
        return evidence
                
    
    def empirical_bayes(self, Phi, t, tol, maxiter, show_message=True):
        """
        This method performs empirical bayes (or evidence approximation), 
        where hyper parameters alpha and beta are chosen in such a way that they maximize the evidence. 
        
        Parameters
        ----------
        Phi : 2-D numpy array
            (N,M) array, representing design matrix
        t : 1-D numpy arra
            (N,) array, representing target values
        tol : float
            The tolerance. 
            If the changes of alpha and beta are smaller than the value, the iteration is judged as converged.
        maxiter : int
            The maximum number of iteration
        show_message : boolean, default True
            If True, the message indicating whether the optimization terminated successfully is shown.
        """
        tmp_lambdas = eigh((Phi.T)@Phi)[0]
        cnt = 0
        while cnt < maxiter:
            lambdas = self.beta * tmp_lambdas
            self.calc_posterior_params(Phi, t)
            
            alpha_old = self.alpha
            beta_old = self.beta
            
            gamma = np.sum( lambdas/ (self.alpha + lambdas) )
            self.alpha = gamma/np.dot(self.m, self.m)
            self.beta = (len(t) - gamma) / ( np.linalg.norm(t -  Phi @ self.m )**2    )
            if (abs(self.alpha - alpha_old) < tol) and ( abs(self.beta - beta_old) < tol ):
                break
            cnt += 1
        if show_message:
            if cnt <= maxiter:
                print(f"Optimization terminated succesfully. The number of iteration : {cnt}")
            else:
                print("Maximum number of iteration exceeded.")
                
    def fit(self, Phi, t, tol=1e-4, maxiter=100, show_message=True, optimize_hyperparams=False):
        """
        This method performs fitting. 
        The user can choose whether or not to perform empirical Bayes.
        
        Parameters
        ----------
        Phi : 2-D numpy array
            (N,M) array, representing design matrix
        t : 1-D numpy arra
            (N,) array, representing target values
        tol : float
            The tolerance. 
            If the changes of alpha and beta are smaller than the value, the iteration is judged as converged.
        maxiter : int
            The maximum number of iteration
        show_message : boolean, default True
            If True, the message indicating whether the optimization terminated successfully is shown.
        optimize_hyperparams : boolean, default False
            If True, the hyper parameters alpha and beta are optimized by empirical Bayes.
        """
        if optimize_hyperparams:
            self.empirical_bayes(Phi, t, tol, maxiter, show_message)
        self.calc_posterior_params(Phi, t)            

3 実験

ここまで延々と書いてきましたが、いよいよ実験です!

3.1 データ

使うデータはこちら
$$
\begin{eqnarray}
&{}& t = f(x) + \varepsilon \
&{}& f(x) = \sin(2x) + 0.2\sin x + 0.1x \
&{}& \varepsilon \sim \mathcal{N}(0, 0.09)
\end{eqnarray}
$$
データ数は50です。

3.2 ベイジアン線型回帰

では、まずはハイパーパラメーターを手で決めてベイジアン線型回帰をやってみましょう。
計画行列のパラメーターは

• params['s'] = 0.5
• params['mus'] = np.reshape(np.linspace(-3,3,16),(16,1))
  で与えておきます。

まず、default値のalpha=1.0, beta=1.0の場合です

水色の領域は、予測分布期待値±予測分布の標準偏差の領域を表します。
まあまあ良い感じにfitしているのが見て取れるかと思います。

ここからalphaの値を一気に大きくするとどうなるか見てみましょう:

予測分布の期待値が全然あわなくなっていますね(その分、credible intervalは広くなっています。)。
$\alpha$を大きくするというのは事前分布でweight $w$が小さくなるようにしているということなので、得られる予測分布の期待値の「振幅」が小さくなるのは直感的に納得できます。

次に、sを小さくするとどうなるかというと...

なんじゃこりゃ、という感じですが、これも直感的には納得しやすい結果です。というのも、$s$を小さくするのは規定関数の幅を小さくすることに対応するからです。

3.3 エビデンス近似

次に、ハイパーパラメーターを自動で調節してもらいましょう。

3.3.1 $\alpha$と$\beta$についての最適化

まず、alphaとbetaについての最適化です。これはfitメソッドでoptimize_hyperparams=Trueとしておけば、自動でやってもらえます。

全体として、予測分布の期待値は真の曲線に近く、また、credible intervalも前の例に比べて小さくなっています。
得られた値は、alpha=4.4, beta=12.9でした。

3.3.2 sについての最適化

次に、計画行列の取り方についての最適化を考えましょう。
原理的には、$s$と$\mu_j$を動かしてエビデンスを最大化するものをとればよいのですが、ここでは簡単のため、$s$のみを考えます。解析的な計算はできないので、$s$についてのgrid searchを行いました。
その結果がこちら

値はおおよそs=0.61、alpha=3.6、beta=13.0でした。
さきほどと大きくは変わらないので、差はぱっと見るではあまりないですかね。。。

4 まとめ

この記事では、ベイジアン線型回帰とそのエビデンス近似を実装しました。
ベイジアン的な手法を用いることにより、予測値だけでなくその分布を求めることができました。また、エビデンス近似を用いることにより、ハイパーパラメーターをデータから決めることができました。

ただ、今回実装した線型回帰のアプローチは、(PRMLの3.6節でも述べられているように)大きな弱点があります。それは、入力の次元$d$が大きくなったとき、基底関数の個数を$d$に対して指数的に増やしていく必要があることです³。この問題を回避できるSupport vector machineやNeural networkなどの手法は、5章以降で取り扱っていきます。

1. ここの添え字の並びを、[データ点の番号, 特徴量の番号]か[特徴量の番号, データ点の番号]どちらにするかは任意性がありますが、scikit-learnなどは前者を前提としているので、この方式を採用します。Andrew NgのCourseraでのDeep Learningの講義では後者だったと思います。) ↩

2. 訓練データと区別するために, $x$ではなく$\xi$と書くことにします。 ↩

3. Gaussian既定関数の例を考えると: 幅$s$のGaussianの「山」が一辺の長さ$L$の$d$次元の立方体を「覆い尽くす」ためには、おおよそ$(L/s)^d$個のGaussianが必要になります。これは$d$に対して指数的に増えますね。 ↩