Pythonを使って経験分布と正規分布をグラフで視覚化してKS検定をしてみた
得られたデータが正規分布などの確率分布に適合しているかどうかを調べるのにKS検定は便利です。
このKS検定はExcelでも比較的簡単にできますが、Pythonのscify.statsに入っているkstestモジュールを使えば実質1行で検定できてしまいます。
実務ではこれでよいのでしょうが、KS検定の原理を知らずにブラックボックスにしておくのはよくありません。
そこで本記事ではあえてkstestモジュールを使わずに、自作関数を作ってKS検定を行ってみました。
そのためには経験分布を理解する必要がありますので、まずは経験分布関数をPythonを使ってグラフにすることから始めました。
データは【Excelで簡単！】KS検定で出荷数が正規分布に従うかどうかを検定するで扱ったのと同じデータを使います。
Pandasで読み込んだあと、Item_Aという名でnumpy.arrayに変換しています。
経験分布をグラフにする関数を作る
【Excelで簡単！】KS検定で出荷数が正規分布に従うかどうかを検定するで解説したように、KS検定では観測データを元に経験分布を作り、それと比較したい確率分布を比較します。
ですのでまずは経験分布を作るのが第一歩ですが、これは簡単です。
観測データを昇順に並べ変えた上で、1から順に通し番号を振り、それを全体のデータ数で割ればよいからです。
なぜこれでいいのかというと、通し番号がそれに対応する観測データ以下のデータ数になるからです。
これをPythonの関数にするのであれば、次のようにできます。
def ed(data_list, element):
n = sum(data_list <= element) # 各要素についてそれ以下の観測点数を調べる
n = n / len(data_list) # 経験分布に変換する
return element, nこの関数edに引数として観測データが入ったリストと変数elementを渡してやれば、変数element以下の観測データ数nを全体の観測データ数で割った値が、変数elementと一緒に返されます。
そしてこの関数の引数としてItem_Aの観測データリストとその要素を一つずつfor文で入れてやれば、各観測データとそれに対応する経験分布の値がリストaに格納されます。
Item_A.sort() # 観測データを昇順に並べ替え（順序統計量を作る）
a = [] # 順序統計量と経験分布を入れる空のリスト
for x in Item_A:
empirical_dist = ed(Item_A, x)
a.append(empirical_dist)
あとは、リストaに入っている観測データ（昇順に並べ替えたデータを順序統計量と呼びます）をリストxに、それに対応する経験分布の値をリストyに入れて、それをmatplotlibでグラフにします。
x = [] # 順序統計量を入れる空のリスト
y = [] # 経験分布を入れる空のリスト
cnt = 0 #
for i in a:
x.append(a[cnt][0])
y.append(a[cnt][1])
cnt += 1
plt.plot(x, y)
plt.showすると、次のように経験分布がグラフで視覚化できます。
正規分布をグラフにする関数を作る
次にItem_Aの観測データを近似する正規分布のグラフをmatplotlibで描いてみましょう。
まずは、観測データと変数elementを引数として渡したら、変数elementとそれに対応する正規分布の累積密度関数の値を返す関数sdを定義します。
def sd(data_list, element):
mean, std = data_list.mean(), data_list.std()
pdf = norm.cdf(element, mean, std)
return element, pdfそして先程の経験分布と同じようにItem_Aの各要素をfor文で順に関数sdに引数として渡すことで、順序統計量とそれに対応する正規分布の累積密度関数の値のリストを作ることができます。
Item_A.sort() # 観測データを昇順に並べ替え（順序統計量を作る）
b = [] # 順序統計量と正規分布を入れる空のリスト
for x in Item_A:
standard_dist = sd(Item_A, x)
b.append(standard_dist)そしてこれをmatplotlibでグラフにします。
x2 = [] # 順序統計量を入れる空のリスト
y2 = [] # 正規分布を入れる空のリスト
cnt = 0
for i in b:
x2.append(b[cnt][0])
y2.append(b[cnt][1])
cnt += 1
plt.plot(x2, y2)
plt.showすると次のような近似正規分布の累積分布関数がグラフで表示されます。
2つのグラフを重ねてKS検定を理解する
あとはPythonで「ヒストグラム＋正規分布の近似曲線」を作成する４通りの方法でやったように、いままでのコードをつなげて書けば2つのグラフを重ねて表示することができます。
# 経験分布を計算する関数を定義
def ed(data_list, element):
n = sum(data_list <= element) # 各要素についてそれ以下の観測点数を調べる
n = n / len(data_list) # 経験分布に変換する
return element, n
# 正規分布を計算する関数を定義
def sd(data_list, element):
mean, std = data_list.mean(), data_list.std()
pdf = norm.cdf(element, mean, std)
return element, pdf
Item_A.sort() # 観測データを昇順に並べ替え（順序統計量を作る）
a = [] # 順序統計量と経験分布を入れる空のリスト
for x in Item_A:
empirical_dist = ed(Item_A, x)
a.append(empirical_dist)
x = [] # 順序統計量を入れる空のリスト
y = [] # 経験分布を入れる空のリスト
cnt = 0
for i in a:
x.append(a[cnt][0])
y.append(a[cnt][1])
cnt += 1
plt.plot(x, y)
b = [] # 順序統計量と正規分布を入れる空のリスト
for x in Item_A:
standard_dist = sd(Item_A, x)
b.append(standard_dist)
x2 = [] # 順序統計量を入れる空のリスト
y2 = [] # 正規分布を入れる空のリスト
cnt = 0
for i in b:
x2.append(b[cnt][0])
y2.append(b[cnt][1])
cnt += 1
plt.plot(x2, y2)
plt.show
KS検定を行う関数を作る
KS検定は、経験分布と調べたい確率分布との最大乖離値を調べます。
これは次のような関数ks_normを作れば求めることができます。
Item_A.sort() # 観測データを昇順に並べ替え（順序統計量を作る）
def ks_norm(test_data):
difference = [] # 正規分布と経験分布の差異の絶対値を入れる空のリスト
for x in test_data:
normal_dist = sd(test_data, x) # 正規分布
empirical_dist = ed(test_data, x) # 経験分布
difference.append(abs(normal_dist - empirical_dist))
ks_value = max(difference) # 差異の絶対値の最大値＝KS値
return ks_value実行すると、最大乖離値ks_valueは0.08486になります。
これをKS値と呼びます。
次にKS分布上でこれに対応するp値を求めます。
KS分布は観測データの数によって一意に決まり、観測データ数が91の場合には次のようになります。
x = np.linspace(0, 0.2, 100)
y = st.kstwo.sf(x, 91)
plt.plot(x, y)
plt.show
横軸で先程求めたKS値0.08486を探して、それに対応する縦軸の値を見ると大体0.5くらいになることがわかりますが、詳しい値まではわかりません。
これを求める関数もScifyに用意されており、st.kstwo.sf(KS値, データ数)で求められます。
従って、KS値とそれに対応するp値を次のように求めることができます。
Item_A.sort() # 観測データを昇順に並べ替え（順序統計量を作る）
# KS値とp値を求める関数を定義
def ks_norm(test_data):
difference = [] # 正規分布と経験分布の差異の絶対値を入れる空のリスト
for x in test_data:
normal_dist = sd(test_data, x) # 正規分布
empirical_dist = ed(test_data, x) # 経験分布
difference.append(abs(normal_dist - empirical_dist))
ks_value = max(difference) # 差異の絶対値の最大値＝KS値
p_value = st.kstwo.sf(ks_value, len(test_data)) # KS値に対応するp値を求める
return [ks_value, p_value]
ks_norm(Item_A)
[0.08486080924359743, 0.5018970577950879]つまりp値は0.5019で遥かに0.05よりも大きいため、観測データが正規分布に従うという帰無仮説を棄却できません。
これはPythonのkstestを使って検定した結果とも、Excelを使って検定した結果とも一致します。