単回帰モデルのベイズ分析
PyStan 3対応
岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版の第11章において、単回帰モデルのPyStanスクリプトについて説明している。これはPyStan 2対応であるので、PyStan 3を利用するときは、少し変更が必要である。以下に変更例を示す。サンプルスクリプトのファイルは圧縮ファイルc11_files.zipにまとめているので、ダウンロード展開すれば利用できる。
なお、WindowsでPyStan 3を利用するときは、WindowsにUbuntuをインストールしてUbuntu上で利用する。WindowsにUbuntuをインストールすることはWSL（Windows Subsystem for Linux）により簡単にできるが、このウェブサイトで説明している。
拙著第11章のリスト11.2は次のリスト１のように変更するとPyStan 3が利用できる。
PyStan 2では、StanModelによる翻訳、samplingによるMCMCサンプリングの実行と進むが、PyStan 3ではbuildによるビルド、sampleによる実行と進む。リスト１では、
sm =
stan.build(f_stan.read(), data = Data)
fit = sm.sample()
となっている。関数buildの第１引数は、Stanスクリプトの文字列であり、パラメータdataにMCMCサンプリングで使用されるデータが設定されている。第１引数の文字列は、上の例ではファイルからの読み込み値f_stan.read()が設定されている。
MCMCサンプリングは関数sampleによって実行される。
MCMCサンプリングの統計は、サンプリング結果fitを次のコード
print(fit)
のように関数printで直接出力すると以下のような出力になる。
<stan.Fit>
Parameters:
b0: ()
b1: ()
sgm: ()
mu: (10,)
Draws: 4000
もっと詳しい統計を得るときは、Arviz InferenceData型に変換してArvizの関数summaryを用いるとよい。例は、このウェブサイトを参照されたい。
リスト１では、fitをPandasのDataFrame型に変換している。コードは次のようである。
fit = fit.to_frame()
事後分布などは、PandasのDataFrame型オブジェクトのデータとして取り出して処理が行われている。
パラメータが配列のときは、Stanでは１から始まる添え字で配列の要素が区別されるが、Pandasでは配列名に「.添え字」が付いた形のキーで区別されることに注意しなければならない。
リスト１における次のコード
for i in range(N):
YsArray.append(fit[f'mu.{i+1}'])
では、「fit[f'mu.{i+1}']」によってパラメータmuのi+1番目の要素が選ばれる。for文ではiは０から始まるが、PandasのDataframeでは、添え字は１から始まることに注意。
リスト１ 拙著第11章リスト11.2をPyStan 3用に変更したもの（L11_2.py）
import stan
# pystan --> stan
import matplotlib.pyplot as plt
import scipy.stats as scst
import numpy as np
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}
#sm = pystan.StanModel(file = 'Regression.stan')
#fit = sm.sampling(data = Data, n_jobs
= 1)
with open('Regression.stan', 'r') as
f_stan:
sm =
stan.build(f_stan.read(), data = Data)
fit = sm.sample()
print(fit)
f = open('fit.txt', 'w')
f.write('fit = \n{}\n'.format(fit))
f.close()
fit = fit.to_frame()
Sampleb0 = fit['b0']
print('len(Sampleb0) = ',
len(Sampleb0))
plt.title('Samples of b0')
plt.xlabel('b0')
plt.ylabel('frequency')
plt.hist(Sampleb0)
plt.show()
Sampleb1 = fit['b1']
nb1 = len(Sampleb1)
print('len(Sampleb1) = ', nb1)
npos = 0
for v in Sampleb1:
if v > 0.0:
npos
+= 1
plt.title(('Samples of b1\n') +
('Prop(b1 > 0) = {0:.3f}'.format(npos/nb1)))
plt.xlabel('b1')
plt.ylabel('frequency')
plt.hist(Sampleb1)
plt.show()
Samplesgm = fit['sgm']
plt.title('Samples of sgm\n')
plt.xlabel('sgm')
plt.ylabel('frequency')
plt.hist(Samplesgm)
plt.show()
r, p = scst.pearsonr(Sampleb0,
Sampleb1)
print('r = ', r)
plt.xlabel('b0')
plt.ylabel('b1')
plt.title('r = {0:.3f}'.format(r))
plt.plot(Sampleb0, Sampleb1, 'b.')
plt.show()
#PredYs = fit['mu']
YsArray = []
for i in range(N):
YsArray.append(fit[f'mu.{i+1}'])
for i in range(10):
YsArray.append([])
#for v in PredYs:
# for i in range(10):
#
YsArray[i].append(v[i])
MedPreds = []
for i in range(10):
MedPreds.append(np.percentile(YsArray[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.show()
リスト１を実行すると、以下のグラフが順番に表示される。
