単回帰モデルによるベイズ分析
StanのPython用モジュールPyStan 3は、Windows上では直接利用できない。WindowsにWSL（Windows Subsystem for Linux）をインストールして用いることができる。WindowsでのPyStan 3の準備については、このウェブサイトなどで説明している。
WindowsでPyStan 3を利用するときは、PyStan 3によるMCMCサンプリングをUbuntu/WSL/Windowsで行うが、WSLでは描画ができなかったので、描画はWindows上で行った。
単回帰モデルによるベイズ分析のスクリプトをリスト１とリスト２に分けて示す。ファイルは、まとめて圧縮ファイルregresfiles.zipとしている。
リスト１のスクリプトは、Miniconda/Ubuntu24/WSL/Windows10に用意したPython 3.12仮想環境において実行した。リスト１のスクリプトの実行によってファイルに保存されたデータは、リスト２のスクリプトで読み出されてグラフ描画が行われる。リスト１のスクリプトファイルとリスト２のスクリプトファイルは同じフォルダに置いておく。
リスト２は、Windows上でPython 3.13.2により実行した。
リスト１におけるデータ
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/から読み取ったものである。
import stan
import pickle
import numpy as np
import scipy.stats as ss
import arviz as az
StanCode = """
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);
}
"""
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,
'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}
sm = stan.build(StanCode, data=Data)
fit = sm.sample()
print(fit)
inf_data =
az.from_pystan(posterior=fit, posterior_model=sm)
print(az.summary(inf_data))
df_sample = fit.to_frame()
"""
データとMCMCをファイルに保存して、次の処理で使用する。
"""
with open('temp.pkl', 'wb') as f:
pickle.dump((ID, Y, X,
inf_data, df_sample), f)
Latitudeを独立変数XとしてTemperatureを従属変数Yとする単回帰モデル
Y=b0+b1X+residual
の係数をMCMCサンプリングで求めるStanスクリプトである。
係数の安定化のために独立変数を中心化している。
X_mean = np.mean(X)
X_c = X - X_mean # Centering
また、PyStan 3では、パラメータの生成分布の分散は大きい場合に警告が出るので、modelブロックでの生成では分散は大きくならないようにして、transformed parametersブロックにおいて必要な大きさにスケーリングしている。
リスト１によるMCMCサンプルの結果とデータは、次のリスト２のスクリプトによる分析で用いるために、次のコードでファイルに保存している。
with open('temp.pkl', 'wb') as f:
pickle.dump((ID, Y, X,
inf_data, df_sample), f)
リスト１のスクリプトは、Miniconda/Ubuntu24/WSL/Windows10に用意したPython 3.12仮想環境において実行した。WSLではグラフ描画ができなかったので、グラフ描画のためのスクリプトをリスト２として作成して、直接Windows上でPython 3.13.2により実行した。リスト１のスクリプトファイルとリスト２のスクリプトファイルは、同じフォルダ内に置かれているものとする。
import pickle
import matplotlib.pyplot as plt
import scipy.stats as scst
import numpy as np
import arviz as az
import seaborn as sb
"""
リスト１の実行において保存されたファイルの読み込み
"""
with open('temp.pkl', 'rb') as f:
(ID, Y, X, inf_data,
df_sample) = pickle.load(f)
for v in zip(ID, Y, X):
print(v)
"""
トレース図の描画
"""
az.plot_trace(inf_data,
var_names=['b0', 'b1', 'sgm'], figsize=(10,7))
plt.tight_layout()
plt.show()
"""
パラメータb0の事後分布
"""
sb.kdeplot(df_sample['b0'])
plt.show()
"""
パラメータb1の事後分布
"""
sb.kdeplot(df_sample['b1'])
plt.show()
"""
パラメータsgmの事後分布
"""
sb.kdeplot(df_sample['sgm'])
plt.show()
"""
(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()
リスト２のスクリプトを実行すると、次のコードによりリスト１のスクリプトの実行時に保存されたファイルが読み込まれる。
with open('temp.pkl', 'rb') as f:
(ID, Y, X, inf_data,
df_sample) = pickle.load(f)
ファイルを読み込むと、まず次のコードによりMCMCサンプリングのトレース図が表示される（図１）。
az.plot_trace(inf_data,
var_names=['b0', 'b1', 'sgm'], figsize=(10,7))
plt.tight_layout()
plt.show()
図１
図１のWindowを閉じると、パラメータb0の事後分布が表示される（図２）。
図２
図２のWindowを閉じると、パラメータb1の事後分布が表示される（図３）。
図３
図３のWindowを閉じると、パラメータsgmの事後分布が表示される（図４）。
図４
図４のWindowを閉じると、パラメータb0とb1の散布図が表示される（図５）。
図５
独立変数の中心化を行って回帰分析を行っているので、強い相関は見られないことが判る。
図５のWindowを閉じると、図６のグラフが表示される。
図６
データ点を青の小円で表している。緑の直線は、パラメータの事後分布の中央値を点推定値としたときの回帰直線である。
独立変数を中心化して回帰分析を行っているので、b0は独立変数の平均値における回帰直線の値である。b1は、緯度1度の差に対する回帰直線の値の変化量である。
なお、モジュールpickleについては、拙著 pp.197-200を参照されたい。