2項分布の階層的ベイズ分析
比率の比較
2025.12：CmdStanPy対応に改訂
図１のExcelデータは、平成３１年度男女別合格者（京都大学ホームページより：合計が200名以上の学部、内数で「認」と「他」を除く）である。
図１
学部ごとの女子学生の比率に対して、２項分布を用いたベイズ分析を行ってみる。すなわち、学部kの男子学生数を
、女子学生数を
で表し、女子学生の比率を
としたとき
とおく。さらに、確率は、当該の大学における女子の比率として共通の母集団からのサンプルとして考える。
ここで、
と
は、ハイパーパラメータである（Lambert, 2018; Kruschke,
2015）。
なお、図１はデータの例であり、学部を表している列は学部以外のカテゴリで２カテゴリ以上であれば、図１の形式でExcelデータを用意すればよい。行カテゴリも性別以外のカテゴリでよいが、２項分布を適用しているので２カテゴリであればよい。
モデル（１）と（２）に基づいてベイズ分析を行うためのStanスクリプトを本ウェブサイト後半に示すリスト１のように用意した。本ウェブサイトのスクリプトファイルなどは圧縮ファイルfilesBinomials.zipにまとめたので、ダウンロード解凍すれば自由に利用できる。
リスト１のStanスクリプトを使って図１の形式のデータを分析するPythonスクリプトを、本ウェブサイトの後半に示すリスト２として用意した。
リスト１のStanスクリプトファイル（BinModel.stan）、リスト２のPythonスクリプトファイル（MainBinomials.py）、および入力データファイル（例えば、Data.xlsx）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、CmdStanPyのインストールされた環境で次のコマンドを実行する。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
(stan) *****/filesBinomials$ python MainBinomials.py
上のコマンドを実行すると、入力データファイル名の設定を求めるダイアログウィンドウが表示される（図２）。
図２
入力データファイルは、図１の形式で用意する。１行目にグループ名（図１の場合は、学部名）を設定し、１列目に２項分布のカテゴリ値（図１の場合は、「男性」と「女性」）を書く。２項分布における該当カテゴリ値（「成功」カテゴリあるいは「１」として説明されることが多い）を３行目に、非該当カテゴリ値（「失敗」カテゴリあるいは「０」として説明されることが多い）を２行目に設定する。各グループにおけるカテゴリ値の数はそれぞれの行・列のセルに設定する。
図２では、図１のExcelデータファイル名「Data.xlsx」が設定されている。設定後、「Open」ボタンをクリックすると、画像保存用ファイル名の設定を求めるダイアログウィンドウが表示される（図３）。
図３
ファイル拡張子が「.png」である任意のファイル名を設定する。グラフは、設定したファイル名の後ろをグラフに合わせた名前に変更したファイル名で保存される。
ファイル名を設定して「Save」ボタンをクリックすると、次に、計算結果の出力用テキストファイル名の設定を求めるダイアログウィンドウが表示される（図４）。
図４
テキストファイル用ファイル名を設定して「save」ボタンをクリックすると、入力データファイルが読み込まれて処理が始まる。Stanスクリプトのコンパイルに多少時間が掛かるが、コンパイル後、MCMCサンプリングが始まる。MCMCサンプリングが終了すると、トレース図が表示される（図５）。
図５
左側に事後分布のグラフ、右側にサンプリングのトレース図が表示されている。
図５のWindowを閉じると、図６のグラフが表示される。
図６
式（２）のベータ分布のグラフの事後分布である。MCMCサンプルからaとbのサンプルをランダムに100対選び、各対をパラメータとするベータ分布のグラフを描いたものである。赤の太い実線で描かれた曲線は、それら100本の曲線を平均したものである。
スクリプトを、以下に示す。
idx_slct = rng.choice(len(fit['a']), size=100, replace=False)
xcoord = np.linspace(0,
1, 1000)
curve_list = []
for idx in idx_slct:
ycoord
= ss.beta.pdf(xcoord, fit['a'][idx],
fit['b'][idx])
plt.plot(xcoord,ycoord, alpha=0.3)
curve_list.append(ycoord)
mn_curve = np.mean(curve_list, axis=0)
plt.plot(xcoord, mn_curve, c='r', label='Mean values')
図６の分布は、0.1から0.4の範囲に集まっているので、学部全体の傾向として、式（２）で与えられる確率（女子学生の比率）は男子よりも小さいと言えそうである。
図６のWindowを閉じると、図７のグラフが表示される。
図７
学部別の女子の比率の事後分布のグラフである。女子の比率に対応する２項分布の確率の事後分布は、工学部と理学部、法学部、農学部と文学部の３グループに分かれていることが表されている。
図７のWindowを閉じると、スクリプトの実行終了である。
実行終了後、図４で設定した出力ファイルを開くと以下のようになっている。
入力データファイル = /mnt/c/WorkSpace/www/LaCoocan/Python/IntroPyStat/BysBinomials/filesBinomials/Data.xlsx
K = 5
['文学部', '法学部', '理学部', '工学部', '農学部']
男性, 女性
[[130 83]
[236
75]
[281
32]
[858
79]
[189 103]]
P(女性)
文学部:
平均 = 0.37744
中央値 = 0.37635
Q1 = 0.35551 Q3 = 0.39894
95%確信区間 = [0.31461, 0.44447]
法学部:
平均 = 0.24023
中央値 = 0.23993
Q1 = 0.22343 Q3 = 0.25655
95%確信区間 = [0.19445, 0.28822]
理学部:
平均 = 0.10901
中央値 = 0.10856
Q1 = 0.09709 Q3 = 0.12004
95%確信区間 = [0.07655, 0.14512]
工学部:
平均 = 0.08691
中央値 = 0.08656
Q1 = 0.08062 Q3 = 0.09292
95%確信区間 = [0.06948, 0.10654]
農学部:
平均 = 0.34543
中央値 = 0.34487
Q1 = 0.32580 Q3 = 0.36416
95%確信区間 = [0.29468, 0.40066]
Kruschke, J.K. (2015) Doing Bayesian data analysis. Academic Press
Lambert, B. (2018) A student’s guide to Bayesian statistics. SAGE.
リスト１ ２項分による分析のStanスクリプト（BinModel.stan）
data {
int K;
array[K, 2] int Data;
}
parameters {
real<lower=0.0001,
upper=1000> a;
// Hyper parameter
real<lower=0.0001,
upper=1000> b;
// Hyper parameter
array[K] real<lower =
0.0, upper = 1.0> p;
}
model {
a ~ uniform(0.0001,
1000); // Hyper distribution
b ~ uniform(0.0001,
1000); // Hyper distribution
for (i
in 1:K)
p[i] ~ beta(a, b);
for (i
in 1:K)
Data[i][2] ~ binomial(Data[i][1]
+ Data[i][2], p[i]);
}
リスト２ リスト１のStanスクリプトによるベイズ分析（MainBinomial.py）
import tkinter
as tk
import tkinter.filedialog
as fd
import pandas as pd
import numpy
as np
import scipy.stats
as ss
from cmdstanpy
import CmdStanModel
import matplotlib.pyplot
as plt
import japanize_matplotlib
import seaborn as sb
import arviz
as az
def askinfilename(msg
= None, types = [('', '*.*')]):
"""
入力用ファイル名の設定
"""
rt = tk.Tk()
rt.withdraw()
flnm
= fd.askopenfilename(title = msg, filetypes = types)
rt.destroy()
return flnm
def askoutfilename(msg
= None, types = [('', '*.*')], deftype = '.txt'):
"""
出力用ファイル名の設定
"""
rt = tk.Tk()
rt.withdraw()
flnm
= fd.asksaveasfilename(title = msg, filetypes =
types,
defaultextension = deftype)
rt.destroy()
return flnm
finnm = askinfilename('Input
Data File', types = [('Excel(*.xlsx)', '*.xlsx')])
foutnm = askoutfilename(msg
= 'Output Image File',
types = [('Image(*.png)', '*.png')], deftype = '.png')
ftxtoutnm = askoutfilename(msg
= 'Output Text File',
types = [('Text (*.txt)', '*.txt')],
deftype = '.txt')
fout = open(ftxtoutnm,
'w')
fout.write('入力データファイル = {}\n'.format(finnm))
xlsx_fl = pd.ExcelFile(finnm)
data_xlsx = pd.read_excel(xlsx_fl, header
= None)
data = data_xlsx.values
K = len(data[0])
- 1
print('K = ', K)
fout.write('\nK = {}\n'.format(K))
VarGrp = []
for i in
range(1, K + 1):
VarGrp.append(data[0][i])
print(VarGrp)
fout.write('\n{}\n'.format(VarGrp))
VarCnd1 = data[1][0]
VarCnd2 = data[2][0]
print(VarCnd1, ', ', VarCnd2)
fout.write('\n{0:}, {1:}\n'.format(VarCnd1,
VarCnd2))
Data = np.empty((K,
2), dtype = int)
for k in range(K):
for i
in range(2):
Data[k][i] = data[i
+ 1][k + 1]
print(Data)
fout.write('\n{}\n'.format(Data))
StanData = {'K': K, 'Data': Data}
model = CmdStanModel(stan_file='BinModel.stan')
fit = model.sample(data=StanData)
print(fit.diagnose())
print(fit.summary())
infdata = az.from_cmdstanpy(fit) # Arviz
InfereceData型へ変換
az.plot_trace(infdata)
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_trace.png')
plt.show()
fit = fit.draws_pd()
# Pandas DataFrame型への変換
print(fit.keys())
rng = np.random.default_rng()
idx_slct = rng.choice(len(fit['a']), size=100, replace=False)
xcoord = np.linspace(0,
1, 1000)
curve_list = []
for idx in idx_slct:
ycoord
= ss.beta.pdf(xcoord, fit['a'][idx],
fit['b'][idx])
plt.plot(xcoord,ycoord, alpha=0.3)
curve_list.append(ycoord)
mn_curve = np.mean(curve_list, axis=0)
plt.plot(xcoord, mn_curve, c='r', label='Mean values')
plt.legend()
plt.title('Beta Distributions', fontsize=16)
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_beta.png')
plt.show()
plt.figure(figsize =
(9, 5))
plt.title('pの事後分布', fontsize = 18)
for k in range(K):
med = np.median(fit[f'p[{k+1}]'])
sb.kdeplot(fit[f'p[{k+1}]'], label = VarGrp[k] +
f'\nMed.={med:.3f}')
plt.yticks([])
plt.xlabel('p', fontsize
= 14)
plt.legend(loc = 'upper right', bbox_to_anchor = (1.25, 1.0), fontsize
= 12)
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_p.png')
plt.show()
print('\n\nP({})\n'.format(VarCnd2))
fout.write('\n\n\nP({})\n\n'.format(VarCnd2))
for k in range(K):
print('\n{}:'.format(VarGrp[k]))
fout.write('\n{}:\n'.format(VarGrp[k]))
p_mean
= fit[f'p[{k+1}]'].mean()
p_025, Q1, Med, Q3, p_975 = np.percentile(fit[f'p[{k+1}]'],
[2.5, 25, 50, 75, 97.5])
print('平均 = {0:.5f}'.format(p_mean))
fout.write('平均 = {0:.5f}\n'.format(p_mean))
print('中央値 = {0:.5f}'.format(Med))
fout.write('中央値 = {0:.5f}\n'.format(Med))
print('Q1 = {0:.5f} Q3 =
{1:.5f}'.format(Q1, Q3))
fout.write('Q1
= {0:.5f} Q3 =
{1:.5f}\n'.format(Q1, Q3))
print('95%確信区間 = [{0:.5f},
{1:.5f}]'.format(p_025, p_975))
fout.write('95%確信区間 = [{0:.5f}, {1:.5f}]\n'.format(p_025, p_975))
fout.close()