はじめに
本記事では標準ベイズ統計学の7章で掲載されている図表やモデルをPythonで実装する方法に関して説明します。
7章：多変量正規モデル
これまでの章では、一つの変数を分析する方法を学びました。しかし、この章では複数の変数（身長、体重、年齢など）を同時に分析する方法に踏み込みます。多変量正規モデルを使うことで、複数の変数の平均、分散（ばらつき）、そして変数間の相関（関係性）を一度に推定できるようになります。これは、データ全体の構造をより深く理解するのに役立ちます。また、現実のデータでは、一部の情報が欠けていることがよくあります。この章では、多変量正規モデルを使って、欠けているデータを統計的に適切な方法で推測する方法も紹介しています。
多変量正規分布からのサンプルと密度
図7.1は三つの異なる二次元正規分布の等高線図と各分布から30個のサンプルを表しています。各ケースにおいて、
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
# 多変量正規分布ベクトル
def rmvnorm(n, mu, Sigma):
return np.random.multivariate_normal(mu, Sigma, n)
# 多変量正規分布の対数密度
def ldmvnorm(y, mu, Sig):
return stats.multivariate_normal.logpdf(y, mean=mu, cov=Sig)
# 設定
mu = [50, 50]
sig = [8, 12]
ng = 50
y = np.linspace(20, 80, ng)
np.random.seed(1)
# プロット設定
fig, axes = plt.subplots(1, 3, figsize=(15, 10))
for idx, rho in enumerate([-0.5, 0, 0.5]):
Sigma = np.array([[1, rho], [rho, 1]]) * np.outer(sig, sig)
Y = rmvnorm(30, mu, Sigma)
LDY = np.zeros((ng, ng))
for i in range(ng):
for j in range(ng):
LDY[i, j] = ldmvnorm([y[i], y[j]], mu, Sigma)
ax = axes[idx]
ax.set_xlim([20, 80])
ax.set_ylim([20, 80])
ax.set_xlabel('y1', fontsize=10)
ax.set_ylabel('y2', fontsize=10)
ax.contourf(y, y, np.exp(LDY), cmap='Greys')
ax.scatter(Y[:, 0], Y[:, 1], marker='x', color='red')
plt.tight_layout()
plt.show()
多変量正規分布からのサンプルと密度
読解力データと事後分布
図7.2の解説に入る前に、多変量正規分布のギブスサンプリングについて解説します。
標本データが単変量正規分布に従う時、母集団平均の事前分布もまた単変量正規分布に従うことは以前の章で解説しました。同様に、多変量正規分布の平均の事前分布も以下のような多変量正規分布に従います。
標本データが単変量正規分布に従う時、母集団分散の事前分布は逆ガンマ分布に従うことは以前の章で解説しました。一方で、多変量正規分布の共分散行列の事前分布は以下のような逆ウィシャート分布
ただし、
となります。
すなわち、
まとめると、ギブスサンプリングの手順は以下のとおりです。
-
\theta^{(s+1)}
a.y_1 \dots y_n \Sigma^{(s)} \mu_n \Lambda_n
b.\theta^{(s+1)} \sim multivariate\ normal(\mu_n, \Lambda_n) -
\Sigma^{(s+1)}
a.y_1 \dots y_n \theta^{(s+1)} S_n
b.\Sigma^{(s+1)} \sim inverse-Wishart(\nu_0+n, S_n^{-1})
図7.2で利用されているデータは22人の子供たちを対象に、指導前と指導後に行われた二つの読解試験の点数から構成されています。各子供の得点は、多変量正規分布からの独立同一分布に従う標本としてモデル化されています。両試験とも平均点は約50点と想定され、これに基づいて事前分布の期待値は(50, 50)に設定されています。事前分布の分散は625とし、相関は0.5と仮定しています。この相関から、共分散は312.5と計算されます。
図7.2は二つの部分から構成されています：
左図は
右図は新たな生徒の得点に対する事後予測分布を表しています。同様に、等高線はHPD領域を示しています。また、元のデータ点も散布図として表示されています。
import numpy as np
import scipy.stats as stats
import pandas as pd
# CSVファイルを読み込む
reading_df = pd.read_csv('./csv/reading.csv')
# 多変量正規分布ベクトル
def rmvnorm(n, mu, Sigma):
return np.random.multivariate_normal(mu, Sigma, n)
# ウィシャート分布
def rwish(n, df, scale):
return stats.wishart.rvs(df=df, scale=scale, size=n)
# データの読み込み
Y = reading_df[['pretest', 'posttest']].values
# 初期パラメータ設定
mu0 = np.array([50, 50])
L0 = np.array([[625, 312.5], [312.5, 625]])
nu0 = 4
S0 = np.array([[625, 312.5], [312.5, 625]])
n = Y.shape[0]
ybar = Y.mean(axis=0)
Sigma = np.cov(Y, rowvar=False)
THETA, SIGMA, YS = [], [], []
np.random.seed(15)
for s in range(5000):
# Thetaの更新
Ln = np.linalg.inv(np.linalg.inv(L0) + n * np.linalg.inv(Sigma))
mun = Ln @ (np.linalg.inv(L0) @ mu0 + n * np.linalg.inv(Sigma) @ ybar)
theta = rmvnorm(1, mun, Ln).flatten()
# Sigmaの更新
Sn = S0 + (Y - theta).T @ (Y - theta)
Sigma = np.linalg.inv(rwish(1, nu0 + n, np.linalg.inv(Sn)))
# YSの更新
YS.append(rmvnorm(1, theta, Sigma)[0])
# 結果の保存
THETA.append(theta)
SIGMA.append(Sigma.flatten())
THETA = np.array(THETA)
SIGMA = np.array(SIGMA)
YS = np.array(YS)
二つ目の部分では以下のことを計算しています。
-
\theta_2 - \theta_1 -
\theta_2 - \theta_1 p(\theta_2 - \theta_1 |y_1 \dots y_n) p(Y_2 > Y_1|y_1 \dots y_n)
# 統計量の計算
theta_diff_quantiles = np.quantile(THETA[:, 1] - THETA[:, 0], [0.025, 0.5, 0.975])
prob_theta2_greater_theta1 = np.mean(THETA[:, 1] > THETA[:, 0])
prob_YS2_greater_YS1 = np.mean(YS[:, 1] > YS[:, 0])
# 統計量の出力
print(f"Theta difference quantiles: {theta_diff_quantiles}")
# [ 1.49778945 6.57558615 11.73849073]
print(f"P(theta_2 > theta_1): {prob_theta2_greater_theta1}")
# 0.9914
print(f"P(y_2 > y_1): {prob_YS2_greater_YS1}")
# 0.7062
三つ目の部分は図7.2をプロットするプログラムです。
def compute_hpd_2d(x, y, prob_levels):
xy = np.vstack([x, y])
z = gaussian_kde(xy)(xy)
sorted_indices = np.argsort(z)
sorted_z = z[sorted_indices]
cumsum_z = np.cumsum(sorted_z)
cumsum_z /= cumsum_z[-1]
hpd_indices = [np.searchsorted(cumsum_z, p) for p in prob_levels]
hpd_values = sorted_z[hpd_indices]
return hpd_values
def plot_hdr2d(data, ax, xlab, ylab, xlim=None, ylim=None):
x, y = data.T
# カーネル密度推定
kde = gaussian_kde(data.T)
# 格子点
xmin, xmax = x.min(), x.max()
ymin, ymax = y.min(), y.max()
X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([X.ravel(), Y.ravel()])
Z = np.reshape(kde(positions).T, X.shape)
# HPD領域の計算
prob_levels = [0.025, 0.25, 0.5, 0.75, 0.975]
hpd_values = compute_hpd_2d(x, y, prob_levels)
# HPDプロット
cmap = plt.get_cmap('Greys')
ax.contourf(X, Y, Z, levels=hpd_values, cmap=cmap)
ax.contour(X, Y, Z, levels=hpd_values, colors='k', linewidths=0.5)
ax.set_xlim(xlim if xlim else [xmin, xmax])
ax.set_ylim(ylim if ylim else [ymin, ymax])
ax.set_xlabel(xlab)
ax.set_ylabel(ylab)
# プロットの作成
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
# THETAデータのプロット
plot_hdr2d(THETA, axes[0], r'$\theta_1$', r'$\theta_2$')
axes[0].plot([35, 60], [35, 60], 'k--', linewidth=0.5)
# YSデータのプロット
plot_hdr2d(YS, axes[1], r'$y_1$', r'$y_2$', xlim=[0, 100], ylim=[0, 100])
axes[1].scatter(Y[:, 0], Y[:, 1], c='black', s=10, alpha=0.5) # 元のデータ
axes[1].plot([0, 100], [0, 100], 'k--', linewidth=0.5)
plt.tight_layout()
plt.show()
読解力データと事後分布
200人の女性の生体情報に関するデータ
図7.3は、アリゾナ州フェニックス近郊に住むピマ・インディアンを先祖に持つ200人の女性の健康関連情報データセットから取得した4つの変数の単変量のヒストグラムと二変量の散布図を示しています。こちらも特にプログラムの解説は不要だと思いますので、解説は省略します。
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
# CSVファイルを読み込む
pima_df = pd.read_csv('./csv/pima.miss.csv')
# データの選択
Y0 = pima_df.iloc[:, 1:5]
# データの欠損値の導入
np.random.seed(1)
O = np.random.binomial(1, 0.9, Y0.shape)
Y0 = Y0.mask(O == 0)
# 散布図行列とヒストグラムの作成
sns.pairplot(Y0)
plt.show()
200人の女性の生体情報に関するデータ
相関係数と回帰係数に関する95%事後信用区間
図7.4が掲載されている7.5節では欠測データに対するギブスサンプリングが解説されています。
欠測しているデータを以下のように表現するとします。
ただし、
この時、欠測データに対するギブスサンプリングは以下のようになります。
-
\theta^{(s+1)} p(\theta|Y_{obs}, Y_{miss}^{(s)}, \Sigma^{(s)}) -
\Sigma^{(s+1)} p(\Sigma|Y_{obs}, Y_{miss}^{(s)}, \theta^{(s+1)}) -
Y_{miss}^{(s+1)} p(Y_{miss}|Y_{obs},\theta^{(s+1)}, \Sigma^{(s+1)})
ただし、
となります。
図7.4では、四つの変数glu, bp, skin, bmiを持つデータに対してギブスサンプリングを行い、その変数の相関係数と回帰係数をプロットしています。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal, wishart
# 相関係数の計算
def compute_correlation(SIGMA):
n_iter, p, _ = SIGMA.shape
COR = np.zeros_like(SIGMA)
for s in range(n_iter):
Sig = SIGMA[s]
COR[s] = Sig / np.sqrt(np.outer(np.diag(Sig), np.diag(Sig)))
return COR
# 回帰係数の計算
def compute_regression(COR):
n_iter, p, _ = COR.shape
REG = np.zeros_like(COR)
for s in range(n_iter):
for j in range(p):
REG[s, j, :j] = COR[s, j, :j] @ np.linalg.inv(COR[s, :j, :j])
REG[s, j, j+1:] = COR[s, j, j+1:] @ np.linalg.inv(COR[s, j+1:, j+1:])
return REG
# 信頼区間の計算
def compute_quantiles(M):
return np.quantile(M, [0.025, 0.5, 0.975], axis=0)
# 相関係数・回帰係数をプロットするための関数
def plot_coefficients(ax, coef, var_names, row_index):
p = coef.shape[1]
y = coef[1, row_index, :]
ci_low = coef[0, row_index, :]
ci_high = coef[2, row_index, :]
x = np.arange(p)
ax.scatter(x, y, color='black', s=20)
ax.vlines(x, ci_low, ci_high, color='black')
ax.axhline(y=0, color='gray', linestyle='--', linewidth=0.5)
ax.axvline(x=row_index, color='gray', linestyle='--', linewidth=0.5)
ax.set_xlim(-0.5, p-0.5)
ax.set_ylim(-1, 1)
ax.set_yticks(np.arange(-1, 1.1, 0.5))
ax.set_xticks(range(p))
ax.set_xticklabels(var_names)
# メインプロット
def plot_figure_7_4(SIGMA):
COR = compute_correlation(SIGMA)
REG = compute_regression(COR)
COR_q = compute_quantiles(COR)
REG_q = compute_quantiles(REG)
var_names = ['glu', 'bp', 'skin', 'bmi']
fig, axes = plt.subplots(4, 2, figsize=(10, 16))
plt.subplots_adjust(hspace=0.4, wspace=0.3)
for i in range(4):
plot_coefficients(axes[i, 0], COR_q, var_names, i)
plot_coefficients(axes[i, 1], REG_q, var_names, i)
axes[i, 0].set_ylabel(var_names[i])
axes[i, 1].set_ylabel(var_names[i])
if i < 3:
axes[i, 0].set_xticklabels([])
axes[i, 1].set_xticklabels([])
axes[0, 0].set_title('Correlation coefficients')
axes[0, 1].set_title('Regression coefficients')
plt.tight_layout()
plt.show()
# データの読み込み
Y = pd.read_csv('./csv/pima.miss.csv', index_col=0)
Y_complete = pd.read_csv('./csv/pima_nomiss.csv', index_col=0)
n, p = Y.shape
O = ~Y.isna()
# 事前分布のパラメータ
mu0 = np.array([120, 64, 26, 26])
sd0 = mu0 / 2
L0 = np.full((p, p), 0.1)
np.fill_diagonal(L0, 1)
L0 *= np.outer(sd0, sd0)
nu0 = p + 2
S0 = L0
# ギブスサンプリング
np.random.seed(1000)
n_iter = 5000
THETA = []
SIGMA = []
Y_MISS = []
# 初期化
Sigma = S0.copy()
Y_full = Y.fillna(Y.mean())
for s in range(n_iter):
# thetaの更新
ybar = Y_full.mean().values
Ln = np.linalg.inv(np.linalg.inv(L0) + n * np.linalg.inv(Sigma))
mun = Ln @ (np.linalg.inv(L0) @ mu0 + n * np.linalg.inv(Sigma) @ ybar)
theta = multivariate_normal.rvs(mean=mun, cov=Ln)
# Sigmaの更新
Sn = S0 + ((Y_full - theta).T @ (Y_full - theta)).values
Sigma = np.linalg.inv(wishart.rvs(df=nu0 + n, scale=np.linalg.inv(Sn)))
# 欠測データの更新
for i in range(n):
b = ~O.iloc[i, :].values
a = O.iloc[i, :].values
if np.any(b):
iSa = np.linalg.inv(Sigma[np.ix_(a, a)])
beta_j = Sigma[np.ix_(b, a)] @ iSa
s2_j = Sigma[np.ix_(b, b)] - beta_j @ Sigma[np.ix_(a, b)]
theta_j = theta[b] + beta_j @ (Y_full.iloc[i, a].values - theta[a])
Y_full.iloc[i, b] = multivariate_normal.rvs(mean=theta_j, cov=s2_j)
THETA.append(theta)
SIGMA.append(Sigma)
Y_MISS.append(Y_full.values[~O.values])
# 結果を numpy 配列に変換
THETA = np.array(THETA)
SIGMA = np.array(SIGMA)
Y_MISS = np.array(Y_MISS)
# プロットの作成
plot_figure_7_4(SIGMA)
相関係数（左）と回帰係数（右）に関する95%事後信用区間
相関係数の計算
def compute_correlation(SIGMA):
n_iter, p, * = SIGMA.shape
COR = np.zeros_like(SIGMA)
for s in range(n_iter):
Sig = SIGMA[s]
COR[s] = Sig / np.sqrt(np.outer(np.diag(Sig), np.diag(Sig)))
return COR
この関数は共分散行列から相関係数を計算します。入力として3次元の共分散行列（SIGMA）を受け取り、出力として3次元の相関行列（COR）を返します。関数は各反復で1つの共分散行列を処理します。次に、相関係数の計算には一般的な公式
回帰係数の計算
def compute_regression(COR):
n_iter, p, * = COR.shape
REG = np.zeros_like(COR)
for s in range(n_iter):
for j in range(p):
REG[s, j, :j] = COR[s, j, :j] @ np.linalg.inv(COR[s, :j, :j])
REG[s, j, j+1:] = COR[s, j, j+1:] @ np.linalg.inv(COR[s, j+1:, j+1:])
return REG
この関数は相関行列から回帰係数を計算します。入力として3次元の相関行列（COR）を受け取り、出力として3次元の回帰係数行列（REG）を返します。
各変数を順番に目的変数として扱います。これにより、すべての変数間の関係性を分析することができます。次に、回帰係数の計算には、相関行列の部分行列を使用する方法
欠測値の真の値と事後期待値
欠測値の真の値と事後期待値をプロットしています。
# 補完値の平均を計算
Y_imputed_mean = np.zeros_like(Y_complete.values)
missing_indices = np.where(~O.values)
for idx in range(len(missing_indices[0])):
i, j = missing_indices[0][idx], missing_indices[1][idx]
Y_imputed_mean[i, j] = np.mean([Y_MISS[k][idx] for k in range(len(Y_MISS))])
# プロットの作成
fig, axes = plt.subplots(2, 2, figsize=(12, 12))
variables = ['glu', 'bp', 'skin', 'bmi']
for i, (ax, var) in enumerate(zip(axes.flatten(), variables)):
true_values = Y_complete.loc[~O.iloc[:, i], var].values
imputed_values = Y_imputed_mean[~O.iloc[:, i], i]
ax.scatter(true_values, imputed_values, alpha=0.6)
min_val = min(min(true_values), min(imputed_values))
max_val = max(max(true_values), max(imputed_values))
ax.plot([min_val, max_val], [min_val, max_val], 'k-', lw=1)
ax.set_xlabel(f'true {var}')
ax.set_ylabel(f'predicted {var}')
# 軸の範囲を調整
ax.set_xlim(min(true_values) * 0.9, max(true_values) * 1.1)
ax.set_ylim(min(imputed_values) * 0.9, max(imputed_values) * 1.1)
plt.tight_layout()
plt.show()
欠測値の真の値と事後期待値
最後に
本ブログでは、標準ベイズ統計学の7章で扱われている多変量正規モデルについて、Pythonを用いて実装し、視覚化を行いました。次回は8章「グループ比較と階層モデリング」の内容をPythonで実装したいと思います。
Discussion