🎃

ベイズ線形回帰とエビデンス近似

2022/03/19に公開約13,000字

はじめに

PRML(Pattern Recognition and Machine Learning)でベイズ線形回帰について学んだ内容をまとめて、実際のデータを使って学習しました。主に3章の内容です。学習アルゴリズムの全体像を示すことを意識したので、導出は端折っている部分が多いです。詳しい解の導出は、PRMLの本文や演習問題にあります。

ベイズ線形回帰のアルゴリズム

変数一覧

  • 訓練データの数 N
  • 基底関数の数 M
  • 訓練データの説明変数 \mathbf{X} = \left\{\mathbf{x_1}, \mathbf{x_2}, \cdots \mathbf{x_N}\right\}, \mathbf{x_n} = (x_{n, 1}, x_{n, 2}, \cdots, x_{n_D})^T
  • 訓練データの目的変数 \mathbf{t} = \left\{t_1, t_2, \cdots t_N\right\}
  • 基底関数 \mathbf{\phi(x_n)} = (\phi_0(\mathbf{x_n}), \phi_1(\mathbf{x_n}), \cdots \phi_{M-1}(\mathbf{x_n}))^T
  • 重み \mathbf{w} = (w_0, w_1, \cdots, w_{M-1})^T
  • 新たな説明変数 \mathbf{x}

また以下の行列を使います。

\mathbf{\Phi} = \begin{pmatrix} \phi_0(\mathbf{x_1}) & \phi_1(\mathbf{x_1}) & \cdots & \phi_{M-1}(\mathbf{x_1}) \\ \phi_0(\mathbf{x_2}) & \phi_1(\mathbf{x_2}) & \cdots & \phi_{M-1}(\mathbf{x_2}) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_0(\mathbf{x_N}) & \phi_1(\mathbf{x_N}) & \cdots & \phi_{M-1}(\mathbf{x_N}) \\ \end{pmatrix}

目的と手段

最終的な目標は以下の予測分布を求めることです。訓練データから、新たな説明変数に対する目的変数の分布を求めます。

\begin{aligned} p(t | \mathbf{x}, \mathbf{t}, \mathbf{X}) \end{aligned}

今回は以下のような線形回帰モデルを使用します。

\begin{aligned} p(t | \mathbf{x}, \mathbf{w}, \beta) = N(t | \mathbf{w}^T\mathbf{\phi(\mathbf{x}}), \beta^{-1}) \end{aligned}

これ以降、見やすさのために\mathbf{x}\mathbf{X}は省略します。

\mathbf{w}の事後分布

\mathbf{w}についての事後分布は次のように表せます。

\begin{aligned} p(\mathbf{w} | \mathbf{t}) \propto p(\mathbf{t} | \mathbf{w})p(\mathbf{w}) \end{aligned}

尤度関数と事前分布は以下の通りです。\alpha, \betaはハイパーパラメータです。

\begin{aligned} p(\mathbf{t} | \mathbf{w}) &= N(\mathbf{t} | \mathbf{\Phi}\mathbf{w}, \beta^{-1}\mathbf{I}) \\\\ p(\mathbf{w}) &= N(\mathbf{w} | \mathbf{0}, \alpha^{-1}\mathbf{I}) \end{aligned}

よって、事後分布は次のようになります。

\begin{aligned} p(\mathbf{w} | \mathbf{t}) &= N(\mathbf{w} | \mathbf{m_N}, \mathbf{S_N}) \\\\ where \quad \mathbf{m_N} &= \beta\mathbf{S_N}\mathbf{\Phi^T}\mathbf{t} \\\\ \mathbf{S_N}^{-1} &= \alpha\mathbf{I} + \beta\mathbf{\Phi^T}\mathbf{\Phi} \end{aligned}

予測分布とエビデンス近似

改めて、予測分布は以下のような形で書けます。

\begin{aligned} p(t | \mathbf{t}) = \int\int\int p(t | \mathbf{w}, \beta) p(\mathbf{w} | \mathbf{t}, \alpha, \beta) p(\alpha, \beta | \mathbf{t}) d\mathbf{w}d\alpha d\beta \end{aligned}

