度数の分析
カイ二乗検定とベイズ分析
度数の検定において適合度検定とクロス表（Contingency Table）の検定のためのPythonスクリプト例を用意した。スクリプトファイルは圧縮ファイルFrequencies.zipにまとめた。Jupyter Notebook用にまとめたものも、圧縮ファイルFrequencies_jn.zipとして用意した。それぞれダウンロード・解凍して自由に利用されたい。
なお、ベイズ分析の例も上げておいた。
観測度数の分布の確率分布モデル（一様分布）への適合度の検定のためのスクリプト例を作成した。
サイコロを50回投げて表A.1に示す結果を得たとする。
|
表A.1 サイコロ投げの結果 |
||||||
|
目 |
１の目 |
２の目 |
３の目 |
４の目 |
５の目 |
６の目 |
|
度数 |
9 |
11 |
9 |
6 |
8 |
7 |
表A.1のデータを図A.1に示すCSVファイルとした。
図A.1
図A.1のCSVファイルTestUniFreq.csvの形式は、１行目にカテゴリ名を書き、その下に度数を書くものである。
図A.1のファイルを入力として適合度検定を行うスクリプトをリストA.1のように用意した。
リストA.1のスクリプトを実行すると、入力データファイル名の設定が求められる。
例として、図A.1のファイルTestUniFreq.csvを設定する。
File Name =
TestUniFreq.csv
ファイル名を設定すると、ファイルが読み込まれて、図A.2のグラフが表示される。
図A.2
グラフのタイトル部にカイ２乗の値、自由度、ｐ値が表示されている。
グラフには、一様分布のときの確率が赤の直線で、データの度数の相対頻度が青の折れ線で表示されている。
グラフのWindowを閉じると、端末には以下の出力が表示されている。
cat_names =
['１の目' '２の目' '３の目' '４の目' '５の目' '６の目']
data =
[ 9 11 9
6 8 7]
chi_squared = 1.84
p-value = 0.870811402015331
df = 5
FigFit.png was saved.
import pandas as pd
import numpy as np
import scipy.stats as ss
import matplotlib. pyplot as plt
import japanize_matplotlib
fin_nm = input('File Name = ')
d_frm = pd.read_csv(fin_nm,
encoding='shift-jis')
cat_names = d_frm.keys().values
data = d_frm.values[0]
print('cat_names =\n', cat_names)
print('data =\n', data)
chi2, p = ss.chisquare(data)
print('chi_squared = ', chi2)
print('p-value = ', p)
df = len(data) - 1
print('df = ', df)
N = np.sum(data)
n_cat = len(cat_names)
p_uni = np.ones(n_cat) / n_cat
plt.title(f'$\\chi^2$={chi2:.3f}, df={df}, p-value={p:.5f}')
plt.plot(cat_names, data / N,
label='Observed')
plt.plot(cat_names, np.ones(n_cat) *
p_uni, label='Expected')
plt.xticks(rotation=45)
plt.yticks([0, 0.5, 1])
plt.ylabel('Probability')
plt.legend()
plt.tight_layout()
plt.savefig('FigFit.png')
plt.show()
print('FigFit.png was saved.')
２×２のクロス表のFisher’s Exact検定については、このウェブサイトで扱っている。
表B.1は京都大学学部入学者の平成31年度の男女別データであある。
|
表B.1 平成３１年度男女別合格者（京都大学ホームページより） 合計が200名以上の学部、内数で「認」と「他」を除く） |
|||||
|
|
文学部 |
法学部 |
理学部 |
工学部 |
農学部 |
|
男性 |
130 |
236 |
281 |
858 |
189 |
|
女性 |
83 |
75 |
32 |
79 |
103 |
表B.1のデータを、図B.1に示すCSVファイルDataCrossTable.csvとした。
図B.1
１行目２列目から列項目のカテゴリ名を書き、１列目２行目から行項目のカテゴリ名を書いている。列カテゴリと行カテゴリの組み合わせのセルに度数を書いている。
図B.1の形式のCSVファイルを読み込んで分析するスクリプトをリストB.1のように作成した。
リストB.1のスクリプトを実行すると、入力データファイル名の設定が求められる。
例として、図B.1のファイル名DataCrossTable.csvを設定する。
File
Name(*.csv) = DataCrossTable.csv
ファイル名を設定すると、ファイルが読み込まれて、図B.2のグラフが表示される。
図B.2
グラフのタイトル部に、カイ2乗の値、自由度、ｐ値が表示されている。
グラフは、
を表し、カイ2乗値を計算するときの和に現れる各項の平方根をデータ値と期待値の差の符号に合わせたものである。クロス表の検定対象である交互作用が視覚的に表されている。
図B.2のWindowを閉じると、プログラムの実行終了となる。端末画面には、以下の出力が表示されている。
File Name(*.csv) = DataCrossTable.csv
row_cat =
['男性' '女性']
col_cat =
['文学部' '法学部' '理学部' '工学部' '農学部']
data =
[[130 236 281 858 189]
[ 83
75 32 79 103]]
chi-squared = 201.2522256170101
p-value = 2.02134193518844e-42
df = 4
Expected <--> Data =
[174.64762827 255.00193611
256.64181994 768.28557599 239.42303969] [130 236 281 858 189]
[ 38.35237173 55.99806389 56.35818006 168.71442401 52.57696031] [ 83 75
32 79 103]
data - expected =
-44.6
-19.0
24.4
89.7
-50.4
44.6
19.0
-24.4
-89.7
50.4
(data - expected) / expected**0.5 =
-0.3
-0.1 0.1
0.1
-0.2
1.2
0.3
-0.4
-0.5
1.0
FigCrossTable.png was saved.
import pandas as pd
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import japanize_matplotlib
fin_nm = input('File Name(*.csv) = ')
d_frm = pd.read_csv(fin_nm,
encoding='shift-jis')
row_cat = d_frm.values[:,0]
col_cat = d_frm.keys().values[1:]
print('row_cat =\n', row_cat)
print('col_cat =\n', col_cat)
data = np.array(d_frm.values[:,1:],
dtype=int)
print('data =\n', data)
chi2, p, df, expt =
ss.chi2_contingency(data)
print('chi-squared = ', chi2)
print('p-value = ', p)
print('df = ', df)
print('Expected <--> Data =\n')
for e, d in zip(expt, data):
print(e, ' ', d)
data = np.array(data)
dif = data - expt
print('data - expected =\n')
for v in dif:
for v1 in v:
print('
{0:<10.1f}'.format(v1), end = '')
print()
dif_r = dif / expt
dif_r0p5 = dif / (expt**0.5)
print('(data - expected) /
expected**0.5 =\n')
for v in dif_r:
for v1 in v:
print('
{0:<10.1f}'.format(v1), end = '')
print()
for i in range(len(row_cat)):
plt.plot(col_cat,
dif_r0p5[i], label=row_cat[i])
plt.ylabel('$(O-E)/\\sqrt{E}$')
plt.title(f'$\\chi^2$={chi2:.3f}, df={df}, p={p:.5f}')
plt.plot(col_cat,
np.zeros(len(col_cat)), '--', c='k', label='O-E=0')
plt.legend()
plt.tight_layout()
plt.savefig('FigCrossTable.png')
plt.show()
print('FigCrossTable.png was saved.')
なお、クロス表作成の簡単なスクリプト例のウェブサイトも用意した。
表B.1のデータに各学部の女子学生の確率を２項分布を適用してベイズ推定してみた。
Stanスクリプトは、以下のようである。
data {
int bun_m;
int bun_f;
int hou_m;
int hou_f;
int ri_m;
int ri_f;
int kou_m;
int kou_f;
int nou_m;
int nou_f;
}
parameters {
real<lower = 0.0, upper =
1.0> bun_p;
real<lower = 0.0, upper =
1.0> hou_p;
real<lower = 0.0, upper =
1.0> ri_p;
real<lower = 0.0, upper =
1.0> kou_p;
real<lower = 0.0, upper =
1.0> nou_p;
}
model {
bun_f ~ binomial(bun_f +
bun_m, bun_p);
hou_f ~ binomial(hou_f +
hou_m, hou_p);
ri_f ~ binomial(ri_f + ri_m, ri_p);
kou_f ~ binomial(kou_f +
kou_m, kou_p);
nou_f ~ binomial(nou_f +
nou_m, nou_p);
}
このStanスクリプト（FreqBin.stan）を用いて、以下のPythonスクリプトで分析した。
import matplotlib.pyplot as plt
import japanize_matplotlib
import pystan
import seaborn as sb
Data = {'bun_m': 130, 'bun_f': 83,
'hou_m': 236, 'hou_f': 75,
'ri_m': 281, 'ri_f': 32,
'kou_m': 858, 'kou_f': 79,
'nou_m': 189, 'nou_f': 103}
sm = pystan.StanModel(file =
'FreqBin.stan')
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
sb.kdeplot(fit['bun_p'], label = '文学部')
sb.kdeplot(fit['hou_p'], label = '法学部')
sb.kdeplot(fit['ri_p'], label = '理学部')
sb.kdeplot(fit['kou_p'], label = '工学部')
sb.kdeplot(fit['nou_p'], label = '農学部')
plt.xlabel('P(女子学生)', fontsize = 14)
plt.yticks([])
plt.legend(fontsize = 14)
plt.show()
スクリプトファイルは圧縮ファイルFreqBinomial.zipにまとめたので、ダウンロード解凍すれば自由に利用できる。
モジュールjapanize_matplotlibは、グラフ表示で日本語が使えるようにするものである。Anacondaのcondaコマンドではインストールできなかった（2020.06.14）が、pipコマンドでインストールできた。
上のスクリプトを実行すると下図のグラフが表示される。
文学部と農学部、法学部、工学部と理学部の３グループに分かれていることが表されている。
＊＊＊＊＊
度数データのをポアッソン回帰モデルによる分析例を＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の第12章、第15章で、度数／比率データのロジスティック回帰モデルによる分析例を＜岡本安晴「心理学データ分析と測定」勁草書房＞の第４章で説明している。分析用プログラムは、それぞれのウェブサイトから自由にダウンロードできる。