単回帰モデルによるベイズ分析
拙著「いまさら聞けないPythonでデータ分析」では、PyStanを用いたStanによるベイズ分析を説明した。ここでは、CmdStanPyによるStanの利用を試みる。CmdStanPyの簡単なインストールの説明のウェブサイトも用意した。
単回帰モデルによる分析を行うデータとして、次のものを用意した。
RawData = [ ['City', 'Temperature',
'Latitude'],
['Sapporo',
15,
43.1],
['Sendai',
20,
38.3],
['Niigata',
17,
37.9],
['Kanazawa',
18,
36.6],
['Tokyo', 24,
35.7],
['Osaka', 22,
34.7],
['Fukuoka',
22,
33.6],
['Kochi', 24,
33.6],
['Kagoshima', 23,
31.6],
['Naha',
26,
26.2]
予想最高気温をhttps://tenki.jp/から、緯度をhttps://www.geocoding.jp/から読み取ったものである。
Latitudeを独立変数XとしてTemperatureを従属変数Yとする単回帰モデル
Y=b0+b1X+residual
の係数をMCMCサンプリングで求めるStanスクリプトとして、リスト１に示すもの（regression.stan）を作成した。
係数の安定化のために独立変数は中心化したもの（平均値が０になるように調整したもの）を用いる。なお、ファイルはregres_files.zipにまとめた。
リスト１ 単回帰モデルのStanスクリプト（regression.stan）
data{
int N;
array[N] real Y;
array[N] real X; // Centered values
real X_min;
real X_max;
real Y_min;
real Y_max;
real R;
}
parameters{
/*
transformed parametersブロックでスケーリングを行う
*/
real<lower=0, upper=1>
vb0;
real vb1;
real<lower=0, upper=1>
vsgm;
}
transformed parameters{
real b0;
real b1;
real sgm;
array[N] real mu;
b0 = Y_min + (Y_max -
Y_min)*vb0;
b1 = (Y_max - Y_min)/(X_max
- X_min) * (R + vb1);
sgm = 0.0001 + vsgm * (Y_max
- Y_min);
for (i in 1:N)
mu[i] = b0 + b1 * X[i]; // 回帰モデル
}
model{
vb0 ~ beta(1,1);
vb1 ~ normal(0.0, 1.0);
vsgm ~ beta(1,1);
Y ~ normal(mu, sgm);
}
リスト１のStanスクリプトでは、modelブロック部でのパラメータの乱数の生成は以下のようになっている。
vb0 ~ beta(1,1);
vb1 ~ normal(0.0, 1.0);
vsgm ~ beta(1,1);
区間
における一様乱数、あるいは平均０，標準偏差１の正規乱数というように、変域が制限されている。これらを基にデータに合わせた変域にtransformed parametersブロックにおいて以下のように調整している。
real b0;
real b1;
real sgm;
array[N] real mu;
b0 = Y_min + (Y_max -
Y_min)*vb0;
b1 = (Y_max - Y_min)/(X_max
- X_min) * (R + vb1);
sgm = 0.0001 + vsgm * (Y_max
- Y_min);
リスト１のStanスクリプトを用いてベイズ分析を行うPythonスクリプトをリスト２（regression.py）のように作成した。
リスト２ リスト１のStanスクリプトを用いる単回帰分析のPythonスクリプト（regression.py）
from cmdstanpy import CmdStanModel
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import arviz as az
RawData = [ ['City', 'Temperature',
'Latitude'],
['Sapporo',
15,
43.1],
['Sendai',
20,
38.3],
['Niigata',
17,
37.9],
['Kanazawa',
18,
36.6],
['Tokyo', 24,
35.7],
['Osaka', 22,
34.7],
['Fukuoka',
22,
33.6],
['Kochi', 24,
33.6],
['Kagoshima', 23,
31.6],
['Naha',
26,
26.2]
]
ID = []
Y = []
X = []
N = len(RawData) - 1
for i in range(N):
ID.append(RawData[i + 1][0])
Y.append(RawData[i + 1][1])
X.append(RawData[i + 1][2])
X_c
= X - X_mean
# Centering
Data = {'N': N, 'Y': Y, 'X': X_c,
'X_min': np.min(X_c), 'X_max': np.max(X_c),
'Y_min': np.min(Y), 'Y_max': np.max(Y),
'R':
ss.pearsonr(X_c, Y).statistic}
model =
CmdStanModel(stan_file="regression.stan") # コンパイル
fit = model.sample(data=Data)
# サンプリング
print(fit.summary())
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
# トレースプロット
az.plot_trace(inf_data,
var_names=['b0', 'b1', 'sgm'], figsize=(10,7))
plt.tight_layout()
plt.show()
df_sample = fit.draws_pd()
"""
(b0, b1)の散布図
"""
b0s = df_sample['b0']
b1s = df_sample['b1']
for s, t in zip(b0s, b1s):
plt.plot(s, t, 'o',
fillstyle='none')
plt.xlabel('b0', fontsize=14)
plt.ylabel('b1', fontsize=14)
plt.show()
"""
データ（X, Y）の散布図と回帰直線
"""
b0_med = np.median(b0s)
b1_med = np.median(b1s)
X_mean = np.mean(X)
min_X = np.min(X)
max_X = np.max(X)
for nm, xv, yv in zip(ID, X, Y):
plt.plot(xv, yv, 'o', c='b')
plt.text(xv+0.1, yv+0.1, nm,
rotation=45)
[b0_med + b1_med * (min_X-X_mean), b0_med + b1_med * (max_X-X_mean)])
plt.xlabel('Latitude', fontsize=14)
plt.ylabel('Temperature', fontsize=14)
plt.title(f'b0={b0_med:.2f}, b1={b1_med:.2f}', fontsize=16)
plt.show()
リスト１のStanスクリプトのコンパイルを、次のコードで行っている。
model =
CmdStanModel(stan_file="regression.stan") # コンパイル
コンパイルしたモデルによるMCMCサンプリングは次のコードで行っている。
fit = model.sample(data=Data)
分析用のデータは、回帰係数とｙ切片の相関を防ぐために中心化（平均値が０になるように調整）を次のコードで行っている。
X_mean = np.mean(X)
X_c = X - X_mean # Centering
中心化を行っているので、ｙ切片の値は独立変数の平均値における回帰直線の値である。
リスト１のファイルとリスト２のファイルを同じフォルダに置いておいて、そのフォルダをカレントディレクトリとしてリスト２のPythonスクリプトを実行する。
例えば、次のコマンドを実行する。
(ckstan)
*****/regres_files$ python regression.py
上のコマンドは、Ubuntu24/wsl/windowsで実行した。
スクリプトが実行されると、コンパイルに多少の時間が掛かるが、その後、サンプリングが行われて、トレース図が表示される（図１）。
図１
Windowの右上のX印アイコンをクリックしてWindowを閉じると、次にｙ切片と傾きのパラメータのサンプルの散布図が表示される（図２）。
の
図２
独立変数を中心化して分析したので、強い相関は認められない。
図２のWindowを閉じると、回帰直線のグラフが表示される（図３）。
図３
ｙ切片と傾きは、事後分布の中央値を点推定値としている。
回帰直線は、パラメータ値が中心化されたデータに基づいていることから、次式で求めている。
plt.plot([min_X, max_X],
[b0_med + b1_med * (min_X-X_mean), b0_med + b1_med * (max_X-X_mean)])
図３のWindowを閉じると、スクリプトの実行終了である。