重回帰モデルのベイズ分析
拙著「Pythonでデータ分析」の「第５章 重回帰分析」のベイズ分析用スクリプトを用意した。
重回帰分析を、変数をすべて標準化して行い、事前分布として、Gelman et al. (2021), p. 124を参考にして、係数の事前分布は平均が０，標準偏差が2.5の正規分布、誤差の事前分布は平均が１の指数分布を採用した（soft constraint; Gelmanら、2021、p.124）。Stanスクリプトをリスト１に示す。スクリプトファイルなどは圧縮ファイルmlregfiles.zipにまとめた。ダウンロード・展開して自由に利用されたい。
リスト１ 重回帰モデルのベイズ分析Stanスクリプト（標準化変数：MLReg.stan）
data {
int N;
int M;
vector[N] y;
array[N] vector[M] X;
}
parameters {
real b0;
vector[M] b;
real<lower = 0.0> sgm;
}
model {
b ~ normal(0.0, 2.5);
sgm
~ exponential(1/1.0);
for (i
in 1:N) {
y[i] ~ normal(b0 + dot_product(b,X[i]), sgm);
}
}
事前分布
b0 ~ normal(0.0, 2.5);
b ~ normal(0.0, 2.5);
sgm
~ exponential(1/1.0);
は、Gelman et al. (2021), p. 124を参考にした。
リスト１のStanスクリプト（CmdStanPy対応）を用いるベイズ分析Pythonスクリプトをリスト２に示す。
リスト１のStanスクリプトのファイル（MLReg.stan）、リスト２のPythonスクリプトのファイル（MLRegr.py）、および入力データファイル（予想最高気温.xlsx）を同じフォルダに置く。そのフォルダに、カレントディレクトリ（フォルダ）を移し、CmdStanPyのインストールされた環境において、次のコマンドを実行する。CmdStanPyの簡単な解説を、このウェブサイトに用意した。
(stan) *****/mlregfiles$ python MLRegr.py
うえのコマンドを実行すると、入力データファイル名の設定が求められr。
(stan) *****/mlregfiles$
python MLRegr.py
入力データファイル(*.xlsx) = 予想最高気温.xlsx
上の例では、入力データファイル名として図１のファイル名を設定している。入力データファイルは、図１の形式のExcelファイルである。
図１
1行目に変数名を設定して、データ値は２行目以降に書く。１列目は、データの識別文字列、２列目に従属変数、３列目以降に独立変数を置く。変数値は、スクリプトの実行時に標準得点化される。
入力データファイル名を設定すると、Stanスクリプトのコンパイルが始まる。コンパイルには多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、MCMCサンプリングが終了すると、図２のトレース図が表示される。
図２
左側に事後分布のグラフ、右側にMCMCサンプリングのトレースが表示されている。
図２のWindowを閉じると、b0の事後分布が表示される（図３）。
図３
変数が標準化されているので、b0の事後分布は原点０を中心としたものとなる。事後分布の代表値として中央値（med）、標準偏差に相当するものとしてmad(median absolute deviation) sdが表示されている。Gelmanら（2021）は、平均値より安定している中央値とそれに基づく統計量を薦めている（ｐ.73）。
図３のWindowを閉じると、b1の事後分布のグラフが表示される（図４）。
図４
標準化回帰係数（standardized regression coefficient）は、標準化独立変数の値が１異なるときの回帰モデルの値の違いを表す（Gelmanら、2021、p.133；岡本、2021、p.55）。
図４のWindowを閉じると、b2の事後分布のグラフが表示される（図５）。
図５
図５のWindowを閉じると、図６の散布図が表示される。
図６
図６のWindowを閉じると、スクリプトの実行終了である。端末には次のメッセージが表示される。
Results.txt
was saved.
出力ファイルResults.txtを開くと、以下のようになっている。
入力データファイル = 予想最高気温.xlsx
y_mean = 20.750, y_sd =
6.193
b1(緯度):
mean = 37.458, sd = 5.784
b2(西経):
mean = 95.342, sd = 17.235
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
b0 0.001 0.210 -0.387 0.392 0.003 0.004 3715.0 2573.0 1.0
b[0] -0.808 0.230 -1.235 -0.366 0.004 0.005 3539.0 2214.0 1.0
b[1] 0.267 0.231 -0.152 0.721 0.004 0.005 3252.0 2267.0 1.0
sgm 0.745 0.200 0.438 1.118 0.004 0.005 2244.0 2393.0 1.0
b0_med = -0.002, b0_mad_sd = 0.200
b1(緯度):
med = -0.808, mad_sd = 0.215
b2(西経):
med = 0.268, mad_sd = 0.214
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.
岡本安晴（2022）．データを使いこなすための統計入門．22世紀アート（電子書籍）
リスト２ ベイズ分析Pythonスクリプト（MLRegr.py）
import pandas as pd
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
import japanize_matplotlib
import seaborn as sb
from cmdstanpy
import CmdStanModel
import arviz
as az
"""
Excelファイルからの読み込み
"""
in_file_name = input('入力データファイル(*.xlsx) = ')
xlsx = pd.ExcelFile(in_file_name)
data_xlsx = pd.read_excel(xlsx)
"""
出力ファイルの設定
"""
f = open('Results.txt', 'w')
f.write('\n入力データファイル = ' + in_file_name)
keys = data_xlsx.keys()
# １行目の変数名の読み込み
print('Keys...\n', keys)
case_name = keys[0]
y_name = keys[1]
x_names = keys[2:]
print(y_name)
print(x_names)
Data = data_xlsx.values
print(Data)
case_id = Data[:,0]
y = Data[:,1]
X = Data[:,2:]
print(case_id)
print(y)
print(X)
y_mean = np.mean(y)
y_sd = np.std(y)
f.write(f'\n\ny_mean
= {y_mean:.3f}, y_sd = {y_sd:.3f}\n')
X = np.array(X,
dtype='float')
x_means = np.mean(X,
axis=0)
print(x_means)
x_sds = np.std(X,
axis=0)
print(x_sds)
for j in range(len(x_names)):
f.write(f'\nb{j+1}({x_names[j]}): ')
f.write(f'mean = {x_means[j]:.3f}, sd = {x_sds[j]:.3f}\n')
y = ss.zscore(np.array(y, dtype ='float'))
X = ss.zscore(np.array(X, dtype = 'float'))
print(y)
print(X)
N, M = np.shape(X)
print(N, M)
model = CmdStanModel(stan_file='MLReg.stan')
fit = model.sample(data={'N':N,
'M':M, 'y':y, 'X':X})
print(fit.summary())
d_frm = fit.draws_pd()
# Pandas DataFram型への変換
i_data = az.from_cmdstanpy(fit) # Arviz InferenceData型への変換
print(az.summary(i_data))
f.write('\n\n')
f.write(az.summary(i_data).__str__())
az.plot_trace(i_data)
plt.tight_layout()
plt.show()
def med_mad_sd(x):
""" 中央値とmad sdの計算
Gelman et al. (2021), p.73
"""
med = np.median(x)
mad = np.median(np.abs(x - med))
mad_sd
= 1.483 * mad
return med, mad_sd
b0_med, b0_mad_sd = med_mad_sd(d_frm['b0'])
print(f'\nb0_med = {b0_med:.3f}, b0_mad_sd =
{b0_mad_sd:.3f}')
f.write(f'\n\nb0_med = {b0_med:.3f}, b0_mad_sd =
{b0_mad_sd:.3f}\n')
sb.kdeplot(d_frm['b0'])
plt.title('Posterior distribution for b0' +
f'\nMed. = {b0_med:.3f}, mad sd =
{b0_mad_sd:.3f}')
plt.show()
b_med = np.empty(M)
b_mad_sd = np.empty(M)
for j in range(M):
b_med[j],
b_mad_sd[j] = med_mad_sd(d_frm[f'b[{j+1}]'])
print(x_names[j])
print(b_med[j],
b_mad_sd[j])
f.write(f'\nb{j+1}({x_names[j]}):' +
f' med = {b_med[j]:.3f}, mad_sd = {b_mad_sd[j]:.3f}\n')
sb.kdeplot(d_frm[f'b[{j+1}]'])
plt.title(f'Posterir distribution for b{j+1}: ' + x_names[j]
+
f'\nMed. = {b_med[j]:.3f}, mad sd =
{b_mad_sd[j]:.3f}')
plt.show()
y_hat = np.empty(N)
for i in
range(N):
y_hat[i] = b0_med + np.dot(b_med, X[i])
for i in
range(N):
plt.plot(y_hat[i], y[i],
'o', c='b')
plt.text(y_hat[i], y[i],
case_id[i])
plt.xlabel(r'$\hat{y}$(model)')
plt.ylabel(r'$y$(data)')
x_min = np.min(y_hat)
x_max = np.max(y_hat)
plt.plot([x_min, x_max], [x_min, x_max], label =r'y(data) =
$\hat{y}$(model)')
plt.legend()
plt.show()
print('Results.txt was saved.')