２項分布の比較／ベイズ分析

2項分布のベイズ分析

比率の比較

図１のExcelデータは、平成３１年度男女別合格者（京都大学ホームページより：合計が200名以上の学部、内数で「認」と「他」を除く）である。

図１

学部ごとの女子学生の比率に対して、２項分布を用いたベイズ分析を行ってみる。すなわち、学部kの男子学生数を、女子学生数をで表し、女子学生の比率をとしたとき

とおく。モデル（１）に基づいてベイズ分析を行うためのStanスクリプトを本ウェブサイト後半に示すリスト１のように用意した。本ウェブサイトのスクリプトファイルなどは圧縮ファイルfilesBinomials.zipにまとめたので、ダウンロード解凍すれば自由に利用できる。

リスト１のStanスクリプトを使って図１の形式のデータを分析するPythonスクリプトを、本ウェブサイトの後半に示すリスト２として用意した。

リスト２のスクリプトを実行すると、入力データファイル名の設定を求めるダイアログウィンドウが表示される（図２）。

図２

入力データファイルは、図１の形式で用意する。１行目にグループ名（図１の場合は、学部名）を設定し、１列目に２項分布のカテゴリ値（図１の場合は、「男性」と「女性」）を書く。２項分布における該当カテゴリ値（「成功」カテゴリあるいは「１」として説明されることが多い）を３行目に、非該当カテゴリ値（「失敗」カテゴリあるいは「０」として説明されることが多い）を２行目に設定する。各グループにおけるカテゴリ値の数はそれぞれの行・列のセルに設定する。

図２では、図１のExcelデータファイル名「Data.xlsx」が設定されている。設定後、「開く」ボタンをクリックすると、画像保存用ファイル名の設定を求めるダイアログウィンドウが表示される（図３）。

図３

ファイル拡張子が「.png」である任意のファイル名を設定する。ファイル名を設定して「保存」ボタンをクリックすると、次に、計算結果の出力用テキストファイル名の設定を求めるダイアログウィンドウが表示される（図４）。

図４

テキストファイル用ファイル名を設定して「保存」ボタンをクリックすると、入力データファイルが読み込まれて処理が始まる。Stanスクリプトのコンパイルとビルドに多少時間が掛かるが、その後、MCMCサンプリングが行われ、サンプリング終了後に事後分布が図５のように表示される。

図５

女子の比率に対応する２項分布の確率の事後分布は、工学部と理学部、法学部、農学部と文学部の３グループに分かれていることが表されている。

図３のウィンドウを閉じると、プログラムの実行終了となる。実行終了後、図４で設定した出力テキストファイルを開くと、以下のようになっている。Python IDLEのエディタで開こうとすると

図６のボックスが表示されることがある。

図６

このときは、「s-jis」を設定して「OK」ボタンをクリックすればよい（図７）。

図７

出力ファイルの内容

入力データファイル = D:/yasuharu/www/LaCoocan/Python/IntroPyStat/BysBinomials/filesBinomials/Data.xlsx

K = 5

['文学部', '法学部', '理学部', '工学部', '農学部']

男性, 女性

[[130 83]

[236 75]

[281 32]

[858 79]

[189 103]]

Inference for Stan model: anon_model_9bbe48d3e7a751424b61594e434755ef.

4 chains, each with iter=2000; warmup=1000; thin=1;

post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

p[1] 0.39 4.5e-4 0.03 0.33 0.37 0.39 0.41 0.46 5181 1.0

p[2] 0.24 3.4e-4 0.02 0.2 0.23 0.24 0.26 0.29 5038 1.0

p[3] 0.1 2.4e-4 0.02 0.07 0.09 0.1 0.12 0.14 5061 1.0

p[4] 0.09 1.3e-4 9.1e-3 0.07 0.08 0.08 0.09 0.1 4812 1.0

p[5] 0.35 4.2e-4 0.03 0.3 0.33 0.35 0.37 0.41 4680 1.0

lp__ -890.0 0.03 1.58 -894.0 -890.8 -889.7 -888.8 -887.9 2251 1.0

Samples were drawn using NUTS at Sun Jul 19 10:37:11 2020.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at

convergence, Rhat=1).

P(女性)

文学部:

平均 = 0.39105

中央値 = 0.39101

Q1 = 0.36830 Q3 = 0.41270

95%確信区間 = [0.32880, 0.45646]

法学部:

平均 = 0.24387

中央値 = 0.24333

Q1 = 0.22752 Q3 = 0.26007

95%確信区間 = [0.19880, 0.29129]

理学部:

平均 = 0.10436

中央値 = 0.10323

Q1 = 0.09233 Q3 = 0.11567

95%確信区間 = [0.07353, 0.14061]

工学部:

平均 = 0.08510

中央値 = 0.08487

Q1 = 0.07867 Q3 = 0.09071

95%確信区間 = [0.06829, 0.10378]

農学部:

平均 = 0.35361

中央値 = 0.35323

Q1 = 0.33429 Q3 = 0.37283

95%確信区間 = [0.29953, 0.41217]

リスト１　２項分による分析のStanスクリプト（BinModel.stan）

data {

int K;

int Data[K, 2];

}

parameters {

real<lower = 0.0, upper = 1.0> p[K];

}

model {

for (i in 1:K)

p[i] ~ uniform(0.0, 1.0);

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 pystan

import matplotlib.pyplot as plt

import seaborn as sb

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('入力データExcelファイル', types = [('Excel(*.xlsx)', '*.xlsx')])

foutnm = askoutfilename(msg = '画像保存ファイル', types = [('画像 (*.png)', '*.png')],

deftype = '.png')

ftxtoutnm = askoutfilename(msg = 'テキスト出力ファイル', 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}

sm = pystan.StanModel(file = 'BinModel.stan')

fit = sm.sampling(data = StanData, n_jobs = 1)

print(fit)

fout.write('\n{}\n'.format(fit))

# 日本語表記の設定

plt.rcParams['font.family'] = ['Yu Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']

plt.figure(figsize = (9, 5))

plt.title('pの事後分布', fontsize = 18)

for k in range(K):

sb.kdeplot(fit['p'][:, k], label = VarGrp[k])

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)

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['p'][:, k].mean()

p_025, Q1, Med, Q3, p_975 = np.percentile(fit['p'][:, k], [2.5, 25, 50, 75, 97.5])