Kolmogorov-Smirnov 이표본 검정¶

from scipy import stats
A = [10.1, 10.3, 10.2, 10.5, 10.4]
B = [10.0, 10.2, 10.6, 10.8, 10.4]
print(stats.ks_2samp(A, B, method='exact').pvalue)   # 0.873016
print(stats.ks_2samp(A, B, method='asymp').pvalue)   # 0.820000

0.873015873015873
0.82

import itertools
import math
from collections import Counter

import numpy as np

n1 = n2 = 5
N = n1 + n2
pool = np.array(sorted(A + B))
print(f"정렬한 합친 표본 {pool.tolist()}")
print(f"동점 {[(k, v) for k, v in sorted(Counter(pool.tolist()).items()) if v > 1]}")

# 252 가지 배치를 전수 열거한다.
Ds = []
for comb in itertools.combinations(range(N), n1):
    step = np.full(N, -1 / n2)
    step[list(comb)] = 1 / n1
    Ds.append(np.abs(np.cumsum(step)).max())
Ds = np.round(Ds, 10)
tally = Counter(Ds)
print(f"\n열거한 배치 수 {len(Ds)}  (C(10,5) = {math.comb(N, n1)})")
print(f"D 가 가질 수 있는 값 {sorted(tally)}")
print(f"배치 수           {[tally[d] for d in sorted(tally)]}")

D_obs = 0.4
hit = int((Ds >= D_obs - 1e-12).sum())
print(f"\nD >= {D_obs} 인 배치 {hit} 개  ->  정확 p = {hit}/{len(Ds)}"
      f" = {hit / len(Ds):.9f}")
print(f"scipy method='exact' : {stats.ks_2samp(A, B, method='exact').pvalue:.9f}")
print("이 검정이 낼 수 있는 p-값 전부:")
for d in sorted(tally):
    k = int((Ds >= d - 1e-12).sum())
    print(f"  D = {d}:  p = {k}/{len(Ds)} = {k / len(Ds):.4f}")

c = np.sqrt(n1 * n2 / N) * D_obs
series = 2 * sum((-1) ** (k - 1) * np.exp(-2 * k ** 2 * c ** 2)
                 for k in range(1, 200))
print(f"\nc = sqrt({n1 * n2}/{N}) * {D_obs} = {c:.6f}")
print(f"급수 1 - K(c)        = {series:.6f}")
print(f"kstwobign.sf(c)      = {stats.kstwobign.sf(c):.6f}")
print(f"n_eff = {n1 * n2 / N} -> round -> {np.round(n1 * n2 / N):.0f}")
print(f"kstwo.sf(D, round(n_eff)) = {stats.kstwo.sf(D_obs, int(np.round(n1 * n2 / N))):.6f}")
print(f"scipy method='asymp'      = {stats.ks_2samp(A, B, method='asymp').pvalue:.6f}")

정렬한 합친 표본 [10.0, 10.1, 10.2, 10.2, 10.3, 10.4, 10.4, 10.5, 10.6, 10.8]
동점 [(10.2, 2), (10.4, 2)]

열거한 배치 수 252  (C(10,5) = 252)
D 가 가질 수 있는 값 [0.2, 0.4, 0.6, 0.8, 1.0]
배치 수           [32, 130, 70, 18, 2]

D >= 0.4 인 배치 220 개  ->  정확 p = 220/252 = 0.873015873
scipy method='exact' : 0.873015873
이 검정이 낼 수 있는 p-값 전부:
  D = 0.2:  p = 252/252 = 1.0000
  D = 0.4:  p = 220/252 = 0.8730
  D = 0.6:  p = 90/252 = 0.3571
  D = 0.8:  p = 20/252 = 0.0794
  D = 1.0:  p = 2/252 = 0.0079

c = sqrt(25/10) * 0.4 = 0.632456
급수 1 - K(c)        = 0.818621
kstwobign.sf(c)      = 0.818621
n_eff = 2.5 -> round -> 2
kstwo.sf(D, round(n_eff)) = 0.820000
scipy method='asymp'      = 0.820000

from scipy import stats
c1 = [72, 78, 85, 90, 65]; c2 = [80, 88, 92, 95, 85, 76]
print(stats.ks_2samp(c1, c2, method='exact'))
# statistic=0.4333, pvalue=0.5909
print(stats.mannwhitneyu(c1, c2, method='exact').pvalue)   # 0.2468

KstestResult(statistic=0.43333333333333335, pvalue=0.5909090909090909, statistic_location=78, statistic_sign=1)
0.24675324675324672

import numpy as np
from scipy import stats
rng = np.random.default_rng(0)
B, n = 3000, 25

def rates(gen):
    pk, pm = [], []
    for _ in range(B):
        x, y = gen()
        pk.append(stats.ks_2samp(x, y).pvalue)
        pm.append(stats.mannwhitneyu(x, y).pvalue)
    return round(np.mean(np.array(pk) < .05), 3), \
           round(np.mean(np.array(pm) < .05), 3)

for d in (0.0, 0.8):
    print("shift", d, rates(lambda: (rng.normal(0, 1, n), rng.normal(d, 1, n))))
for s in (1.0, 2.0, 3.0):
    print("scale", s, rates(lambda: (rng.normal(0, 1, n), rng.normal(0, s, n))))

shift 0.0 (0.039, 0.047)
shift 0.8 (0.614, 0.757)
scale 1.0 (0.037, 0.047)
scale 2.0 (0.153, 0.058)
scale 3.0 (0.414, 0.069)

import numpy as np
from scipy import stats
rng = np.random.default_rng(1)
for name, gen in [("normal", lambda: rng.normal(0, 1, 30)),
                  ("exp",    lambda: rng.exponential(1, 30)),
                  ("cauchy", lambda: rng.standard_cauchy(30)),
                  ("unif",   lambda: rng.random(30))]:
    D = [stats.ks_2samp(gen(), gen()).statistic for _ in range(20000)]
    print(name, round(np.mean(D), 4), round(np.std(D), 4))

normal 0.2086 0.0667
exp 0.2075 0.0663
cauchy 0.2086 0.0676
unif 0.2079 0.0667

import numpy as np
from scipy import stats
rng = np.random.default_rng(7)
n = 500

cases = {
    "위치이동":   (rng.normal(0, 1, n),  rng.normal(1, 1, n)),
    "척도차이":   (rng.normal(0, 1, n),  rng.normal(0, 2.5, n)),
    "오른쪽 꼬리": (rng.normal(0, 1, n),  rng.standard_t(3, n)),
}
for name, (x, y) in cases.items():
    r = stats.ks_2samp(x, y)
    print(name, round(r.statistic, 3), round(r.statistic_location, 3))

위치이동 0.478 0.217
척도차이 0.22 -1.57
오른쪽 꼬리 0.078 1.199