KStwo 分布#

這是從 \(n\) 個樣本或觀測值計算出的經驗分布函數,與一個比較(或目標)累積分布函數(假設為連續函數)之間的最大絕對差異的分布。(名稱中的「two」是因為這是雙邊差異。ksone 是正差異的分布,\(D_n^+\),因此它關注單邊差異。kstwobign標準化最大絕對差異 \(\sqrt{n} D_n\) 的極限分布。)

寫成 \(D_n = \sup_t \left|F_{empirical,n}(t)-F_{target}(t)\right|\)kstwo\(D_n\) 值的分布。

kstwo 也可用於兩個經驗分布函數之間的差異,適用於分別具有 \(m\)\(n\) 個樣本的觀測值集合。寫成 \(D_{m,n} = \sup_t \left|F_{1,m}(t)-F_{2,n}(t)\right|\),其中 \(F_{1,m}\)\(F_{2,n}\) 是兩個經驗分布函數,則在適當條件下,\(Pr(D_{m,n} \le x) \approx Pr(D_N \le x)\),其中 \(N = \sqrt{\left(\frac{mn}{m+n}\right)}\)

有一個形狀參數 \(n\),為正整數,且 support 為 \(x\in\left[0,1\right]\)

此實作遵循 Simard & L’Ecuyer 的方法,結合 Durbin 和 Pomeranz 的精確演算法與 Li-Chien、Pelz 和 Good 的漸近估計,以計算具有 5-15 位有效數字的 CDF。

範例#

>>> import numpy as np
>>> from scipy.stats import kstwo

顯示大小為 5 的樣本,間隙至少為 0、0.5 和 1.0 的機率

>>> kstwo.sf([0, 0.5, 1.0], 5)
array([1.   , 0.112, 0.   ])

將從來源 N(0.5, 1) 分布中抽取的 5 個樣本,與目標 N(0, 1) CDF 進行比較。

>>> from scipy.stats import norm
>>> n = 5
>>> gendist = norm(0.5, 1)       # Normal distribution, mean 0.5, stddev 1
>>> x = np.sort(gendist.rvs(size=n, random_state=np.random.default_rng()))
>>> x
array([-1.59113056, -0.66335147,  0.54791569,  0.78009321,  1.27641365])  # may vary
>>> target = norm(0, 1)
>>> cdfs = target.cdf(x)
>>> cdfs
array([0.0557901 , 0.25355274, 0.7081251 , 0.78233199, 0.89909533])   # may vary
>>> # Construct the Empirical CDF and the K-S statistics (Dn+, Dn-, Dn)
>>> ecdfs = np.arange(n+1, dtype=float)/n
>>> cols = np.column_stack([x, ecdfs[1:], cdfs, cdfs - ecdfs[:n], ecdfs[1:] - cdfs])
>>> np.set_printoptions(precision=3)
>>> cols
array([[-1.591,  0.2  ,  0.056,  0.056,  0.144],     # may vary
       [-0.663,  0.4  ,  0.254,  0.054,  0.146],
       [ 0.548,  0.6  ,  0.708,  0.308, -0.108],
       [ 0.78 ,  0.8  ,  0.782,  0.182,  0.018],
       [ 1.276,  1.   ,  0.899,  0.099,  0.101]])
>>> gaps = cols[:, -2:]
>>> Dnpm = np.max(gaps, axis=0)
>>> Dn = np.max(Dnpm)
>>> iminus, iplus = np.argmax(gaps, axis=0)
>>> print('Dn- = %f (at x=%.2f)' % (Dnpm[0], x[iminus]))
Dn- = 0.246201 (at x=-0.14)
>>> print('Dn+ = %f (at x=%.2f)' % (Dnpm[1], x[iplus]))
Dn+ = 0.224726 (at x=0.19)
>>> print('Dn  = %f' % (Dn))
Dn  = 0.246201
>>> probs = kstwo.sf(Dn, n)
>>> print(chr(10).join(['For a sample of size %d drawn from a N(0, 1) distribution:' % n,
...      ' Kolmogorov-Smirnov 2-sided n=%d: Prob(Dn >= %f) = %.4f' % (n, Dn, probs)]))
For a sample of size 5 drawn from a N(0, 1) distribution:
 Kolmogorov-Smirnov 2-sided n=5: Prob(Dn >= 0.246201) = 0.8562

繪製經驗 CDF 與目標 N(0, 1) CDF 的圖

>>> import matplotlib.pyplot as plt
>>> plt.step(np.concatenate([[-3], x]), ecdfs, where='post', label='Empirical CDF')
>>> x3 = np.linspace(-3, 3, 100)
>>> plt.plot(x3, target.cdf(x3), label='CDF for N(0, 1)')
>>> plt.ylim([0, 1]); plt.grid(True); plt.legend();
>>> plt.vlines([x[iminus]], ecdfs[iminus], cdfs[iminus], color='r', linestyle='solid', lw=4)
>>> plt.vlines([x[iplus]], cdfs[iplus], ecdfs[iplus+1], color='m', linestyle='solid', lw=4)
>>> plt.annotate('Dn-', xy=(x[iminus], (ecdfs[iminus]+ cdfs[iminus])/2),
...              xytext=(x[iminus]+1, (ecdfs[iminus]+ cdfs[iminus])/2 - 0.02),
...              arrowprops=dict(facecolor='white', edgecolor='r', shrink=0.05), size=15, color='r');
>>> plt.annotate('Dn+', xy=(x[iplus], (ecdfs[iplus+1]+ cdfs[iplus])/2),
...             xytext=(x[iplus]-2, (ecdfs[iplus+1]+ cdfs[iplus])/2 - 0.02),
...             arrowprops=dict(facecolor='white', edgecolor='m', shrink=0.05), size=15, color='m');
>>> plt.show()

參考文獻#

  • 「Kolmogorov-Smirnov 檢定」,維基百科 https://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test

  • Durbin J. 「樣本分布函數位於兩條平行直線之間的機率。」Ann. Math. Statist., 39 (1968) 39, 398-411.

  • Pomeranz J. 「Kolmogorov-Smirnov 統計量在小樣本下的精確累積分布(演算法 487)。」Communications of the ACM, 17(12), (1974) 703-704.

  • Li-Chien, C. 「關於 A. N. Kolmogorov 統計量的精確分布及其漸近展開。」Acta Matematica Sinica, 6, (1956) 55-81.

  • Pelz W, Good IJ. 「逼近 Kolmogorov-Smirnov 單樣本統計量的下尾區域。」Journal of the Royal Statistical Society, Series B, (1976) 38(2), 152-156.

  • Simard, R., L’Ecuyer, P. 「計算雙邊 Kolmogorov-Smirnov 分布」,Journal of Statistical Software, Vol 39, (2011) 11.

實作: scipy.stats.kstwo