Up

ポアッソン分布による2つの独立なデータの平均値のベイズ分析

 

度数データの分析にポアッソン分布が用いられることがある。Spiegelhalter (2019)は、EnglandWalesにおける1998年〜2016年の各年ごとの殺人事件の件数を取り上げているが、このデータをポアッソン分布を用いた回帰モデルでベイズ分析すると、2005年が特別な年であったことが示された。

警視庁のホームページ「各種交通事故発生状況(令和6年上半期)」からダウンロードしたPDFファイル「小学生の交通人身事故発生状況(令和6年上半期)更新日:2024718日」に掲載されている月単位のデータを読み取りポアッソン分布による項式回帰モデルによるベイズ分析を行うと、3次の多項式が得られ、2月が事故が最低で5月6月が多くなっていることが示された。

 

ポアッソン分布は次式で与えられる(Gelman et al., 2014)。

 

 

ポアッソン分布による2つの独立なデータのベイズ分析を行うスクリプトをリスト1リスト2のように用意した。スクリプトは、リスト1とリスト2に分割した。これは、PyStan3Windows上で直接実行できないので、WindowsWSL(Windows Subsystem for Linux)をインストールしてWSL上で実行する必要があり、グラフィック(GUI)がWSLでは利用できなかった(2025.02)のでグラフ描画はWindows上で行ったからである。ファイル類は、まとめて圧縮ファイルbapoissonfiles.zipとした。

リスト1のスクリプトを、miniconda/Ubuntu24/WSL/Windows10Python3.12仮想環境で実行した。Minicondaの使用は、Anacondaと同様である。このウェブサイトを参照されたい。

リスト1のスクリプトを実行すると、入力データファイル名の設定が求められる。

 

(py312) *****/bapoissonfiles$ python ba_poisson.py

Input data file(*.csv) =

 

入力データファイルは、CSV形式で用意する(図1)

図1 CSV形式のファイル(dataN30.csv

 

1行目に変数名を設定する。1列目は空欄でもよい。データは2行目から設定する。1列目はデータの識別用であり、分析では用いられない。2列目と3列目にデータ値を設定する。データ数が2つのグループで異なるときは、少ない方のデータは空欄にしておけばよい。空欄はスクリプトの実行時に無視される。

作成したデータファイルは、CSVファイルとして保存する。ファイル名の拡張子として.csvを選べばよい。データファイルとリスト1およびリスト2のファイルは、同じフォルダ内に置いておく。

 

次の例では図1のファイルを入力データとして設定している。

 

(py312) *****/bapoissonfiles$ python ba_poisson.py

Input data file(*.csv) = dataN30.csv

 

入力ファイル名を設定すると、ファイルが読み込まれ、計算が始まる。

MCMCサンプリングが終了すると、MCMCサンプルは保存されてスクリプトの実行終了となる。

 

fit.pkl was saved.

(py312) *****/bapoissonfiles$

 

 

リスト1の実行によるMCMCサンプルのファイルを読み込んで分析を行うために、リスト2のスクリプトを実行する。WSLではグラフ描画ができなかったので、リスト2のスクリプトは直接Windows上でPython3.13.2により実行した。

リスト1の実行後、リスト2のスクリプトを実行すると、リスト1の実行により生成されたMCMCサンプルのファイルが読み込まれて、パラメータlmbdの事後分布のグラフが表示される(図2)。

図2

 

タイトル部には、事後分布の中央値が表示されている(Gelman et al., 2021)。

図2のWindowを閉じると、図3のグラフが表示される。

図3

 

パラメータの差lmbd1-lmbd2の事後分布である。タイトル部にlmbd1>lmbd2の事後確率が0.00000であること、lmbsd1<lmbd2の事後確率が1.00000ことが表示されている。ほぼ確実に、パラメータlmbd1の方が小さいと言えそうである。パラメータlmbdは、ポアッソン分布の平均値であるので、data1の平均値の方が小さいといえそうである

図3のWindowを閉じると、図4のグラフが表示される。

図4

 

データのヒストグラムと、事後分布の中央値をパラメータ値とするポアッソン分布のprobability mass functionである。

 

 

 

リスト1 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.')

 

 

 

リスト2 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'])

"""

       パラメータlmbd1lmbd2の事後分布

"""

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()

 

 

 

 

ベイズ分析の基本的考え方の説明を岡本(2019pp.169-177岡本(2014pp.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

岡本安晴(2014)「心理学データ分析と測定」勁草書房

岡本安晴(2019)「いまさら聞けないPythonでデータ分析」丸善出版

Spiegelhalter, D. (2019).The art of statistics: Learning from data. Pelican

 

 

 

Home