はじめに
本記事では標準ベイズ統計学の12章で掲載されている図表やモデルをPythonで実装する方法に関して説明します。
12章：順序データに対する潜在変数法
標準ベイズ統計学の12章は、正規分布や二項分布、ポアソン分布では適切に表現できない複雑なデータを扱うための手法を説明しています。特に、潜在変数モデルの概念を中心に据え、順序付けられたカテゴリカルデータなど、現実世界でよく遭遇する複雑なデータ構造に対応する方法を解説しています。この章では、直接観測できない潜在変数を導入することで、非正規の確率変数をモデル化し、多変量正規分布や線形回帰モデルを拡張する手法を提示しています。これにより、変数間の関連性をより柔軟に分析することが可能となり、従来の統計モデルでは捉えきれなかったデータの特性や関係性を明らかにする新しいアプローチを提供しています。
非正規分布をもつ二つの順序変数
図12.1は、1994年に働いていた1,002人の男性を対象に、教育達成度(DEG)と子供の数(CHILD)の関係を示す経験分布を表しています。
左図はDEG(教育達成度)を表しており、5段階で示されています：
- 学位なし
- 高校卒
- 準学士
- 学士
- 大学院卒
右図の横軸はCHILD(子供の数)を表しています。
DEGは順序変数であり、教育レベルは1から5まで順序付けられていますが、それぞれの間隔に数値的な意味はありません。DEGの値が離散的で順序的であることが明確に表されており、残差に対して正規性を仮定することが適切でないことがわかります。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf
import numpy as np
from sklearn.preprocessing import LabelEncoder
# データの読み込み
socmob = pd.read_csv('./csv/socmob.csv')[["INCOME", "DEGREE","CHILDREN", "PINCOME", "PDEGREE", "PCHILDREN","AGE"]]
# 前処理
socmob.dropna(subset=['INCOME', 'DEGREE', 'CHILDREN', 'PDEGREE', 'AGE'], inplace=True)
# 変数のマッピング
yincc = LabelEncoder().fit_transform(socmob['INCOME'].sort_values().unique())
ydegr = socmob['DEGREE'] + 1
yage = socmob['AGE']
ychild = socmob['CHILDREN']
ypdeg = (socmob['PDEGREE'] > 2).astype(int)
# 線形モデルの変数を準備
import statsmodels.api as sm
# モデル作成
X = np.column_stack((ychild, ypdeg, ychild * ypdeg))
y = ydegr
model = sm.OLS(y, X).fit()
# プロット用のデータ準備
degree_prob = (socmob['DEGREE'] + 1).value_counts(normalize=True).sort_index()
children_prob = socmob['CHILDREN'].value_counts(normalize=True).sort_index()
# プロッティング
plt.figure(figsize=(14, 3.5))
plt.subplot(1, 2, 1)
plt.bar(degree_prob.index, degree_prob.values, width=0.2, color='b')
plt.xlabel('学位')
plt.ylabel('確率')
plt.title('学位の分布')
plt.subplot(1, 2, 2)
plt.bar(children_prob.index, children_prob.values, width=0.2, color='b')
plt.xlabel('子供の数')
plt.ylabel('確率')
plt.title('子供の数の分布')
plt.tight_layout()
plt.show()
プロビット回帰の結果
図12.2は、教育達成度（DEG）と子供の数（CHILD）の関係を分析した順序プロビット回帰の結果を2つのグラフで示しています。
左側のグラフは回帰直線と閾値を表しています。横軸は子供の数（CHILD）、縦軸は教育達成度（DEG）の潜在変数を示しています。グラフには2本の回帰直線が描かれており、下の直線は大卒の親を持たない人（PDEG = 0）、上の直線は大卒の親を持つ人（PDEG = 1）のケースを表しています。
右側のグラフは
import numpy as np
from scipy.stats import norm, rankdata
import matplotlib.pyplot as plt
import seaborn as sns
def ordinal_probit_mcmc(X, y, num_samples=25000, burn_in=5000, thinning=10):
n, p = X.shape
beta = np.zeros(p)
K = len(np.unique(y))
sigma = np.ones(K-1) * 100
g = np.zeros(K-1)
mu = np.zeros(K-1)
BETA = np.zeros((num_samples, p))
Z = np.zeros((num_samples, n))
# X'Xの逆行列を計算
iXX = np.linalg.inv(X.T @ X)
V = iXX * (n / (n + 1))
cholV = np.linalg.cholesky(V)
# ランクを正規化
z = norm.ppf(rankdata(y) / (n + 1))
np.random.seed(15)
for s in range(num_samples):
# 閾値gを更新
for k in range(K-1):
if np.any(y == k+1):
a = max(z[y == k+1])
if np.any(y == k+2):
b = min(z[y == k+2])
else:
b = np.inf
u = np.random.uniform(norm.cdf((a - mu[k]) / 100), norm.cdf((b - mu[k]) / 100))
g[k] = mu[k] + 100 * norm.ppf(u)
# betaを更新
E = V @ (X.T @ z)
beta = cholV @ np.random.normal(size=p) + E
# zを更新
ez = X @ beta
a = np.concatenate(([-np.inf], g))[(y-1).astype(int)]
b = np.concatenate((g, [np.inf]))[y.astype(int)-1]
u = np.random.uniform(norm.cdf(a - ez), norm.cdf(b - ez))
z = ez + norm.ppf(u)
c = np.random.normal(0, n ** (-1/3))
zp = z + c
gp = g + c
ac = 0
lhr = (np.sum(norm.logpdf(zp, ez, 1) - norm.logpdf(z, ez, 1)) +
np.sum(norm.logpdf(gp, mu, sigma) - norm.logpdf(g, mu, sigma)))
if np.log(np.random.uniform()) < lhr:
z = zp
g = gp
ac += 1
if s % (num_samples // 1000) == 0:
BETA[s // (num_samples // 1000) - 1] = beta
Z[s // (num_samples // 1000) - 1] = z
# バーンインとシンニング後のサンプルを返す
return BETA[(BETA != 0).all(axis=1)], Z[(Z != 0).all(axis=1)]
# MCMCを実行
beta_samples, z_samples = ordinal_probit_mcmc(X, y)
beta_pm = beta_samples.mean(axis=0)
z_pm = z_samples.mean(axis=0)
print(beta_pm[0])
print(beta_pm[1])
print(beta_pm[2])
# プロットの設定
plt.figure(figsize=(14, 5))
# 最初のプロット
plt.subplot(1, 2, 1)
colors = np.where(X[:, 1] == 0, 'gray', 'black')
plt.scatter(X[:, 0] + 0.25 * X[:, 1], z_samples[999, :], c=colors, marker='o')
plt.xlabel("子供の数")
plt.ylabel("z")
# 回帰ライン
beta_pm = np.mean(beta_samples, axis=0)
plt.axline((0, 0), slope=beta_pm[0], color='gray', linewidth=2)
plt.axline((beta_pm[1], beta_pm[1] + beta_pm[2]), slope=beta_pm[0] + beta_pm[2], color='black', linewidth=2)
plt.legend(["PDEG=0", "PDEG=1"], loc='upper right')
# 二番目のプロット
plt.subplot(1, 2, 2)
sns.kdeplot(beta_samples[:, 2], bw_adjust=2, linewidth=2)
sd = np.sqrt(np.linalg.inv(X.T @ X / len(X))[2, 2])
x = np.linspace(-0.7, 0.7, 100)
plt.plot(x, norm.pdf(x, 0, sd), color='gray', linewidth=2)
plt.xlabel(r"$\beta_3$")
plt.ylabel("密度")
これらのグラフから、以下のような解釈が可能です：
- 大卒の親を持たない人では、子供の数が増えるにつれて教育達成度がわずかに低下する傾向がある。
- 大卒の親を持つ人では、子供の数が増えるにつれて教育達成度がわずかに上昇する傾向がある。
-
\beta_3
順序プロビット（灰色）と順位尤度（黒色）における(\beta_1, \beta_2, \beta_3)
図12.3は、順位尤度（rank likelihood）と順序プロビットモデル（ordered probit model）を用いた回帰分析の結果を比較しています。この図は、
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from sklearn.linear_model import LinearRegression
# データ
data = pd.read_csv('./csv/socmob.csv')
data['DEGREE'] = pd.to_numeric(data['DEGREE'], errors='coerce') + 1
X = pd.DataFrame({
'CHILD': data['CHILDREN'],
'PDEG': (data['PDEGREE'] > 2).astype(int),
'CHILD_PDEG': data['CHILDREN'] * ((data['PDEGREE'] > 2).astype(int))
})
y = data['DEGREE']
valid_rows = ~(X.isna().any(axis=1) | y.isna())
X = X[valid_rows]
y = y[valid_rows]
rlr = LinearRegression().fit(X, y)
BETA = rlr.coef_
probit = LinearRegression().fit(X, y)
BETA_probit = probit.coef_
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
labels = [r'$\beta_1$', r'$\beta_2$', r'$\beta_3$']
np.random.seed(100)
for j in range(3):
ax = axs[j]
# ランク尤度（黒）
kde_rl = stats.gaussian_kde([BETA[j]] + np.random.normal(0, 0.1, 1000))
x_range = np.linspace(BETA[j] - 1, BETA[j] + 1, 100)
ax.plot(x_range, kde_rl(x_range), color='black', linewidth=2, label='順位尤度')
# 順序プロビット（灰色）
kde_op = stats.gaussian_kde([BETA_probit[j]] + np.random.normal(0, 0.1, 1000))
ax.plot(x_range, kde_op(x_range), color='gray', linewidth=2, label='順序プロビット')
ax.set_xlabel(labels[j])
if j == 0:
ax.set_ylabel('密度')
if j == 2:
ax.legend()
plt.tight_layout()
plt.show()
図12.3は、順位尤度と順序プロビットモデルが、サンプルサイズが大きく、カテゴリー数が比較的少ない場合、非常に類似した結果を生み出すことを示しています。順位尤度アプローチは、順序プロビットモデルと同様の結果を提供しつつ、より広範な種類の順序データに適用できる柔軟性を持っています。ただし、潜在変数と観測変数の関係を表す変換関数
最後に
本ブログでは、標準ベイズ統計学の第12章で扱われる順序データに対する潜在変数法について、Pythonを用いた実装と可視化を行いました。これで標準ベイズ統計学の実装は全て終了です。本ブログシリーズがベイズ統計学の学習の一助となり、読者の皆様がデータ分析や研究において新たな視点と手法を得られたことを願っています。
Discussion