Up

ポアッソン分布によるクロス表のCmdStanPyを用いたベイズ分析

拙著「いまさら聞けないPythonでデータ分析」では、PyStanを用いてStanスクリプトの説明をした。ここでは、CmdStanPyを用いてStanによるベイズ分析を試みる。CmdStanPyの簡単なインストール法の説明のウェブサイトも用意した。

表１の形式のデータを分散分析のデザインでポアッソン分布を用いたベイズ分析を行う。

表１
	心理学≒臨床心理学	その他	心理学＞臨床心理学	計
小学生	1	1	10	12
中学生	3	2	10	15
高校生	6	2	4	12
計	10	5	24	39

岡本（2019）では、期待度数をで、をで、をで、をで表したとき、

と置き、観測度数が、パラメータのポワッソン分布から得られたとしている。すなわち、

である。

上のモデルに基づいてベイズ分析を行うStanスクリプトを、リスト１（poisson_reg.stan）に示す。ファイルは、poissonfiles.zipにまとめている。

リスト１　Stanスクリプト（poisson_reg.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];

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

}

リスト１では、２次元配列データである表１のデータを、ダミー変数を使った重回帰モデル（表２）で表している。

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

表２　表１のデータをダミー変数を使った重回帰モデルの形式で表現
データ番号
１	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

リスト１のStanスクリプトを用いて表１のデータをベイズ分析するPythonスクリプトを、リスト２（poisson_reg.py）のように用意した。ファイルはpoissonfiles.zipにまとめている。

リスト２　リスト１のStanスクリプトにより表１のデータを分析するPythonスクリプト（poisson_reg.py）

from cmdstanpy import CmdStanModel

import numpy as np

import scipy.stats as ss

import matplotlib.pyplot as plt

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}

model = CmdStanModel(stan_file="poisson_reg.stan") # コンパイル

fit = model.sample(data=Data) # サンプリング

print(fit.summary())

inf_data = az.from_cmdstanpy(fit) # Arviz(az)用のデータに変換

# トレースプロット

az.plot_trace(inf_data,

var_names=['mu','a2','a3','b2','b3','g22','g23','g32','g33'],

figsize=(10,8)) #, var_names=['b0', 'b1', 'sgm'], figsize=(10,7))

plt.tight_layout()

plt.show()

df_sample = fit.draws_pd() # PandasのDataFrame型に変換

"""

　　　　　　主効果パラメータの事後分布

"""

sb.kdeplot(df_sample['a2'], label='a2')

sb.kdeplot(df_sample['a3'], label='a3')

sb.kdeplot(df_sample['b2'], label='b2')

sb.kdeplot(df_sample['b3'], label='b3')

plt.xlabel('Main Effects', fontsize=14)

plt.legend()

plt.tight_layout()

plt.show()

"""

交互作用パラメータの事後分布

"""

sb.kdeplot(df_sample['g22'], label='g22')

sb.kdeplot(df_sample['g23'], label='g23')

sb.kdeplot(df_sample['g32'], label='g32')

sb.kdeplot(df_sample['g33'], label='g33')

plt.xlabel('Interaction', fontsize=14)

plt.legend()

plt.tight_layout()

plt.show()

a3 = df_sample['a3']

cnt_a3_pos = (a3 > 0.0).sum()

print('\nP(a3 > 0.0) = {0:.3f}'.format(cnt_a3_pos / len(a3)))

b3 = df_sample['b3']

cnt_b3_pos = (b3 > 0.0).sum()

print('P(b3 > 0.0) = {0:.3f}'.format(cnt_b3_pos / len(b3)))

g33 = df_sample['g33']

cnt_g33_neg = (g33 < 0.0).sum()

print('P(g33 < 0.0) = {0:.3f}'.format(cnt_g33_neg / len(g33)))

sb.kdeplot(a3, label=f'P(a3>0)={cnt_a3_pos/len(a3):.3f}')

sb.kdeplot(b3, label=f'P(b3>0)={cnt_b3_pos/len(b3):.3f}')

sb.kdeplot(g33, label=f'P(g33<0)={cnt_g33_neg/len(g33):.3f}')

plt.plot([0,0],[0.0,0.3], c='r', lw=2)

plt.legend()

plt.show()

次のコード

model = CmdStanModel(stan_file="poisson_reg.stan")

により、リスト１のStanスクリプトのコンパイルが行われ、次のコード

fit = model.sample(data=Data)

により、MCMCサンプリングが行われる。

サンプリングの統計は次のコードで表示される。

print(fit.summary())

トレース図をarviz（az）の関数で表示するために、サンプリングの結果fitをarvizのInferenceData型に次のコードで変換する。

inf_data = az.from_cmdstanpy(fit)

このinf_dataを用いて、次のスクリプトでトレース図を描いている。

az.plot_trace(inf_data,

var_names=['mu','a2','a3','b2','b3','g22','g23','g32','g33'],

figsize=(10,8)) #, var_names=['b0', 'b1', 'sgm'], figsize=(10,7))

plt.tight_layout()

plt.show()

MCMCサンプルの値を取り出すために、次のコードでサンプルをPandasのDataFrame型に変換している。

df_sample = fit.draws_pd()

パラメータの変数名は、DataFrameのキーを調べればわかる。パラメータ変数名を表す文字列、あるいは配列であれば配列の要素を表す文字列がキーになっている。

リスト２のスクリプトを、次のコマンドによりUbuntu24/WSL/Windowsで実行した。

(stan) *****/poissonfiles$ python poisson_reg.py

実行すると、コンパイルに多少に時間が掛かるが、コンパイル後、MCMCサンプリングが実行され、トレース図（図１）が表示される。

図１

図１のWindowの右上のX印アイコンをクリックして閉じると、次に主効果パラメータの事後分布のグラフが表示される（図２）。

図２

図２のWindowを閉じると、次に交互作用パラメータの事後分布が表示される（図３）。

図３

図２および図３の事後分布において、分布の主要な領域が０から離れていると思われる分布を取り出し、事後分布が０より大きい（正である）事後確率、あるいは事後分布が０より小さい（負である）事後分布を示すグラフが表示される（図４）。

図４

これらの事後分布では、正である領域、あるいは負である領域の事後確率が95％より大きいので、これらのパラメータ値は０でないと考えられる。

図４のWindowを閉じると、スクリプトの実行終了である。

Home