しかし、この積分を計算するのは困難です。なので、まず\mathbf{w}することで周辺化して得られた周辺尤度における\alpha, \betaの最尤推定量\hat{\alpha}, \hat{\beta}を求めて、その値で固定することで近似的に予測分布を求めます。この手法をエビデンス近似といいます。

\begin{aligned} p(t | \mathbf{t}) &\simeq p(t | \mathbf{t}, \hat{\alpha}, \hat{\beta}) \\\\ &= \int p(t | \mathbf{w}, \hat{\beta}) p(\mathbf{w} | \mathbf{t}, \hat{\alpha}, \hat{\beta}) d\mathbf{w} \\\\ &= \int N(t | \mathbf{w}^T\mathbf{\phi(\mathbf{x}}), \hat{\beta}^{-1}) N(\mathbf{w} | \mathbf{m_N}, \mathbf{S_N}) d\mathbf{w} \\\\ &= N(t | \mathbf{m_N}^T\mathbf{\phi(x)}, {\sigma_N(\mathbf{x})}^2) \\\\ where \quad \sigma_N(\mathbf{x})^2 &= \frac{1}{\hat{\beta}} + \mathbf{\phi(x)}^T\mathbf{S_N}\mathbf{\phi(x)} \end{aligned}

\alpha, \betaの最尤推定

\alpha, \betaの事後分布と、対数事後分布を求めます。ここでは、事前分布は一定と考えます。

\begin{aligned} p(\alpha, \beta | \mathbf{t}) &\propto p(\mathbf{t} | \alpha, \beta) p(\alpha, \beta) \\\\ &\propto p(\mathbf{t} | \alpha, \beta) \\\\ &= \int p(\mathbf{t} | \mathbf{w}, \beta)p(\mathbf{w} | \alpha) d\mathbf{w} \\\\ &= \int N(\mathbf{t} | \mathbf{\Phi}\mathbf{w}, \beta^{-1}\mathbf{I}) N(\mathbf{w} | \mathbf{0}, \alpha^{-1}\mathbf{I}) d\mathbf{w} \\\\ &= \left(\frac{\beta}{2\pi}\right)^{\frac{N}{2}} \left(\frac{\alpha}{2\pi}\right)^{\frac{M}{2}} \int exp\{-\frac{\beta}{2}||\mathbf{t} - \mathbf{\Phi}\mathbf{w}||^2 - \frac{\alpha}{2}\mathbf{w}^T\mathbf{w}\} d\mathbf{w} \\\\ &= \alpha^{\frac{M}{2}} \beta^{\frac{N}{2}} exp\{-E(\mathbf{m_N})\} |\mathbf{A}|^{-\frac{1}{2}} (2\pi)^{-\frac{N}{2}} \\\\ lnp(\mathbf{t} | \alpha, \beta) &= \frac{M}{2}ln\alpha + \frac{N}{2}ln\beta - E(\mathbf{m_N}) - \frac{1}{2}ln|\mathbf{A}| - \frac{N}{2}ln(2\pi) \\\\ where \quad \mathbf{A} &= \alpha\mathbf{I} + \beta\mathbf{\Phi^T}\mathbf{\Phi} \\\\ \mathbf{m_N} &= \beta^{-1}\mathbf{A}^{-1}\mathbf{\Phi}^T\mathbf{t} \\\\ E(\mathbf{m_N}) &= \frac{\beta}{2}||\mathbf{t} - \mathbf{\Phi}\mathbf{m_N}||^2 + \frac{\alpha}{2}\mathbf{m_N}^T\mathbf{m_N} \end{aligned}

この対数事後分布を\alpha, \betaそれぞれで微分して、\alpha, \betaの最尤推定による解を求めた結果、以下のようになります。ただし、\lambda_i\beta\mathbf{\Phi}^T\mathbf{\Phi}の固有値です。

