ポアッソン分布による２つの独立なデータの平均値のベイズ分析
度数データの分析にポアッソン分布が用いられることがある。Spiegelhalter (2019)は、EnglandとWalesにおける1998年〜2016年の各年ごとの殺人事件の件数を取り上げているが、このデータをポアッソン分布を用いた回帰モデルでベイズ分析すると、2005年が特別な年であったことが示された。
警視庁のホームページ「各種交通事故発生状況（令和6年上半期）」からダウンロードしたPDFファイル「小学生の交通人身事故発生状況（令和6年上半期）更新日：2024年7月18日」に掲載されている月単位のデータを読み取りポアッソン分布による多項式回帰モデルによるベイズ分析を行うと、３次の多項式が得られ、2月が事故が最低で5月６月が多くなっていることが示された。
ポアッソン分布は次式で与えられる（Gelman et al., 2014）。
ポアッソン分布による２つの独立なデータのベイズ分析を行うスクリプトをリスト１とリスト２のように用意した。スクリプトは、リスト１とリスト２に分割した。これは、PyStan3はWindows上で直接実行できないので、WindowsにWSL(Windows Subsystem for Linux)をインストールしてWSL上で実行する必要があり、グラフィック（GUI）がWSLでは利用できなかった（2025.02）のでグラフ描画はWindows上で行ったからである。ファイル類は、まとめて圧縮ファイルbapoissonfiles.zipとした。
リスト１のスクリプトを、miniconda/Ubuntu24/WSL/Windows10のPython3.12仮想環境で実行した。Minicondaの使用は、Anacondaと同様である。このウェブサイトを参照されたい。
リスト１のスクリプトを実行すると、入力データファイル名の設定が求められる。
(py312) *****/bapoissonfiles$
python ba_poisson.py
Input data file(*.csv) =
入力データファイルは、CSV形式で用意する（図１）
・
・
・
図１ CSV形式のファイル（dataN30.csv）
1行目に変数名を設定する。１列目は空欄でもよい。データは２行目から設定する。１列目はデータの識別用であり、分析では用いられない。２列目と３列目にデータ値を設定する。データ数が２つのグループで異なるときは、少ない方のデータは空欄にしておけばよい。空欄はスクリプトの実行時に無視される。
作成したデータファイルは、CSVファイルとして保存する。ファイル名の拡張子として.csvを選べばよい。データファイルとリスト１およびリスト２のファイルは、同じフォルダ内に置いておく。
次の例では図１のファイルを入力データとして設定している。
(py312) *****/bapoissonfiles$
python ba_poisson.py
Input data file(*.csv) = dataN30.csv
入力ファイル名を設定すると、ファイルが読み込まれ、計算が始まる。
MCMCサンプリングが終了すると、MCMCサンプルは保存されてスクリプトの実行終了となる。
fit.pkl was saved.
(py312) *****/bapoissonfiles$
リスト１の実行によるMCMCサンプルのファイルを読み込んで分析を行うために、リスト２のスクリプトを実行する。WSLではグラフ描画ができなかったので、リスト２のスクリプトは直接Windows上でPython3.13.2により実行した。
リスト１の実行後、リスト２のスクリプトを実行すると、リスト１の実行により生成されたMCMCサンプルのファイルが読み込まれて、パラメータlmbdの事後分布のグラフが表示される（図２）。
図２
タイトル部には、事後分布の中央値が表示されている（Gelman et al., 2021）。
図２のWindowを閉じると、図３のグラフが表示される。
図３
パラメータの差lmbd1-lmbd2の事後分布である。タイトル部にlmbd1>lmbd2の事後確率が0.00000であること、lmbsd1<lmbd2の事後確率が1.00000ことが表示されている。ほぼ確実に、パラメータlmbd1の方が小さいと言えそうである。パラメータlmbdは、ポアッソン分布の平均値であるので、data1の平均値の方が小さいといえそうである
図３のWindowを閉じると、図４のグラフが表示される。
図４
データのヒストグラムと、事後分布の中央値をパラメータ値とするポアッソン分布のprobability mass functionである。
リスト１ PyStan3によるMCMCサンプリング（ba_poisson.py）
import numpy
as np
import pandas as pd
import stan
import pickle
inf_nm = input('Input data file(*.csv) = ')
dframe = pd.read_csv(inf_nm)
# データの読み込み
with open('in_dframe.pkl',
'wb') as f:
pickle.dump(dframe, f)
# 読み込みデータの保存
df_keys = dframe.keys()
x_name = df_keys[1]
y_name = df_keys[2]
print(x_name,
y_name)
x = dframe[x_name].values
x = x[np.isnan(x)==False] # 空セルデータの削除
x = np.array(x,
dtype=int)
print('X =\n', x)
y = dframe[y_name].values
y = y[np.isnan(y)==False] # 空セルデータの削除
y = np.array(y,
dtype=int)
print('Y =\n', y)
mean1 = np.mean(x)
sd1 = np.std(x)
print('mean1 =', mean1, ' sd1 =', sd1)
mean2 = np.mean(y)
sd2 = np.std(y)
print('mean2 =', mean2, ' sd2 =', sd2)
max1 = np.max(x)
max2 = np.max(y)
print('max1 =',max1)
print('max2 =',max2)
stan_code = """
data {
int n1;
array[n1] int x;
int n2;
array[n2] int y;
real
max1;
real
max2;
}
parameters {
real<lower=0, upper=1> vlmbd1;
real<lower=0, upper=1> vlmbd2;
}
transformed parameters {
/*
Parameters rescaled
*/
real
lmbd1;
real
lmbd2;
lmbd1 = 0.001 + vlmbd1*max1;
lmbd2 = 0.001 + vlmbd2*max2;
}
model {
vlmbd1 ~ beta(1,1);
vlmbd2 ~ beta(1,2);
x ~ poisson(lmbd1);
y ~ poisson(lmbd2);
}
generated quantities {
real
diff_lmbd;
diff_lmbd = lmbd1 - lmbd2;
}
"""
Data = {'n1':len(x), 'x':x, 'max1':max1,
'n2':len(y), 'y':y, 'max2':max2}
sm = stan.build(stan_code, data=Data)
fit = sm.sample()
pd_fit_frame = fit.to_frame()
with open('fit.pkl',
'wb') as f:
pickle.dump(pd_fit_frame,f) # MCMC samples saved
print('fit.pkl
was saved.')
リスト２ MCMCサンプルの分析（draw_poisson.py）
import numpy
as np
import pandas as pd
import pickle
import matplotlib.pyplot
as plt
import seaborn
as sb
import scipy.stats
as ss
with open('in_dframe.pkl',
'rb') as f:
dframe
= pickle.load(f)
df_keys = dframe.keys()
x_name = df_keys[1]
y_name = df_keys[2]
print(x_name,
y_name)
x = dframe[x_name].values
x = x[np.isnan(x)==False] # 空セルデータの削除
x = np.array(x,
dtype=int)
print('X =\n', x)
y = dframe[y_name].values
y = y[np.isnan(y)==False] # 空セルデータの削除
y = np.array(y,
dtype=int)
print('Y =\n', y)
mean1 = np.mean(x)
sd1 = np.std(x)
print('mean1 =', mean1, ' sd1 =', sd1)
mean2 = np.mean(y)
sd2 = np.std(y)
print('mean2 =', mean2, ' sd2 =', sd2)
max1 = np.max(x)
max2 = np.max(y)
print('max1 =',max1)
print('max2 =',max2)
with open('fit.pkl',
'rb') as f:
pd_fit_frame
= pickle.load(f)
lmbd1 = np.median(pd_fit_frame['lmbd1'])
lmbd2 = np.median(pd_fit_frame['lmbd2'])
"""
パラメータlmbd1とlmbd2の事後分布
"""
sb.kdeplot(pd_fit_frame['lmbd1'],
label='lmbd1')
sb.kdeplot(pd_fit_frame['lmbd2'],
label='lmbd2')
plt.title('Posterior distributions' +
f'\nlmbd1={lmbd1:.1f},
lmbd2={lmbd2:.1f}')
plt.xlabel('lmbd', fontsize=14)
plt.legend()
plt.show()
"""
パラメータの差lmbd1-lmbd2の事後分布
"""
diff_lmbd = pd_fit_frame['diff_lmbd']
pos_diff = diff_lmbd
> 0
neg_diff = diff_lmbd
< 0
prob_pos = np.mean(pos_diff)
prob_neg = np.mean(neg_diff)
plt.title('Posterior distribution' +
f'\nP(lmbd1>lmbd2)={prob_pos:.5f}, P(lmbd1<lmbd2)={prob_neg:.5f}')
sb.kdeplot(diff_lmbd)
plt.xlabel('lmbd1-lmbd2', fontsize=14)
plt.show()
"""
データのヒストグラムと
事後分布の中央値をパラメータ値とする
Probabilituy Mass Functions
"""
maxxy = max1 if max1 > max2 else max2
minxy = np.min(x)
if np.min(x) < np.min(y)
else np.min(y)
xy_coord = np.arange(minxy, maxxy+1) # np.linspace(minxy, maxxy, 1000)
y1_coord = ss.poisson.pmf(xy_coord, mu=lmbd1)
y2_coord = ss.poisson.pmf(xy_coord, mu=lmbd2)
plt.hist(x, density=True, color='b',
alpha=0.3, label=x_name)
plt.hist(y, density=True, color='orange',
alpha=0.3, label=y_name)
plt.title(f'lmbd1={lmbd1:.1f}, lmbd2={lmbd2:.1f}')
plt.plot(xy_coord,
y1_coord, 'o', color='b', label=x_name)
plt.plot(xy_coord,
y2_coord, 'o', color='orange', label=y_name)
plt.ylabel('Prob.Mass/Density')
plt.legend()
plt.tight_layout()
plt.show()
ベイズ分析の基本的考え方の説明を岡本（2019）、pp.169-177、岡本（2014）、pp.227-232で行っている。
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data
analysis, 3rd. Ed. CRC Press.
Gelman, A., Hill, J., & Vehtari,
A. (2021). Regression and other stories. Cambridge University
岡本安晴（2019）「いまさら聞けないPythonでデータ分析」丸善出版
Spiegelhalter, D. (2019).The art of statistics: Learning from data. Pelican