単回帰モデルによるベイズ分析
拙著「いまさら聞けないPythonでデータ分析」では、PyStan2を用いた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+b1*X+residual
の係数をMCMCサンプリングで求めるStanスクリプトとして、リスト１に示すもの（regression.stan）を作成した。
係数の安定化のために独立変数は中心化したもの（平均値が０になるように調整したもの）を用いる。なお、ファイルはregres_files.zipにまとめた。
リスト１ 単回帰モデルのStanスクリプト（regression.stan）
data{
int N;
array[N] real Y;
array[N] real X; // Centered values
}
parameters{
real b0;
real b1;
real<lower=0.0001,
upper=1000.0> sgm;
}
transformed parameters{
array[N] real mu;
for (i
in 1:N)
mu[i] = b0 + b1 * X[i]; // 回帰モデル
}
model{
b0 ~ normal(0.0, 1000.0);
b1 ~ normal(0.0, 1000.0);
sgm
~ uniform(0.0001, 1000.0);
Y ~ normal(mu, sgm);
}
パラメータの事前分布を、
b0 ~ normal(0.0, 1000.0);
b1 ~ normal(0.0, 1000.0);
sgm
~ uniform(0.0001, 1000.0);
と設定している。
Gelman et al. (2021)は、分析の第１段階において、事前分布として一様分布を採用している（p. 123）。Gelman et al. (2014)は、事後分布が無意味にならない程度に実社会の情報を含む弱い事前分布が望ましいとしている（p. 55）。Lambert (2018)は、一様分布のメリットを、「尤度関数によって、データをして語らしめる」ことであるとして、このことは多くの分析において平らな分布が用いられることの正当化とされていると説明している（p. 97）。
リスト１のStanスクリプトを用いるPythonスクリプトをリスト２のように用意した。
リスト２ リスト１の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_mean = np.mean(X)
X_c = X - X_mean # Centering
Data = {'N': N, 'Y': Y, 'X': X_c}
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()
# Pandas DataFrame型に変換
"""
(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)
plt.plot([min_X, max_X],
[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スクリプトをCmdStanPyのインストールされている環境で実行する。CmdStanPyの簡単な説明は、このウェブサイトで行った。
例えば、次のコマンドを実行する。
(stan) *****/regres_files$ python regression.py
上のスクリプトを実行すると、Stanスクリプトがコンパイルされる。コンパイルに多少の時間が掛かるが、その後、サンプリングが行われて、トレース図が表示される（図１）。
図１
Windowの右上のX印アイコンをクリックしてWindowを閉じると、次にｙ切片と傾きのパラメータのサンプルの散布図が表示される（図２）。
の
図２
独立変数を中心化して分析したので、強い相関は認められない。
図２のWindowを閉じると、回帰直線のグラフが表示される（図３）。
図３
ｙ切片と傾きは、事後分布の中央値を点推定値としている。Gelman et al. (2021)は、事後分布の代表値として、平均値より安定した中央値を勧めている。
回帰直線は、パラメータ値が中心化されたデータに基づいていることから、次式で求めている。
plt.plot([min_X, max_X],
[b0_med + b1_med
* (min_X-X_mean), b0_med + b1_med * (max_X-X_mean)])
独立変数の中心化のみ行っているので、傾きb1_medは緯度1度あたりの気温℃の違いである。ｙ切片b0_medは、緯度の平均値における気温である。
図３のWindowを閉じると、スクリプトの実行終了である。
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014) Bayesian data analysis, third edition. CRC Press
Gelman, A., Hill, J., & Vehtari,A. (2021) Regression and other sories. Cambridge University Press.
Lambert, B. (2018) A student’s guide to Bayesian statistics. SAGE
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版