度数の分析
カイ二乗検定とベイズ分析
度数の検定において適合度検定とクロス表（Contingency Table）の検定のためのPythonスクリプト例を用意した。スクリプトファイルは圧縮ファイルFrequencies.zipにまとめた。Jupyter Notebook用にまとめたものも、圧縮ファイルFrequencies_jn.zipとして用意した。それぞれダウンロード・解凍して自由に利用されたい。
なお、ベイズ分析の例も上げておいた。
観測度数の分布の確率分布モデルへの適合度の検定のためのスクリプト例を作成した。
サイコロを50回投げて表A.1に示す結果を得たとする。
|
表A.1 サイコロ投げの結果 |
||||||
|
目 |
１の目 |
２の目 |
３の目 |
４の目 |
５の目 |
６の目 |
|
度数 |
9 |
11 |
9 |
6 |
8 |
7 |
表A.1に示されている度数が、等確率の分布に従うものであると見なせるかどうかの検定のために次のスクリプトを用意した。
import scipy.stats as ss
data = [ 9, 11, 9,
6, 8, 7]
chi2, p = ss.chisquare(data)
print('chi_squared = ', chi2)
print('p-value = ', p)
df = len(data) - 1
print('df = ', df)
度数データをリストdataに設定して、検定のための関数を
chi2, p = ss.chisquare(data)
と呼び出している。関数の戻り値chi2に
の値、pにｐ値が返される。
上のスクリプトを実行すると、以下の結果を得る。
chi_squared = 1.84
p-value = 0.870811402015331
df = 5
クロス表の分析のためのプログラムを作成した。
表B.1は京都大学学部入学者の平成31年度の男女別データであある。
|
表B.1 平成３１年度男女別合格者（京都大学ホームページより） 合計が200名以上の学部、内数で「認」と「他」を除く） |
|||||
|
|
文学部 |
法学部 |
理学部 |
工学部 |
農学部 |
|
男性 |
130 |
236 |
281 |
858 |
189 |
|
女性 |
83 |
75 |
32 |
79 |
103 |
表B.1のカイ二乗検定のためのスクリプトを以下のように作成した。
import numpy as np
import scipy.stats as ss
data = [[130, 236, 281, 858, 189],
[
83, 75, 32,
79, 103]]
chi2, p, df, expt =
ss.chi2_contingency(data)
print('chi-squared = ', chi2)
print('p-value = ', p)
print('df = ', df)
print('Expected <--> Data')
for e, d in zip(expt, data):
print(e, ' ', d)
data = np.array(data)
dif = data - expt
print('data - expected')
for v in dif:
for v1 in v:
print('
{0:<10.1f}'.format(v1), end = '')
print()
dif_r = dif / expt
print('(data - expected) / expected')
for v in dif_r:
for v1 in v:
print('
{0:<10.1f}'.format(v1), end = '')
print()
データをリストdataに設定して、カイ二乗検定のための関数を
chi2, p, df, expt =
ss.chi2_contingency(data)
と呼び出している。
戻り値は、chi2, p,
df, exptに値、ｐ値、自由度、期待度数が返される。
上のスクリプトを実行すると、以下の結果を得る（端末出力からコピー・貼り付け）。
= RESTART:
D:\yasuharu\www\LaCoocan\Python\IntroStat\Frequencies\FrequenciesFiles\TestCrossTable.py
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.3
-0.1
0.1
0.1
-0.2
1.2
0.3
-0.4
-0.5
1.0
>>>
ｐ値2.02134193518844e-42、数値
を表している。
期待度数とデータの度数は、左に期待度数、右にデータ値である。それに続いて、データから期待度数を引いたもの、それらを期待度数で割ったもの出力されている。
なお、クロス表作成の簡単なスクリプト例のウェブサイトも用意した。
表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章で、度数／比率データのロジスティック回帰モデルによる分析例を＜岡本安晴「心理学データ分析と測定」勁草書房＞の第４章で説明している。分析用プログラムは、それぞれのウェブサイトから自由にダウンロードできる。