🎃
ベイズ線形回帰とエビデンス近似
はじめに
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{w}
尤度関数と事前分布は以下の通りです。
よって、事後分布は次のようになります。
予測分布とエビデンス近似
改めて、予測分布は以下のような形で書けます。
しかし、この積分を計算するのは困難です。なので、まず
\alpha, \beta
この対数事後分布を
トイデータで学習
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個の基底関数を作ります。ガウス基底関数やシグモイド基底関数も使っています。今回は
\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}
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}
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を計算
今回は二乗誤差を考えます。この誤差を最小にする値は、
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回帰があります。これは今回扱ったモデルに近いですが、以下の式を最小にする
しかし、適切な
ハイパーパラメータの適切さを検証
様々な\alpha, \beta
様々な
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()
この図を見ると周辺尤度は単峰性で、分散もある程度小さいことが見てとれるので、エビデンス近似によって本来の予測分布をうまく近似できているでしょう。
様々な\lambda = \frac{\alpha}{\beta}
様々な
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()
赤い線が
Discussion