拙著「Pythonでデータ分析」の第12章「PyStanによるポアッソン回帰モデル分析」更新
StanスクリプトをPyStan 3対応とするとともに、グラフ描画を加えた。PyStan 3対応のStanスクリプトをリスト１に示し、リスト１のStanスクリプトを使用するPythonスクリプトをリスト２に示した。これらのファイルは圧縮ファイルc12files.zipにまとめたので、ダウンロード・展開して自由に利用されたい。
リスト２のスクリプトをWindowsで実行するときは、WindowsにUnbuntuすなわちWSL（Windows Subsystem for Linux）をインストールして、WindowsのLinux上で実行すればよい。WindowsにUbuntu/WSLをインストールして、Python 3.11/Anacondaにおいてリスト２を実行する方法については、このウェブサイトを参照されたい。
リスト２のスクリプトをコマンド
(py311) ・・・/c12/c12files$ python
PoissonRegStan.py
により実行すると、MCMCサンプリング後、トレース図が表示される（図１）。
図１
図１のフォームを閉じると、図２のグラフが表示される。
図２
主効果、a2、a3、b2、b3の事後分布のグラフである。拙著の表8.2.2において、セル（小学生の時に知った、心理学≒臨床心理学）を基準に主効果などが表されているが（拙著の表12.2）、a3とb3の事後分布が正の領域に位置していることは、「高校生のとき」および「心理学＞臨床心理学」の条件での度数が多いことに対応している。
図２のフォームを閉じると、図３のフォームが表示される。
図３
交互作用g32の事後分布は負の領域に寄っており、g33の分布はほぼ負の領域にある。主効果のみで表される度数に比べて、セル（高校生、その他）および（高校生、心理学＞臨床心理学）の度数が少ないことが表されている。
図３のフォームを閉じると、図４のグラフが表示される。
図４
ポアッソン分布（拙著、p.203）
におけるパラメータ
の推定値（事後分布の中央値）とデータ値
の散布図である。ポアッソン分布の平均値は、そのパラメータに等しいことに注意。
図４のフォームを閉じると、端末に次のメッセージが表示されて、スクリプトの実行終了となる。
Results.txt
was saved.
出力ファイルを開くと、以下のような内容である。
Summary...
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu
-0.315 1.020 -2.266 1.435 0.037 0.026 829.0 911.0 1.0
a2
1.240 1.185 -0.879 3.556 0.040 0.030 880.0 1057.0 1.0
a3
2.003 1.106 0.040 4.173 0.038 0.028 909.0 957.0 1.0
b2
-0.299 1.640 -3.369 2.731 0.055 0.040 916.0 1132.0 1.0
b3
2.565 1.057 0.636 4.526 0.037 0.027 855.0 950.0 1.0
g22
-0.207 1.898 -3.677 3.402 0.061 0.043 963.0 1280.0 1.0
g23
-1.243 1.260 -3.570 1.177 0.042 0.031 906.0 1196.0 1.0
g32
-0.952 1.849 -4.587 2.415 0.058 0.041 1027.0 1261.0 1.0
g33 -3.000 1.245 -5.321 -0.730 0.040 0.029 985.0 1411.0 1.0
Lmbd[0] 1.119 1.010 0.020 3.061 0.031 0.022 829.0 911.0 1.0
Lmbd[1] 0.980 0.992 0.000 2.786 0.020 0.014 1551.0 1135.0 1.0
Lmbd[2] 9.957 3.061 4.530 15.718 0.044 0.031 4728.0 3115.0 1.0
Lmbd[3] 3.011 1.733 0.307 6.144 0.030 0.021 3008.0 2244.0 1.0
Lmbd[4] 1.989 1.403 0.066 4.566 0.021 0.015 4208.0 3113.0 1.0
Lmbd[5] 9.957 3.160 4.397 15.698 0.047 0.033 4502.0 3360.0 1.0
Lmbd[6] 5.907 2.431 1.759 10.315 0.035 0.025 4637.0 3285.0 1.0
Lmbd[7] 2.024 1.414 0.091 4.608 0.022 0.016 4130.0 2702.0 1.0
Lmbd[8] 3.997 2.008 0.721 7.678 0.030 0.022 4074.0 3111.0 1.0
LogLmbd[0] -0.315 1.020 -2.266 1.435 0.037 0.026 829.0 911.0 1.0
LogLmbd[1] -0.614 1.308 -3.075 1.536 0.038 0.033 1551.0 1135.0 1.0
LogLmbd[2] 2.250 0.318 1.662 2.834 0.005 0.003 4728.0 3115.0 1.0
LogLmbd[3] 0.926 0.630 -0.251 2.065 0.012 0.008 3008.0 2244.0 1.0
LogLmbd[4] 0.419 0.794 -1.097 1.775 0.013 0.010 4208.0 3113.0 1.0
LogLmbd[5] 2.247 0.327 1.605 2.806 0.005 0.004 4502.0 3360.0 1.0
LogLmbd[6] 1.688 0.433 0.889 2.466 0.006 0.005 4637.0 3285.0 1.0
LogLmbd[7] 0.437 0.796 -1.076 1.799 0.013 0.010 4130.0 2702.0 1.0
LogLmbd[8] 1.252 0.541 0.292 2.249 0.009 0.006 4074.0 3111.0 1.0
P(a3 > 0.0) = 0.979
P(b3 > 0.0) = 1.000
P(g33 < 0.0) = 0.998
リスト１．ポアッソン回帰分析Stanスクリプト(Poisson.stan)
data{
int N;
array[N] int F;
array[N] int Xr1;
array[N] int Xr2;
array[N] int Xr3;
array[N] int Xc1;
array[N] int Xc2;
array[N] int Xc3;
}
parameters{
real mu;
real a2;
real a3;
real b2;
real b3;
real g22;
real g23;
real g32;
real g33;
}
transformed parameters{
array[N] real Lmbd;
array[N] real LogLmbd;
for (i in 1:N)
LogLmbd[i] = mu + a2 * Xr2[i] + a3 * Xr3[i]
+ b2 * Xc2[i] + b3 * Xc3[i]
+
g22 * Xr2[i] * Xc2[i] + g23 * Xr2[i] * Xc3[i]
+ g32 * Xr3[i] * Xc2[i] + g33 * Xr3[i] * Xc3[i];
for (i in 1:N)
Lmbd[i] = exp(LogLmbd[i]);
}
model{
mu ~ normal(0.0, 10.0);
a2 ~ normal(0.0, 10.0);
a3 ~ normal(0.0, 10.0);
b2 ~ normal(0.0, 10.0);
b3 ~ normal(0.0, 10.0);
g22 ~ normal(0.0, 10.0);
g23 ~ normal(0.0, 10.0);
g32 ~ normal(0.0, 10.0);
g33 ~ normal(0.0, 10.0);
for (i in 1:N)
F[i]
~ poisson(Lmbd[i]);
}
リスト２．ポアッソン回帰分析のPythonスクリプト（PoissonRefSta.py）
import stan
import matplotlib.pyplot as plt
import numpy as np
import arviz as az
import seaborn as sb
RawData = [[1, 1, 0, 0, 1, 0, 0],
[1, 1, 0, 0, 0, 1, 0],
[10, 1, 0, 0, 0, 0, 1],
[3, 0, 1, 0, 1, 0, 0],
[2, 0, 1, 0, 0, 1, 0],
[10, 0, 1, 0, 0, 0, 1],
[6, 0, 0, 1, 1, 0, 0],
[2, 0, 0, 1, 0, 1, 0],
[4, 0, 0, 1, 0, 0, 1]]
F = []
Xr1 = []
Xr2 = []
Xr3 = []
Xc1 = []
Xc2 = []
Xc3 = []
for v in RawData:
F.append(v[0])
Xr1.append(v[1])
Xr2.append(v[2])
Xr3.append(v[3])
Xc1.append(v[4])
Xc2.append(v[5])
Xc3.append(v[6])
for i in range(9):
print(' ', F[i], ' ', Xr1[i], ' ', Xr2[i], ' ', Xr3[i],
' ', Xc1[i], ' ', Xc2[i], ' ', Xc3[i])
N = 9
Data = {'N': N, 'F': F, 'Xr1': Xr1,
'Xr2': Xr2, 'Xr3': Xr3,
'Xc1': Xc1, 'Xc2': Xc2, 'Xc3': Xc3}
with open('Poisson.stan', 'r') as
f_stan:
sm =
stan.build(f_stan.read(), data = Data)
fit = sm.sample()
f = open('Results.txt', 'w')
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data)
print('Summary =\n',smry)
f.write('\nSummary...\n')
f.write(smry.__str__())
f.write('\n')
d_fm = fit.to_frame()
az.plot_trace(i_data, var_names =
('a2', 'a3', 'b2', 'b3'))
plt.show()
a3 = d_fm['a3']
cnt_a3_pos = (a3 > 0.0).sum()
print('\nP(a3 > 0.0) =
{0:.3f}'.format(cnt_a3_pos / len(a3)))
f.write('\nP(a3 > 0.0) =
{0:.3f}\n'.format(cnt_a3_pos / len(a3)))
b3 = d_fm['b3']
cnt_b3_pos = (b3 > 0.0).sum()
print('P(b3 > 0.0) =
{0:.3f}'.format(cnt_b3_pos / len(b3)))
f.write('\nP(b3 > 0.0) =
{0:.3f}\n'.format(cnt_b3_pos / len(b3)))
sb.kdeplot(d_fm['a2'], label = 'a2')
sb.kdeplot(d_fm['a3'], label = 'a3')
sb.kdeplot(d_fm['b2'], label = 'b2')
sb.kdeplot(d_fm['b3'], label = 'b3')
plt.legend()
plt.show()
g33 = d_fm['g33']
print(len(g33))
print(np.shape(g33))
cnt_g33_neg = (g33 < 0.0).sum()
print('P(g33 < 0.0) =
{0:.3f}'.format(cnt_g33_neg / len(g33)))
f.write('\nP(g33 < 0.0) =
{0:.3f}\n'.format(cnt_g33_neg / len(g33)))
sb.kdeplot(d_fm['g22'], label = 'g22')
sb.kdeplot(d_fm['g23'], label = 'g23')
sb.kdeplot(d_fm['g32'], label = 'g32')
sb.kdeplot(d_fm['g33'], label = 'g33')
plt.legend()
plt.show()
lmbd_hat = np.empty(N)
for i in range(N):
lmbd_hat[i] =
np.median(d_fm[f'Lmbd.{i+1}'])
plt.plot(lmbd_hat[i], F[i],
'o', c='b')
min_lmbd_hat = np.min(lmbd_hat)
max_lmbd_hat = np.max(lmbd_hat)
plt.plot([min_lmbd_hat, max_lmbd_hat],
[min_lmbd_hat, max_lmbd_hat],
label = '$F = \hat{\lambda}$')
plt.xlabel('$\hat{\lambda}$', fontsize
= 16)
plt.ylabel('F', fontsize = 14)
plt.legend(fontsize = 16)
plt.tight_layout()
plt.show()
f.close()
print('Results.txt was saved.')