PyStan 3による単回帰分析
拙著＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の第11章で扱ったPyStan 2による単回帰分析StanスクリプトをPyStan 3で実行した。Stanスクリプトは変更なしで用いることが出来たが、PythonスクリプトはPyStan 3用に書き換えた。書き換えたスクリプトを、リスト１に、Stanスクリプト（拙著、リスト11.1）をリスト２に示す。スクリプトファイルは、圧縮ファイルfiles_reg.zipにまとめた。
リスト１では、PyStan 3によるビルドにおいてStanスクリプトがファイルから読み込んだ文字列で与えられている。すなわち
with open('Regression.stan', 'r') as
f:
sm = stan.build(f.read(),
data = Data)
である。PyStan 3では、Stanスクリプトは文字列として関数stan.build()の第１パラメータとして設定する。また、データも関数stan.build()のパラメータとして設定する。
なお、PyStan 3のインストールは、AnacondaのPython 3.10仮想環境では
pip install
pystan
であるが、Pythonスクリプトにおけるimportでは
import stan
である。
リスト１ PyStan 3によるMCMCサンプリングを行うスクリプト（Regression_PS3.py）
import stan
import pickle
import numpy as np
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])
for v in zip(ID, Y, X):
print(v)
Data = {'N': N, 'Y': Y, 'X': X}
with open('Regression.stan', 'r') as
f:
sm = stan.build(f.read(),
data = Data)
fit = sm.sample()
print('fit =\n', fit)
infrnc_d = az.from_pystan(posterior =
fit, posterior_model = sm )
print('\nSummary =\n',
az.summary(infrnc_d))
with open('infrnc_d.pkl', 'wb') as f:
pickle.dump(infrnc_d, f)
d_frame = fit.to_frame()
print(d_frame.keys())
with open('d_frame.pkl', 'wb') as f:
pickle.dump(d_frame, f)
リスト１では、次のコード
with open('Regression.stan', 'r') as
f:
sm = stan.build(f.read(),
data = Data)
によりStanスクリプトのビルドを行い、MCMCサンプリングは次のコードで行っている。
fit = sm.sample()
MCMCサンプリングの結果は、コード
infrnc_d = az.from_pystan(posterior =
fit, posterior_model = sm )
print('\nSummary =\n',
az.summary(infrnc_d))
with open('infrnc_d.pkl', 'wb') as f:
pickle.dump(infrnc_d, f)
により、FitオブジェクトをInferenceDataオブジェクトにまとめ、pickleでファイルに保存して、次のリスト３で用いる。
また、InferenceDataオブジェクトに纏められたものから関数az.summary()によってMCMCサンプリングの統計量を求めている。
次のコードでは、MCMCサンプリングの結果をPandas DataFrameオブジェクトにまとめて、pickleで保存している。Stanスクリプトにおけるパラメータ名とDataFrameでのキーとの対応は、print文の出力で確認できる。
d_frame = fit.to_frame()
print(d_frame.keys())
with open('d_frame.pkl', 'wb') as f:
pickle.dump(d_frame, f)
リスト２ Stanスクリプト（Regression.stan）
data{
int N;
real Y[N];
real X[N];
}
parameters{
real b0;
real b1;
real<lower = 0> sgm;
}
transformed parameters{
real mu[N];
for (i in 1:N)
mu[i] = b0 + b1 * X[i];
}
model{
for (i in 1:N)
Y[i]
~ normal(mu[i], sgm);
}
リスト１の実行後、リスト３を実行する。ただし、Windowsの場合は、リスト１は直接Windows上では実行できない。WindowsにWSL（Windows Subsystem for Linux）を用意して、その上にインストールされたUbuntu上でPyStan3を用いる。詳しくはこのウェブサイトを参照されたい。PyStan 3は、直接Windows上で実行できないので、リスト１とリスト３に分けたが、Windows以外では分ける必要はない。
リスト３ リスト１のスクリプトに続いて実行するスクリプト（DrawRegression_PS3.py）
import pickle
import matplotlib.pyplot as plt
import scipy.stats as scst
import numpy as np
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])
for v in zip(ID, Y, X):
print(v)
with open('infrnc_d.pkl', 'rb') as f:
infrnc_d = pickle.load(f)
print('Summary...\n',
az.summary(infrnc_d))
with open('d_frame.pkl', 'rb') as f:
d_frame = pickle.load(f)
print('Keys...\n', d_frame.keys())
plt.hist(d_frame['b0'])
plt.title('MCMC sample of b0',
fontsize = 20)
plt.xlabel('b0', fontsize = 16)
plt.savefig('Figb0.png')
plt.show()
plt.hist(d_frame['b1'])
plt.title('MCMC sample of b1',
fontsize = 20)
plt.xlabel('b1', fontsize = 16)
plt.savefig('Figb1.png')
plt.show()
plt.hist(d_frame['sgm'])
plt.title('MCMC sample of sgm',
fontsize = 20)
plt.xlabel('sgm', fontsize = 16)
plt.savefig('Figsgm.png')
plt.show()
r = scst.pearsonr(d_frame['b0'],
d_frame['b1'])[0]
print('r = ', r)
plt.xlabel('b0')
plt.ylabel('b1')
plt.title('r = {0:.3f}'.format(r),
fontsize = 24)
plt.plot(d_frame['b0'], d_frame['b1'],
'b.')
plt.xlabel('b0', fontsize = 16)
plt.ylabel('b1', fontsize = 16)
plt.savefig('Figb0b1.png')
plt.show()
k_names = [f'mu.{i}' for i in
range(1,11)]
print(k_names)
MedPreds = []
for i in range(10):
MedPreds.append(np.percentile(d_frame[k_names[i]], 50))
plt.xlabel('Data')
plt.ylabel('Prediction')
plt.plot(Y, MedPreds, 'bo')
for nm, y, pred in zip(ID, Y,
MedPreds):
plt.text(y + 0.1, pred +
0.1, nm)
plt.plot([14, 28], [14, 28], 'b-')
plt.savefig('FigPrdctn.png')
plt.show()
リスト３は、拙著のリスト11.2の後半に対応する処理を行うものである。リスト３を実行すると、まず図１のグラフが表示される。
図１
これは、次のコード
with open('d_frame.pkl', 'rb') as f:
d_frame = pickle.load(f)
で、DataFrameのデータを読み込んで、次のコードで描画したものである。
plt.hist(d_frame['b0'])
plt.title('MCMC sample of b0',
fontsize = 20)
plt.xlabel('b0', fontsize = 16)
plt.savefig('Figb0.png')
plt.show()
図１のフォームを閉じると、図２のグラフが表示される。
図２
図２のフォームを閉じると図３のグラフが表示される。
図３
図３のフォームを閉じると、図４のグラフが表示される。
図４
図４は、以下のコードで描画されている。
r = scst.pearsonr(d_frame['b0'],
d_frame['b1'])[0]
print('r = ', r)
plt.xlabel('b0')
plt.ylabel('b1')
plt.title('r = {0:.3f}'.format(r),
fontsize = 24)
plt.plot(d_frame['b0'], d_frame['b1'],
'b.')
plt.xlabel('b0', fontsize = 16)
plt.ylabel('b1', fontsize = 16)
plt.savefig('Figb0b1.png')
plt.show()
図４のフォームを閉じると、図５のグラフが表示される。
図５
横軸がデータ値、縦軸が予測値である。
図５のフォームを閉じると、スクリプトの実行終了である。