「2つの分布が一致している」を棄却できるかどうか。
はじめに
- 今回は説明可能性についての4回目となります。今回は、モデルの評価をするための1つである、コルモゴロフースミルノフ検定（Kolmogorov–Smirnov test）、いわるゆKS検定について説明します。
- モデルの評価とPoCを作成し説明していく為に必要な情報のまとめはこちら。
- KS検定についてはすでに多くの解説がなされているので簡単に下記にまとめました。
KS検定とは
- KS検定とは、
統計学における仮説検定の一種であり、有限個の標本に基づいて、二つの母集団の確率分布が異なるものであるかどうか、あるいは母集団の確率分布が帰無仮説で提示された分布と異なっているかどうかを調べるために用いられる。
二つの標本を比較する最も有効かつ一般的なノンパラメトリック手法の一つである。
wikiより
- 上記の通り2つの分布が異なっているかどうかを調べるために用いられますが、2つの分布が一致していることを証明することはできません。
- ただし、2つの分布が一致していることを棄却することがきます。
- クレジットスコアの例に当てはめると、「通常の債権の分布と延滞の債権の分布は一致しているとは言えない」となります。
少し回りくどい言い回しとなりますが、クレジットスコアにおいて通常の債権の分布と延滞の債権の分布が一致している場合は、とうてい役に立つものではなくなってしまいます。そのため、必ずと言っていいほど使われる指標の1つではありますが、「一致しているとは言えない」だけであり、それ以上でもそれ以下でもありません。
- 以前にジニ係数・AUCなどで述べたのと同様に、KS検定も順序性尺度を使って評価する手法の1つとなります。
理論と計算
-
それぞれの分布は以下の計算からなります。
F_n(x)=\frac{ \#\{ 1\leq i \leq n | y_i \leq x\}} {n} -
それぞれの分布におけるKS検定統計量は以下の計算となります。
\;\;\;\;\;\;\;\;\;\;\;\;\; D^+_n = sup(F_n(x)-F(x))\\ x\\ \;\;\;\;\;\;\;\;\;\;\;\;\; D^-_n = sup(F_n(x)-F(x))\\ x -
通常は、上記の差が最大値となる値を利用します。また、2つの分布の出現頻度の差を見ているのですが、それに関しては、後ほど図と一緒に説明したいと思います。
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; D_n=sup|F_n(x) - F(x)| = max(D^+_n, D^-_n)\\ x
Pythonでの利用
-
Pythonではscipyを利用すれば簡単に利用できます。以前と同じデータを利用して試してみましょう。 まずはサンプルモデルを作成します。サンプルはこちらから。
# Import Library import numpy as np import pandas as pd import matplotlib.pyplot as plt from matplotlib.patches import ArrowStyle from scipy import stats from sklearn.ensemble import RandomForestClassifierX_train = pd.read_csv('X_train.csv') y_train = pd.read_csv('y_train.csv')# Building sample model clf= RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini', max_depth=25, max_features='auto', max_leaf_nodes=None, min_impurity_decrease=0.0, min_impurity_split=None, min_samples_leaf=1, min_samples_split=15, min_weight_fraction_leaf=0.0, n_estimators=50, n_jobs=4, oob_score=False, random_state=0, verbose=0, warm_start=False) clf.fit(X_train, y_train) -
Predictionを取得しラベルを付与します。
pre=clf.predict_proba(X_train) pre_pd = pd.DataFrame(pre[:,0]) pre_pd['y'] = y_train print(pre_pd.head())
>out
0 y
0 0.891448 0
1 0.042604 1
2 0.467200 1
3 0.016613 1
4 0.989407 0
-
Good（y=0）とBad（y=1）の2つに分類し、statsに渡すと検定統計量とP値が簡単に取得できます。
pre_good=pre_pd[pre_pd.y==0] pre_bad=pre_pd[pre_pd.y==1] ks=stats.ks_2samp(pri_good[0].values,pri_bad[0].values) print (ks) print(ks[0]) print(ks[1])
>out
Ks_2sampResult(statistic=0.7843447416355096, pvalue=5.132844304914007e-115)
0.7843447416355096
5.132844304914007e-115
PoCでの利用
通常のKS検定
-
上記のGoodとBadに2分したものを利用し、グラフを作成してみましょう。
X軸にスコア（Probability）をY軸に累積出現率（累積の件数）をプロットします。そのために、それぞれのスコアを昇順で並びかえてscore_good, score_badとし、その累積件数の割合をcumulative_good, cumulative_badとします。# Score and Accumulate cumulative_good = np.add.accumulate(pri_good) / pri_good.sum() cumulative_good = cumulative_good[0] score_good = pri_good.sort_values(by=0) score_good = score_good[0] cumulative_bad = np.add.accumulate(pri_bad) / pri_bad.sum() cumulative_bad = cumulative_bad[0] score_bad = pri_bad.sort_values(by=0) score_bad = score_bad[0] -
グラフに出力してみましょう。
fig = plt.figure() ax = fig.add_subplot(111) ax.plot(score_good,cumulative_good, color='r',label='Good') ax.plot(score_bad ,cumulative_bad, color='b',label='Bad') arrow_style_list = [ ArrowStyle('<|-|>', head_length=1, head_width=1) ] point = { 'start': [0.5, 0.1], 'end': [0.5, 0.8] } ax.annotate('', xy=point['end'], xytext=point['start'], xycoords='data', arrowprops=dict(arrowstyle=ArrowStyle('<|-|>', head_length=0.5, head_width=0.5), connectionstyle='arc3', facecolor='y', edgecolor='y' ) ) ax.set_title('Kolmogorov–Smirnov test') ax.set_xlabel('Score(Probability)') ax.set_ylabel('Accumlate %') ax.legend() -
下図のようにGoodとBadがそれぞれのスコアでどのくらい出現しているかを表しています。KS統計検定量はYの値の差が最大値になる値のことで、一番出現割合の差が大きいときを評価します。さきほどのスコアでは0.7843...となります。
今回は、X軸にスコア（Goodになる確率＝Badにならない確率）をプロットしているので、確率が低いとき（0に近い所）にはBadが多数発生し、確率が高いとき（1に近い所）Goodが多数発生します。
この特性のため、中央値あたりはボリュームの影響を受けやすくなり、裾の方は影響が小さくなりやすくなります。
10分位数でのKS検定
- 上記の理由から、Goodを10分位数などで区切りその区切りの中でBadとどれだけ差があるかを見る方法もあります。
-
まずは10分位数のスコアを取得します。（便宜上最小値は0、最大値は1にしています）
# Decile decile = np.arange(0.0, 1.1, 0.1) score_good_decile = pd.DataFrame(score_good.quantile(decile)) score_bad_decile= pd.DataFrame(score_bad.quantile(decile)) score_good_decile.iloc[0,0] =0 score_good_decile.iloc[10,0] =1 print(score_good_decile)
>out
0
0.0 0.000000
0.1 0.598406
0.2 0.752527
0.3 0.810827
0.4 0.854746
0.5 0.894150
0.6 0.914538
0.7 0.935213
0.8 0.948343
0.9 0.974172
1.0 1.000000
-
10分位数のそれぞれの値までの累積出現率を求めます。少し値の変換の仕方が綺麗ではありませんが下記のようなコードになります。
good_df = pd.DataFrame() for i in range(0,11): xd1 = pri_good[pri_good[0] <= score_good_decile.iloc[i,0]] a = xd1[0].count() / score_good.count() tmp_se = pd.Series([a]) good_df = good_df.append( tmp_se, ignore_index=True ) good_df = good_df.rename(columns={0: 'good'}) bad_df = pd.DataFrame() for i in range(0,11): xd1 = pri_bad[pri_bad[0] <= score_good_decile.iloc[i,0]] a = xd1[0].count() / score_bad.count() tmp_se = pd.Series([a]) bad_df = bad_df.append( tmp_se, ignore_index=True ) bad_df = bad_df.rename(columns={0: 'bad'}) # Add difference total_df = good_df total_df['bad'] = bad_df['bad'] total_df['dif'] = bad_df['bad'] - good_df['good'] total_df さきほど作成したリストに差分を加えると以下のようになります。
| good | bad | dif | |
|---|---|---|---|
| 0 | 0.000000 | 0.000000 | 0.000000 |
| 1 | 0.100182 | 0.880117 | 0.779935 |
| 2 | 0.200364 | 0.947368 | 0.747004 |
| 3 | 0.300546 | 0.979532 | 0.678986 |
| 4 | 0.400729 | 0.991228 | 0.590499 |
| 5 | 0.500911 | 0.994152 | 0.493241 |
| 6 | 0.599271 | 1.000000 | 0.400729 |
| 7 | 0.699454 | 1.000000 | 0.300546 |
| 8 | 0.799636 | 1.000000 | 0.200364 |
| 9 | 0.899818 | 1.000000 | 0.100182 |
| 10 | 1.000000 | 1.000000 | 0.000000 |
-
ご覧の通り、下記のように一番差が大きくなるのは、Index=1の時で差が0.77934...
つまり、KS統計検定量（KS max）= 0.77934...になります。ks = bad_df - good_df ks_max=ks.max()[0] ks_index=ks.idxmax()[0] print(ks_max) print(ks_index)
>out
0.7799348097018503
1
-
最初と同様にstatsを利用しても同じ結果になります。
good_df2 = pre_good good_df2['score'] = 0 bad_df2 = pre_bad bad_df2['score'] = 0 for i in range(0,11): good_df2.loc[(good_df2[0] <=score_good_decile.iloc[i,0]) & (good_df2['score']==0) , 'score'] = i/10 for i in range(0,11): bad_df2.loc[(bad_df2[0] <=score_good_decile.iloc[i,0]) & (bad_df2['score']==0) , 'score'] = i/10 ks2=stats.ks_2samp(good_df2['score'].values,bad_df2['score'].values) print (ks2) print(ks2[0]) print(ks2[1])
>out
Ks_2sampResult(statistic=0.7799348097018503, pvalue=9.8922304125642e-114)
0.7799348097018503
9.8922304125642e-114
-
こちらをプロットすると下記のようになります。
fig = plt.figure() ax = fig.add_subplot(111) ax.plot(decile,good_df, color='r', marker='o') ax.plot(decile,bad_df, color='b', marker='o') arrow_style_list = [ ArrowStyle('<|-|>', head_length=1, head_width=1) ] point = { 'start': [x_start, y_start], 'end': [x_start, y_end] } ax.annotate('', xy=point['end'], xytext=point['start'], xycoords='data', arrowprops=dict(arrowstyle=ArrowStyle('<|-|>', head_length=0.5, head_width=0.5), connectionstyle='arc3', facecolor='y', edgecolor='y' ) ) ax.text(x_start + 0.05, y_end /2, 'KS max = ' + str(f'{ks_max:.4g}')) ax.set_xticks(decile) ax.set_title('Cumulative distribution of bads and goods per score decile') ax.set_xlabel('Score(decile)') ax.set_ylabel('Accumlate %') fig.show() -
Goodを基準に10分位数としているので、当然Goodは1対1の比例で累積割合が増加しています。
これに対してBadは第1分位の時点で大きく上昇し、KSもMAX値となっています。さきほどは中央あたりで最大と見えていたのが、スコアが低い時点、つまりBadの確率が高い時点で最大の差になっていることがわかります。
まとめ
- KS検定は古くから使われており、すでに計算ツールも揃っている手法なので今更まとめる事でもないかもしれませんが、今回は、あえて少し違う計算の仕方とグラフにしてみました。
普段は簡単に計算して終わりだったのですが、改めて文章にすることで理解が深まると思います。 - クレジットスコアは、2つの分布が一致してはいけないので、一致していないことを評価するための指標として利用しています。
- 順序性尺位なので、以前のジニ係数・AUCと同様に理解しやすい指標の1つになると思います。
- サンプルが綺麗に分類されているものなので結果がこのような形でしたが、実際にはもっと差が小さくなることが一般的となります。
参考文献
追記
-
今回始めてグラフに矢印の線を入れてみたのですが、なかなか苦労しました。
この部分（'ここ'）に文字を入れると線がずれてしまうため少々苦労しました。ax.annotate('ここ', xy=point['end'], xytext=point['start'], xycoords='data' 同じX座標にしても、文字あり・なしで下記のようになります。
結局解決策がみつからずにplt.textを利用し対処しました。