PyStan 2.19.1.1のWindowsでの使用例
PyStan 2.19.1.1をWindowsで使ってみた。PyStan 3は、WSL上で使える。
PyStan 2.19.1.1を使うために、WSL上にPython 3.8の仮想環境をAnaconda/Ubuntu/WSL/Windowsに作成した。Python 3.8仮想環境上でリスト１のスクリプトを実行した。リスト１のスクリプトは拙著の第11章の例をMCMCサンプリングの部分とサンプリングの結果の処理の部分に分け、MCMCサンプリングの部分を書き直したものである。サンプリングの結果の処理はリスト３のスクリプトによって行った。リスト１のスクリプトはPyStan 2.19.1.1と組み合わされるC++との関係でWSL上で実行したが、サンプリングの結果の処理（リスト３）は、Windows上にインストールしたPython 3.11で行った。
なお、スクリプトファイルなどは圧縮ファイルRegressWSL.zipにまとめたので、自由にダウンロードして利用されたい。
リスト１ MCMCサンプリングのPythonスクリプト（リスト２のStanスクリプトRegression.stanを使用）
import pystan
import arviz as az
import pandas as pd
import pickle
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)
print(fit)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
with open('i_data.pkl', 'wb') as f:
pickle.dump(i_data, f)
d_frame = {'b0':fit['b0'], 'b1':
fit['b1'], 'sgm': fit['sgm']}
for i in range(10):
d_frame[f'mu.{i+1}'] =
fit['mu'][:,i]
d_frame = pd.DataFrame(d_frame)
print(d_frame)
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);
}
MCMCサンプリングの結果は、2つのファイルに保存している。
１つは、次のコードによるものである。
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
with open('i_data.pkl', 'wb') as f:
pickle.dump(i_data, f)
上のコードにより、サンプリングの結果は、InferenceDataオブジェクトとして保存され、MCMCサンプリングのトレースの描画に用いる。
他の１つは、次のコードで保存される。
d_frame = {'b0':fit['b0'], 'b1': fit['b1'],
'sgm': fit['sgm']}
for i in range(10):
d_frame[f'mu.{i+1}'] =
fit['mu'][:,i]
d_frame = pd.DataFrame(d_frame)
print(d_frame)
with open('d_frame.pkl', 'wb') as f:
pickle.dump(d_frame, f)
上のコードにより、MCMCサンプリングの結果が、PandasのDataFrameオブジェクトとして保存され、MCMCサンプルの統計量の計算に利用する。
リスト１によるMCMCサンプリングにより保存されたサンプルの分析は、次のリスト３のスクリプトにより行った。リスト３の実行は、WindowsにインストールされたPython 3.11により行った。
import pandas as pd
import arviz as az
import pickle
import matplotlib.pyplot as plt
import scipy.stats as scst
import numpy as np
pd.set_option('display.max_columns',
None)
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('i_data.pkl', 'rb') as f:
i_data = pickle.load(f)
az.plot_trace(i_data)
plt.show()
with open('d_frame.pkl', 'rb') as f:
d_frame = pickle.load(f)
print('d_frame...\n', d_frame)
print('Statistics...')
print(d_frame.describe(percentiles =
[0.025, 0.25, 0.5, 0.75, 0.975]))
Sampleb0 = d_frame['b0']
print('len(Sampleb0) = ',
len(Sampleb0))
plt.title('Samples of b0')
plt.xlabel('b0')
plt.ylabel('frequency')
plt.hist(Sampleb0)
plt.show()
Sampleb1 = d_frame['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 = d_frame['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()
YsArray = []
for i in range(10):
YsArray.append([])
for i in range(10):
for v in
d_frame[f'mu.{i+1}']:
YsArray[i].append(v)
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()
リスト３を実行すると、図１のグラフが表示される。
図１
これは、次のコードにより描かれている。
with open('i_data.pkl', 'rb') as f:
i_data = pickle.load(f)
az.plot_trace(i_data)
plt.show()
図１のフォームを閉じると、図２のグラフが描かれる。
図２
これは、次のコードで描かれている。
with open('d_frame.pkl', 'rb') as f:
d_frame = pickle.load(f)
print('d_frame...\n', d_frame)
print('Statistics...')
print(d_frame.describe(percentiles =
[0.025, 0.25, 0.5, 0.75, 0.975]))
Sampleb0 = d_frame['b0']
print('len(Sampleb0) = ',
len(Sampleb0))
plt.title('Samples of b0')
plt.xlabel('b0')
plt.ylabel('frequency')
plt.hist(Sampleb0)
plt.show()
また、上のコードのprint文により、pandasのDataFrameの関数describe()により算出された統計量が、Python IDLEのコンソール画面に以下のように出力される。
Statistics...
b0
b1
sgm
mu.1
mu.2 \
count 4000.000000 4000.000000 4000.000000 4000.000000 4000.000000
mean 44.928811 -0.678369 2.225930 15.691092 18.947265
std
6.172339
0.173885
0.701283
1.550733
0.906232
min
7.070231
-1.440583
1.023010
9.242012
14.866760
2.5% 32.192090 -1.012165 1.312580 12.657339 17.104291
25% 41.375493 -0.786014 1.735884 14.757127 18.403561
50% 44.950409 -0.679924 2.082504 15.678865 18.946082
75% 48.777279 -0.578110 2.541439 16.563783 19.492584
97.5% 56.836867 -0.322316 3.987685 18.866075 20.767380
max 72.229269 0.419692 8.009790 25.431570 24.964179
mu.3
mu.4
mu.5
mu.6
mu.7 \
count 4000.000000 4000.000000 4000.000000 4000.000000 4000.000000
mean 19.218613 20.100493 20.711026 21.389395 22.135601
std
0.867287
0.771062
0.738418
0.740242
0.787891
min 15.251271 16.389973 17.178305 18.013080 18.586718
2.5% 17.486566 18.560711 19.214789 19.860236 20.524563
25% 18.686603 19.627943 20.264578 20.950071 21.660662
50% 19.215353 20.117354 20.719818 21.394601 22.147361
75% 19.745631 20.555926 21.145212 21.855381 22.631996
97.5% 20.934134 21.642002 22.174498 22.851576 23.651405
max 24.925229 24.798644 24.711008 24.743498 25.997744
mu.8
mu.9
mu.10
count 4000.000000 4000.000000 4000.000000
mean 22.135601 23.492340 27.155534
std
0.787891
0.969857
1.735643
min 18.586718 18.603963 16.854722
2.5% 20.524563 21.529580 23.550024
25% 21.660662 22.900000 26.130414
50% 22.147361 23.512339 27.176282
75% 22.631996 24.104748 28.231530
97.5% 23.651405 25.356349 30.568339
max 25.997744 28.278191 34.486004
len(Sampleb0) = 4000
図２のフォームを閉じると、図３のグラフが表示される。
図３
これは、次のコードによるものである。
Sampleb1 = d_frame['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 = d_frame['sgm']
plt.title('Samples of sgm\n')
plt.xlabel('sgm')
plt.ylabel('frequency')
plt.hist(Samplesgm)
plt.show()
図４のフォームを閉じると、図５のグラフが表示される。
図５
これは、以下のコードによって描かれている。また、b0とb1のサンプルの相関係数も計算され、PythonのIDLEのコンソール画面に表示さる。
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()
図５のフォームを閉じると、図６のデータと予測値の関係を表す散布図が表示される。図中の直線はＸ＝Ｙを表すもので、データと予測値が一致すれば、この直線上に点（データ、予測値）が載ることになる。
図６
これは、以下のコードによるものである。
YsArray = []
for i in range(10):
YsArray.append([])
for i in range(10):
for v in
d_frame[f'mu.{i+1}']:
YsArray[i].append(v)
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()
図６のフォームを閉じると、スクリプトの実行終了である。