\begin{aligned} \hat{\alpha} &= \frac{\gamma}{\mathbf{m_N}^T\mathbf{m_N}} \\\\ \frac{1}{\hat{\beta}} &= \frac{1}{N-\gamma}\sum^N_{n=1}\{t_n - \mathbf{m_N}^T\mathbf{\phi(x_n)}\}^2 \\\\ where \quad \gamma &= \sum_i\frac{\lambda_i}{\hat{\alpha} + \lambda_i} \end{aligned}

トイデータで学習

import文

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split

np.random.seed(28)

今回使用するデータセット

scikit-learnの「糖尿病データセット」というトイデータセットを使います。このデータセットでは、年齢、性別、BMI、平均血圧、6つの血清検査測定値という10の説明変数から、1年後の糖尿病の進行状況を予測します。データ数は442です。このうちの8割に当たる353個を訓練用データとし、残りをテスト用データとします。また、データは標準化しておきます。

input_variables = ['age', 'sex', 'bmi', 'map', 'tc', 'ldl', 'hdl', 'tch', 'ltg', 'glu']

data_url = 'https://www4.stat.ncsu.edu/~boos/var.select/diabetes.tab.txt'
raw_df = pd.read_csv(data_url, sep="\s+")
diabetes_df = raw_df.set_axis(input_variables + ['y'], axis=1)

y_origin_df = diabetes_df.loc[:, 'y']
y_df = (y_origin_df - y_origin_df.mean()) / y_origin_df.std()
x_origin_df = diabetes_df.loc[:, input_variables]
x_df = (x_origin_df - x_origin_df.mean()) / x_origin_df.std()

train_x_df, test_x_df, train_y_df, test_y_df = train_test_split(x_df, y_df, test_size=0.2)

train_x = train_x_df.to_numpy(copy=True)
train_y = train_y_df.to_numpy(copy=True)
test_x = test_x_df.to_numpy(copy=True)
test_y = test_y_df.to_numpy(copy=True)

diabetes_df.head()
index age sex bmi map tc ldl hdl tch ltg glu y
0 59 2 32.1 101.0 157 93.2 38.0 4.0 4.8598 87 151
1 48 1 21.6 87.0 183 103.2 70.0 3.0 3.8918 69 75
2 72 2 30.5 93.0 156 93.6 41.0 4.0 4.6728 85 141
3 24 1 25.3 84.0 198 131.4 40.0 5.0 4.8903 89 206
4 50 1 23.0 101.0 192 125.4 52.0 4.0 4.2905 80 135

基底関数

296個の基底関数を作ります。ガウス基底関数やシグモイド基底関数も使っています。今回はs=0.5としています。

  • \phi_0(\mathbf{x_n}) = 1
  • \phi_j(\mathbf{x_n}) = x_{n, k}x_{n, l} \quad(1 \leq j \leq 45)
  • \phi_j(\mathbf{x_n}) = x_{n, k} \quad(46 \leq j \leq 55)
  • \phi_j(\mathbf{x_n}) = x_{n, k}^2 \quad(56 \leq j \leq 65)
  • \phi_j(\mathbf{x_n}) = x_{n, k}^3 \quad(66 \leq j \leq 75)
  • \phi_j(\mathbf{x_n}) = exp(-\frac{(x_{n, k} - \mu_{l})^2}{2s^2}) \quad(76 \leq j \leq 185)
  • \phi_j(\mathbf{x_n}) = \frac{1}{1 + exp(-\frac{x_{n, k} - \mu_{l}}{s})} \quad(186 \leq j \leq 295)
def init_base_func_list():
    base_func_list = [0 for i in range(45)]
    ind = 0
    for i in range(10):
        for j in range(i+1, 10):
            base_func_list[ind] = (i, j)
            ind += 1

    return base_func_list

M = 296
base_func_list = init_base_func_list()
s = 0.5

