拙著「Pythonでデータ分析」の第15章「PyMCによるポアッソン回帰モデル分析」更新
スクリプトをリスト１に示すように更新した。スクリプトファイルは、圧縮ファイルc15file.zipとして用意した。ダウンロード展開すれば自由に利用できる。
スクリプトの実行は、Python
3.10/Anaconda/Ubuntu/WSL/Windowsで行った。すなわち、WindowsにWSL (Windows Subsystem for Linux)をインストールするとUbuntuもインストールされるので、そのUbuntuにAnacondaをインストールしてPython 3.10仮想環境を用意した。Python 3.10仮想環境においてコマンド
conda install pymc3
によりPyMC3をインストールした。詳しくは、このウェブサイトを参照されたい。
リスト１のスクリプトをPython 3.10仮想環境において実行した。
(py310) ・・・/c15file$ python
PssnRegPyMC3.py
スクリプトを実行すると、MCMCサンプリングが始まり、終了後、トレースグラフが表示される（図１）。
図１
図１のフォームを閉じると、図２のグラフが表示される。
図２
主効果、a2、a3、b2、b3の事後分布のグラフである。拙著の表8.2.2において、セル（小学生の時に知った、心理学≒臨床心理学）を基準に主効果などが表されているが（拙著の表12.2）、a3とb3の事後分布が正の領域に位置していることは、「高校生のとき」および「心理学＞臨床心理学」の条件での度数が多いことに対応している。
図２のフォームを閉じると、図３のグラフが表示される。
図３
交互作用g32の事後分布は負の領域に寄っており、g33の分布はほぼ負の領域にある。主効果のみで表される度数に比べて、セル（高校生、その他）および（高校生、心理学＞臨床心理学）の度数が少ないことが表されている。
図３のフォームを閉じると、図４のグラフが表示される。
図４
ポアッソン分布（拙著、p.203）
におけるパラメータ
の推定値（事後分布の中央値）とデータ値
の散布図である。ポアッソン分布の平均値は、そのパラメータに等しいことに注意。
図４のフォームを閉じると、端末に次のメッセージが表示されて、スクリプトの実行終了となる。
summary.txt
was saved.
出力ファイルの内容は、以下のようである。
Summary...
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu
-0.383 1.105 -2.464 1.527 0.027 0.021 1927.0 1750.0 1.0
a2
1.294 1.250 -0.993 3.744 0.029 0.023 2009.0 1938.0 1.0
a3
2.083 1.181 0.096 4.466 0.028 0.022 1993.0 2013.0 1.0
b2
-0.207 1.667 -3.399 2.972 0.034 0.024 2454.0 2771.0 1.0
b3
2.626 1.151 0.716 4.951 0.028 0.021 1911.0 1832.0 1.0
g22 -0.277 1.931 -3.823 3.495 0.037 0.026 2731.0 3772.0 1.0
g23 -1.281 1.323 -3.874 1.121 0.030 0.024 2064.0 2353.0 1.0
g32 -1.059 1.900 -4.775 2.478 0.036 0.026 2801.0 3528.0 1.0
g33 -3.062 1.335 -5.698 -0.746 0.030 0.023 2142.0 2124.0 1.0
Lmbd[0] 1.094 1.015 0.004 2.960 0.019 0.013 1927.0 1750.0 1.0
Lmbd[1] 0.995 0.991 0.000 2.802 0.012 0.008 4533.0 3470.0 1.0
Lmbd[2] 9.917 3.135 4.416 15.648 0.027 0.020 12453.0 9673.0 1.0
Lmbd[3] 2.964 1.716 0.441 6.127 0.017 0.012 8779.0 6773.0 1.0
Lmbd[4] 1.997 1.400 0.074 4.516 0.012 0.009 12448.0 9437.0 1.0
Lmbd[5] 10.042 3.146 4.580 16.057 0.028 0.020 12732.0 9698.0 1.0
Lmbd[6] 5.969 2.454 1.671 10.363 0.022 0.016 12411.0 9630.0 1.0
Lmbd[7] 2.011 1.416 0.030 4.496 0.012 0.009 12534.0 8864.0 1.0
Lmbd[8] 4.032 2.032 0.790 7.700 0.018 0.013 12339.0 9818.0 1.0
95%CI for a3 = [0.13184, 4.73770]
P(a3 > 0.0) = 0.983
P(b3 > 0.0) = 0.999
P(g33 < 0.0) = 0.997
リスト１．ポアッソン回帰モデル（PssnRegPyMC3.py）
import matplotlib.pyplot
as plt
import seaborn
as sb
import numpy
as np
import pymc3 as pm
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])
with pm.Model()
as PoiModel:
mu = pm.Normal('mu',
mu = 0.0, sd = 10.0)
a2 = pm.Normal('a2',
mu = 0.0, sd = 10.0)
a3 = pm.Normal('a3',
mu = 0.0, sd = 10.0)
b2 = pm.Normal('b2',
mu = 0.0, sd = 10.0)
b3 = pm.Normal('b3',
mu = 0.0, sd = 10.0)
g22 = pm.Normal('g22',
mu = 0.0, sd = 10.0)
g23 = pm.Normal('g23',
mu = 0.0, sd = 10.0)
g32 = pm.Normal('g32',
mu = 0.0, sd = 10.0)
g33 = pm.Normal('g33',
mu = 0.0, sd = 10.0)
LogLmbd
= mu + a2 * Xr2 + a3 * Xr3 + b2 * Xc2 + b3 * Xc3 \
+ g22 * Xr2 * Xc2 + g23 * Xr2 * Xc3 \
+ g32 * Xr3 * Xc2 + g33 * Xr3 * Xc3
#Lmbd
= np.exp(LogLmbd)
Lmbd
= pm.Deterministic('Lmbd', np.exp(LogLmbd))
F_obs
= pm.Poisson('F_obs', mu = Lmbd, observed = F)
trace = pm.sample(draws
= 3000, tune = 3000, target_accept = 0.95)
pm.plot_trace(trace,
var_names = ['a2', 'a3', 'b2', 'b3'])
plt.show()
summary = pm.summary(trace)
print(summary)
f = open('summary.txt', 'w')
f.write('Summary...\n{}'.format(summary))
smpl_a3 = trace['a3']
a3_p2p5 = np.percentile(smpl_a3, 2.5)
a3_p97p5 = np.percentile(smpl_a3, 97.5)
print('\n95%CI for a3 =
[{0:.5f}, {1:.5f}]'.format(a3_p2p5, a3_p97p5))
f.write('\n\n95%CI
for a3 = [{0:.5f}, {1:.5f}]\n'.format(a3_p2p5,
a3_p97p5))
cnt_a3_pos = (smpl_a3 >
0.0).sum()
print('\nP(a3
> 0.0) = {0:.3f}'.format(cnt_a3_pos / len(smpl_a3)))
f.write('\nP(a3 > 0.0) = {0:.3f}\n'.format(cnt_a3_pos
/ len(smpl_a3)))
sb.kdeplot(trace['a2'],
label = 'a2')
sb.kdeplot(trace['a3'],
label = 'a3')
sb.kdeplot(trace['b2'],
label = 'b2')
sb.kdeplot(trace['b3'],
label = 'b3')
plt.legend()
plt.show()
smpl_b3 = trace['b3']
cnt_b3_pos = (smpl_b3 >
0.0).sum()
print('P(b3 > 0.0) =
{0:.3f}'.format(cnt_b3_pos / len(smpl_b3)))
f.write('P(b3
> 0.0) = {0:.3f}\n'.format(cnt_b3_pos / len(smpl_b3)))
smpl_g33 = trace['g33']
cnt_g33_neg = (smpl_g33 <
0.0).sum()
print('P(g33 < 0.0) =
{0:.3f}'.format(cnt_g33_neg / len(smpl_g33)))
f.write('P(g33
< 0.0) = {0:.3f}\n'.format(cnt_g33_neg / len(smpl_g33)))
sb.kdeplot(trace['g22'],
label = 'g22')
sb.kdeplot(trace['g23'],
label = 'g23')
sb.kdeplot(trace['g32'],
label = 'g32')
sb.kdeplot(trace['g33'],
label = 'g33')
plt.legend()
plt.show()
N= 9
lmbd_hat
= np.empty(N)
for i
in range(N):
lmbd_hat[i] = np.median(trace['Lmbd'][:,i])
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()