ポアッソン分布によるクロス表のベイズ分析
本ウェブサイトでは、クロス表を分散分析の形式でポアッソン分布を用いて回帰モデルを設定してベイズ分析を行う。ポアッソン分布を用いて多項式による回帰分析を行う例は、小学生の交通人身事故発生状況の月ごとの変化をポアッソン分布を用いた回帰分析を行い、２月の事故が少なく、５月の事故が多いことを最尤法を用いた分析で示したもの、ロンドンの殺人事件の2005年の件数が異常に少ないことをポアッソン分布を用いたベイズ分析による回帰分析で示したものも用意している。
本ウェブサイトにおいては、表１の形式のデータを分散分析のデザインでポアッソン分布を用いたベイズ分析を行う。
|
表１ |
|
|
|
|
|
|
心理学≒臨床心理学 |
その他 |
心理学＞臨床心理学 |
計 |
|
小学生 |
1 |
1 |
10 |
12 |
|
中学生 |
3 |
2 |
10 |
15 |
|
高校生 |
6 |
2 |
4 |
12 |
|
計 |
10 |
5 |
24 |
39 |
表１のデータは、岡本安晴（2018）最終講義「『「計量心理学」落穂拾い』（http://y-okamoto-psy1949.la.coocan.jp/booksetc/Lec201803/Lecture201803.pdfの表を岡本（2019）表8.2.2に合わせて編集したものである。
岡本（2019）では、期待度数
を
で、
を
で、
を
で、
を
で表したとき、
と置き、観測度数
が、パラメータのポワッソン分布から得られたとしている。すなわち、
である。
上のモデルに基づいてベイズ分析を行うスクリプトを、リスト１とリスト２に分けて示す。スクリプトファイルは、まとめて圧縮ファイルpoissonfiles.zipとした。
StanのPython用モジュールPyStan 3は、Windows上では直接利用できない。WindowsにWSL（Windows Subsystem for Linux）をインストールすれば用いることができる。WindowsでのPyStan 3の準備については、このウェブサイトなどで説明している。WindowsでPyStan 3を利用するときは、PyStan 3によるMCMCサンプリングをUbuntu/WSL/Windowsで行うが、WSLでは描画ができなかったので、描画はWindows上で行った。リスト１のスクリプトはMiniconda/Ubuntu24/WSL/Windows10で実行して計算結果をファイルに保存し、リスト１の出力ファイルをリスト２で読み込み分析した。リスト２のスクリプトは、Python3.13.2/Windows10で実行した。
表１のデータは、リスト１では次のコードで与えられている。
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]]
リストRawDataの各要素の第１要素に表１の数値が置かれ、続く６つの要素は行の指標値（ダミー変数）と列の指標値である（岡本、2019、表12.3）。
MCMCサンプリングは、次のコードで行っている。
sm = stan.build(StanCode, data=Data)
fit = sm.sample()
MCMCサンプリングの結果は、次のコードでファイルtemp.pklに保存される。
inf_data =
az.from_pystan(posterior=fit, posterior_model=sm)
print(az.summary(inf_data))
df_data = fit.to_frame()
with open('temp.pkl', 'wb') as f:
pickle.dump((inf_data,
df_data), f)
保存されたファイルは、リスト２のスクリプトにより読み込まれて分析される。リスト２は、Python 3.13.2/Windows10で実行した。
リスト２を実行すると、次のコードによりリスト１の実行によって保存されたファイルが読み込まれる。
with open('temp.pkl', 'rb') as f:
(inf_data, df_data) =
pickle.load(f)
ファイルが読み込まれると、トレース図が表示される（図１）。
図１
図１のWindowを閉じると、主効果パラメータの事後分布のグラフが表示される（図２）。
図２
図２のWindowを閉じると、交互作用パラメータの事後分布が表示される（図３）。
図３
図３のWindowsを閉じると、以下の出力が端末に表示されて、スクリプトの実行終了となる。
mean sd hdi_3% hdi_97%
mu -0.304 0.992 -2.167 1.449
a2 1.223 1.153 -1.036 3.294
a3 2.003 1.061 0.132 4.052
b2 -0.204 1.477 -3.109 2.456
b3 2.557 1.033 0.715 4.542
g22 -0.308 1.770 -3.517 3.226
g23 -1.224 1.232 -3.656 1.003
g32 -1.069 1.701 -4.188 2.219
g33 -3.012 1.214 -5.284 -0.790
P(a3 > 0.0) = 0.987
P(b3 > 0.0) = 0.999
P(g33 < 0.0) = 0.998
94％HDIが０を含まないパラメータの箇所を赤字で示した。パラメータa3が正であることは、「高校生」の行の数値が他の行より大きいことを、パラメータb3が正であることは、「心理学が臨床心理学を含むものであること」の列の数値が他の列より大きいことを表し、パラメータg33の数値が負であることは、セル（高校生、心理学が臨床心理学を含む）の数値が、主効果で説明される値より小さいことを表している。
import stan
import arviz as az
import pickle
StanCode = """
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];
Lmbd = exp(LogLmbd);
}
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);
F ~
poisson(Lmbd);
}
"""
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}
sm = stan.build(StanCode, data=Data)
fit = sm.sample()
print(fit)
inf_data =
az.from_pystan(posterior=fit, posterior_model=sm)
print(az.summary(inf_data))
df_data = fit.to_frame()
with open('temp.pkl', 'wb') as f:
pickle.dump((inf_data,
df_data), f)
import matplotlib.pyplot as plt
import seaborn as sb
import numpy as np
import arviz as az
import pickle
"""
リスト１のスクリプトで保存されたファイルの読み込み
"""
with open('temp.pkl', 'rb') as f:
(inf_data, df_data) =
pickle.load(f)
"""
トレース図の描画
"""
az.plot_trace(inf_data,
var_names=['mu', 'a2', 'a3', 'b2', 'b3',
'g22', 'g23', 'g32', 'g33'],
figsize=(10,7))
plt.tight_layout()
plt.show()
"""
主効果パラメータの事後分布
"""
sb.kdeplot(df_data['a2'], label='a2')
sb.kdeplot(df_data['a3'], label='a3')
sb.kdeplot(df_data['b2'], label='b2')
sb.kdeplot(df_data['b3'], label='b3')
plt.xlabel('Main Effects',
fontsize=14)
plt.legend()
plt.tight_layout()
plt.show()
"""
交互作用パラメータの事後分布
"""
sb.kdeplot(df_data['g22'], label='g22')
sb.kdeplot(df_data['g23'],
label='g23')
sb.kdeplot(df_data['g32'],
label='g32')
sb.kdeplot(df_data['g33'],
label='g33')
plt.xlabel('Interaction', fontsize=14)
plt.legend()
plt.tight_layout()
plt.show()
"""
パラメータの統計量
"""
print(az.summary(inf_data,
kind='stats',
var_names=['mu', 'a2', 'a3', 'b2', 'b3',
'g22', 'g23', 'g32', 'g33']))
a3 = df_data['a3']
cnt_a3_pos = (a3 > 0.0).sum()
print('\nP(a3 > 0.0) =
{0:.3f}'.format(cnt_a3_pos / len(a3)))
b3 = df_data['b3']
cnt_b3_pos = (b3 > 0.0).sum()
print('P(b3 > 0.0) =
{0:.3f}'.format(cnt_b3_pos / len(b3)))
g33 = df_data['g33']
cnt_g33_neg = (g33 < 0.0).sum()
print('P(g33 < 0.0) =
{0:.3f}'.format(cnt_g33_neg / len(g33)))