def base_func(x_n, m):
    if m == 0:
        return np.ones(len(x_n))
    elif m < 46:
        return x_n[:, base_func_list[m-1][0]] * x_n[:, base_func_list[m-1][1]]
    elif m < 56:
        return x_n[:, m-46]
    elif m < 66:
        return x_n[:, m-56]**2
    elif m < 76:
        return x_n[:, m-66]**3
    elif m < 186:
        m -= 76
        x_ind = m // 11
        gauss_mean = -1.5 + 0.3 * (m % 11)
        return np.exp(-(x_n[:, x_ind] - gauss_mean)**2 / (2 * s**2))
    elif m < 296:
        m -= 186
        x_ind = m // 11
        sigmoid_mean = -1.5 + 0.3 * (m % 11)
        return 1 / (1 + np.exp((sigmoid_mean - x_n[:, x_ind]) / s))

    raise ValueError()

\mathbf{\Phi}を計算

\mathbf{\Phi}は冒頭で示した通り、基底関数に説明変数を入力したものです。

def Phi(x, M):
    x_nm = np.zeros((len(x), M))
    
    for m in range(M):
        x_nm[:, m] = base_func(x, m)
    return x_nm

phi_x_nm = Phi(train_x, M)
phi_x_test = Phi(test_x, M)
phit_phi = np.dot(phi_x_nm.T, phi_x_nm)

\hat{\alpha}, \hat{\beta}を計算

\hat{\alpha}, \hat{\beta}の値を直接求めることはできないです。なので、\alpha, \betaに初期値を与えて、\gamma, \alpha, \betaの順に繰り返し値を更新することで推定します。

v_m, _ = LA.eig(phit_phi)
v_m = v_m.real
alpha = 1.0
beta = 1.0

for i in range(50):
    lambda_m = v_m * beta
    gamma = sum(lambda_m / (alpha + lambda_m))
    A_mm = alpha * np.identity(M) + beta * phit_phi
    m_m = beta * np.dot(np.dot(LA.inv(A_mm), phi_x_nm.T), train_y)
    alpha = gamma / np.dot(m_m, m_m)
    y_phi_m = train_y - np.dot(phi_x_nm, m_m)
    beta = (N - gamma) / np.dot(y_phi_m, y_phi_m)
    if (i+1) % 10 == 0:
        print("iteration, alpha, beta | {0:2d}, {1:.5f}, {2:.5f}".format(i+1, alpha, beta))

alpha_hat = alpha
beta_hat = beta
SN_inv = alpha_hat * np.identity(M) + beta_hat * phit_phi
mN = np.dot(LA.inv(SN_inv), np.dot(beta_hat * phi_x_nm.T, train_y))
iteration, alpha, beta | 10, 470.62081, 2.22515
iteration, alpha, beta | 20, 470.77612, 2.22506
iteration, alpha, beta | 30, 470.77612, 2.22506
iteration, alpha, beta | 40, 470.77612, 2.22506
iteration, alpha, beta | 50, 470.77612, 2.22506

しっかり収束していることがわかります。

RMSEを計算

今回は二乗誤差を考えます。この誤差を最小にする値は、E[t | \mathbf{x}] = \mathbf{m_N}^T\mathbf{\phi(x)}、つまり予測分布であるガウス分布の平均パラメータとなります。二乗平均平方根誤差(RMSE)を計算すると、以下のようになります。また、目的変数の分散をかけることで、標準化する前のRMSEを求めています。

def RMSE(phi, y, w):
    error_sum = sum((y - np.dot(phi, w)) ** 2) * y_origin_df.std() ** 2
    return np.sqrt(error_sum / len(y))

print(RMSE(phi_x_nm, train_y, mN))
print(RMSE(phi_x_test, test_y, mN))
print(mN[:5])
48.50338074738126
57.08358371053724
[-0.00518     0.04444021  0.00333912  0.02872613 -0.02135499]

テストデータのRMSEは57という結果になりました。

ベイズの特長

最尤法による線形回帰モデルの代表的なものにRidge回帰があります。これは今回扱ったモデルに近いですが、以下の式を最小にする\mathbf{w}を求めます。

\begin{aligned} E(\mathbf{w}) = \frac{1}{2} \sum^N_{n=1}\{t_n - \mathbf{w}^T\mathbf{\phi(x_n)}\}^2 + \frac{\lambda}{2}\mathbf{w}^T\mathbf{w} \end{aligned}

しかし、適切な\lambdaの値は分からないため、様々な\lambdaの値を実際に試してみて、検証用データの誤差によって適切な\lambdaの値を探る必要があります。Ridge回帰に限らないですが、適切なハイパーパラメータを探ることで計算量が増加したり、検証用データを学習に用いることができなくなったりします。ですが、今回はエビデンス近似を使ったことで、ハイパーパラメータである\alpha, \betaの適切な値をベイズの枠組みの中で決めることができました。

ハイパーパラメータの適切さを検証

様々な\alpha, \betaに対する対数周辺尤度の図

様々な\alpha, \betaに対して、対数周辺尤度の値を計算して図示しました。途中でlog|\mathbf{A}| = Tr(ln\mathbf{A})を使っていますが、これは|\mathbf{A}|の値が大きすぎて直接計算できないからです。

from scipy.linalg import logm

GRIDS = 31

alpha_line = np.linspace(alpha_hat - 300, alpha_hat + 300, num=GRIDS)
beta_line = np.linspace(beta_hat - 0.3, beta_hat + 0.3, num=GRIDS)
alpha_grid, beta_grid = np.meshgrid(alpha_line, beta_line)
ps = np.empty((GRIDS, GRIDS))

for a in range(GRIDS):
    for b in range(GRIDS):
        alpha = alpha_line[a]
        beta = beta_line[b]

        mN = beta * np.dot(LA.inv(SN_inv), np.dot(phi_x_nm.T, train_y))
        square_error = sum((train_y - np.dot(phi_x_nm, mN)) ** 2)
        E_mN = (beta/2) * square_error + (alpha/2) * np.dot(mN, mN)
        A = alpha * np.identity(M) + beta * phit_phi
        log_det_A = np.trace(logm(A))
        ps[b][a] = (M/2) * np.log(alpha) + (N/2) * np.log(beta) - E_mN - (1/2) * log_det_A
levels = [-77 - i for i in range(17)]
levels.reverse()

plt.figure(figsize=(6 * np.sqrt(2), 6))
plt.contour(alpha_grid, beta_grid, ps, levels=levels)  真の分布
plt.colorbar()
plt.scatter(x=alpha_hat, y=beta_hat, color='red', s=150, marker='x') # 学習したパラメータ
plt.xlabel('alpha')
plt.ylabel('beta')
plt.show()

likelihood

この図を見ると周辺尤度は単峰性で、分散もある程度小さいことが見てとれるので、エビデンス近似によって本来の予測分布をうまく近似できているでしょう。

様々な\lambda = \frac{\alpha}{\beta}に対するRMSEの図

様々な\alpha, \betaに対するテストデータのRMSEを求めます。ここで、\mathbf{m_N}^T\mathbf{\phi(x)}\lambda = \frac{\alpha}{\beta}の値に依存するので、\lambdaとRMSEの関係を図にしました。

NUM = 1000
lambda_line = np.linspace(1, 1000, num=NUM)
train_RMSEs = np.empty(NUM)
test_RMSEs = np.empty(NUM)

for lam in range(NUM):
    w_ML = np.dot(np.dot(LA.inv(lambda_line[lam] * np.identity(M) + phit_phi), phi_x_nm.T), train_y)
    train_RMSEs[lam] = RMSE(phi_x_nm, train_y, w_ML)
    test_RMSEs[lam] = RMSE(phi_x_test, test_y, w_ML)
plt.figure(figsize=(6*np.sqrt(2), 6))
plt.plot(lambda_line, train_RMSEs, label='train_RMSE')
plt.plot(lambda_line, test_RMSEs, label='test_RMSE')
plt.axvline(x=alpha/beta, color='red')
plt.xlabel('$\lambda$')
plt.ylabel('$RMSE$')
plt.legend()
plt.show()

RMSE

赤い線が\frac{\hat{\alpha}}{\hat{\beta}}の値です。赤い線の部分でテストデータのRMSEが、ある程度最小値に近い値をとっていることがわかります。ハイパーパラメータの学習はうまくいきました。

Discussion

ログインするとコメントできます