Kolmogorov-Smirnov 검정¶

(1) 손계산이 끝자리까지 맞고, 점근식은 \(1.2\%\) 어긋난다. 경험분포함수는 각 관측값에서 뛰는 계단이므로 두 방향의 틈을 따로 잰다. 모수를 못박았으니 표준화 없이 \(F_0 = \Phi\)를 바로 쓴다.

이고 \(D = \max = 0.0354148220\)으로 scipy의 값과 열 자리까지 같다.

\(p\)값은 콜모고로프 분포에서 온다. \(t = \sqrt{n}\,D = \sqrt{250}\times0.0354148 = 0.559958\)에 대해 점근 공식은

이고, 이 급수는 매우 빨리 수렴하므로 100항만 더해도 충분하다. 값은 \(0.912469\)다. 그런데 scipy가 돌려주는 것은

로 \(1.2\%\) 작다. 까닭은 scipy가 점근식 대신 \(n\)에 대한 정확 분포(scipy.stats.kstwo)를 쓰기 때문이다. \(n = 250\)은 점근식이 쓸 만한 크기이지만 아직 완전히 같지는 않다. 급수식으로 직접 계산한 값을 그대로 보고하면 셋째 자리가 틀린다.

(2) 정확히 지킨다. 진짜 표준정규 자료에 같은 절차를 4000번 적용하면

이다. \(p\)값이 \([0,1]\)에 균등하다. 귀무가설이 참일 때 \(p\)값이 균등해야 한다는 교과서의 약속이 그대로 지켜진다.

이 표를 KS 와 릴리에포르 보기 1의 표와 나란히 놓는 것이 이 보기의 요점이다. 같은 stats.kstest 함수인데

로 백 배 갈린다. 함수가 틀린 것이 아니다. 함수가 전제하는 바를 지켰는가 아닌가의 차이이고, 그 전제는 "\(F_0\)가 자료를 보기 전에 완전히 정해져 있다"는 것이다.

실무에서 이 전제가 성립하는 경우는 생각보다 드물다. 정규성을 검정하겠다고 할 때 \(\mu\)와 \(\sigma\)를 미리 알고 있는 일은 거의 없다. 그러므로 kstest로 정규성을 검정하는 코드를 보면 거의 언제나 틀렸다고 보아야 하고, 맞는 경우는 균등성 검정처럼 귀무분포가 모수 없이 정해지는 때다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
x = rng.normal(0.0, 1.0, size=250)

# 귀무가설이 N(0,1) 로 완전히 못박혀 있는 경우다. args 로 모수를 직접
# 넘겨 주었고 자료에서 추정하지 않았으므로, 이 p-값은 그대로 믿을 수 있다.
# 모수를 자료에서 뽑아 쓰면 Lilliefors 로 가야 한다.
# D 는 경험분포함수와 이론분포함수의 최대 세로 거리다.
D, p = stats.kstest(x, 'norm', args=(0.0, 1.0))

print(f"n = {x.size}")
print(f"KS one-sample vs N(0,1): D = {D:.4f}, p = {p:.4f}")
if p < 0.05:
    print("=> Reject H0: data may not follow N(0,1).")
else:
    print("=> Fail to reject H0 at alpha = 0.05.")

출력:

n = 250
KS one-sample vs N(0,1): D = 0.0354, p = 0.9013
=> Fail to reject H0 at alpha = 0.05.

손계산과 실제 크기를 함께 확인한다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
x = rng.normal(0.0, 1.0, size=250)
n = x.size

# D 를 정의대로 손계산한다.
z = np.sort(x)
F = stats.norm.cdf(z)
i = np.arange(1, n + 1)
d_plus = np.max(i / n - F)
d_minus = np.max(F - (i - 1) / n)
print(f"D+ = {d_plus:.10f},  D- = {d_minus:.10f},  D = {max(d_plus, d_minus):.10f}")
print(f"scipy D = {stats.kstest(x, 'norm', args=(0.0, 1.0))[0]:.10f}")

# p 값: 모수를 못박았으므로 콜모고로프 분포를 쓸 수 있다.
D = max(d_plus, d_minus)
t = np.sqrt(n) * D
p_asym = 2 * sum((-1)**(k - 1) * np.exp(-2 * k**2 * t**2) for k in range(1, 101))
print(f"\nsqrt(n) D = {t:.6f}")
print(f"점근 p = 2 sum (-1)^(k-1) exp(-2k^2 t^2) = {p_asym:.6f}")
print(f"scipy 의 점근 분포 kstwobign.sf(t)        = {stats.kstwobign.sf(t):.6f}")
print(f"정확 분포 kstwo.sf(D, n)                 = {stats.kstwo.sf(D, n):.6f}")
print(f"kstest p                                 = {stats.kstest(x, 'norm', args=(0.0, 1.0))[1]:.6f}")

# 모수를 못박은 KS 는 명목 크기를 지키는가
rng2 = np.random.default_rng(11)
R = 4000
ps = np.array([stats.kstest(rng2.standard_normal(n), 'norm', args=(0.0, 1.0))[1]
               for _ in range(R)])
print(f"\n모수를 못박았을 때의 실제 기각률 (R = {R})")
for a in (0.05, 0.10):
    print(f"  명목 {a:.2f}: {np.mean(ps < a):.4f}  (MC SE {np.sqrt(a * (1 - a) / R):.4f})")
print(f"  p 값의 평균 = {ps.mean():.4f}  (균등이면 0.5)")

출력:

D+ = 0.0354148220,  D- = 0.0302637669,  D = 0.0354148220
scipy D = 0.0354148220

sqrt(n) D = 0.559958
점근 p = 2 sum (-1)^(k-1) exp(-2k^2 t^2) = 0.912469
scipy 의 점근 분포 kstwobign.sf(t)        = 0.912469
정확 분포 kstwo.sf(D, n)                 = 0.901284
kstest p                                 = 0.901284

모수를 못박았을 때의 실제 기각률 (R = 4000)
  명목 0.05: 0.0500  (MC SE 0.0034)
  명목 0.10: 0.1010  (MC SE 0.0047)
  p 값의 평균 = 0.5014  (균등이면 0.5)

세 가지가 맞아떨어진다. 손계산한 \(D\)가 scipy와 열 자리까지 같고, 손으로 더한 급수가 scipy의 점근 분포 함수와 여섯 자리까지 같으며(둘 다 \(0.912469\)), 그 점근값이 정확값 \(0.901284\)에서 \(1.2\%\) 떨어져 있다. 그리고 모수를 못박은 검정은 명목 크기를 몬테카를로 오차 안에서 지킨다. \(\square\)

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
x = rng.normal(0, 1, size=300)

D1, p1 = stats.kstest(x, 'norm', args=(0.0, 1.0))
D2, p2 = stats.kstest(x, 'norm', args=(0.5, 1.0))

print(f"vs N(0,1):   D = {D1:.4f}, p = {p1:.4g}")
print(f"vs N(0.5,1): D = {D2:.4f}, p = {p2:.4g}")

출력:

vs N(0,1):   D = 0.0507, p = 0.4089
vs N(0.5,1): D = 0.2264, p = 5.45e-14

참 분포인 \(\mathcal{N}(0,1)\)에 대해서는 \(D = 0.051\), \(p = 0.409\)로 기각하지 않는다. 평균이 틀린 \(\mathcal{N}(0.5,1)\)에 대해서는 EDF가 체계적으로 이동해 있어 \(D = 0.226\)으로 네 배 넘게 커지고 \(p = 5.5 \times 10^{-14}\)로 강하게 기각한다.

\(D\)의 크기를 직관적으로 이해해 보자. \(\Phi(x)\)와 \(\Phi(x - 0.5)\) 사이의 최대 수직 거리는 \(x = 0.25\)에서 \(\Phi(0.25) - \Phi(-0.25) = 0.1974\)이다. 표본 \(D = 0.226\)이 이 이론값 근처에 있다.

위치 이동은 KS 검정이 잘 잡는다. 다음 연습문제들에서 볼 꼬리 이탈과 대조된다. \(\square\)

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
x = rng.normal(0, 1, size=200)

mu_hat, sigma_hat = x.mean(), x.std(ddof=1)
D, p = stats.kstest(x, 'norm', args=(mu_hat, sigma_hat))

print(f"Estimated mu = {mu_hat:.4f}, sigma = {sigma_hat:.4f}")
print(f"KS vs N(mu_hat, sigma_hat): D = {D:.4f}, p = {p:.4f}")

출력:

Estimated mu = 0.0153, sigma = 0.9636
KS vs N(mu_hat, sigma_hat): D = 0.0306, p = 0.9892

\(p = 0.989\)는 믿을 수 없다. 적합된 정규분포가 미리 지정된 정규분포보다 EDF에 더 가까우므로 \(D\)가 작아지고 \(p\)값이 부풀려진다.

얼마나 심각한지 모의실험으로 확인해 보자. \(\mathcal{N}(0,1)\)에서 \(n = 200\)인 표본 5,000개를 뽑아 매번 모수를 추정하고 이 순진한 KS 검정을 적용하면 \(\alpha = 0.05\)에서의 경험적 기각률은

Plug-in KS empirical size: 0.0002

0.0002이다. 명목값 0.05의 250분의 1에 불과하다. 5,000회 중 단 한 번만 기각했다는 뜻이다.

이 정도로 극단적인 과소기각은 검정을 사실상 무용하게 만든다. 크기가 0.0002인 검정은 참된 이탈이 있어도 거의 탐지하지 못한다. 타당한 추론을 하려면 반드시 Lilliefors 보정이 필요하다. \(\square\)

EDF \(F_n(x)\)는 각 \(X_{(i)}\)에서 \(1/n\)만큼 뛰는 계단함수이다. 연속된 두 순서통계량 사이에서 \(F_n\)은 상수인 반면 \(F_0\)은 단조증가한다. 따라서 차이 \(F_n(x) - F_0(x)\)는 그 구간에서 단조감소하고, 상한 \(|F_n(x) - F_0(x)|\)는 도약 직전이나 직후에 달성된다.

• \(i\)번째 도약 직후: \(F_n(X_{(i)}) = i/n\)이므로 \(|i/n - F_0(X_{(i)})|\).
• \(i\)번째 도약 직전: \(F_n(X_{(i)}^-) = (i-1)/n\)이므로 \(|(i-1)/n - F_0(X_{(i)})|\).

따라서

\(n\)개의 순서통계량만 평가하면 충분하다. 무한한 실직선 위의 상한을 유한한 최댓값 계산으로 바꿔주는 결과이며, KS 검정을 실용적으로 만드는 핵심이다. \(\square\)

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n, reps, alpha = 200, 5000, 0.05
rejections = 0

for _ in range(reps):
    x = rng.standard_t(df=5, size=n)
    _, p = stats.kstest(x, 'norm', args=(0.0, 1.0))
    if p < alpha:
        rejections += 1

print(f"KS power vs t(5): {rejections / reps:.4f}")

출력:

KS power vs t(5): 0.0910

검정력이 0.091에 불과하다. 크기 0.05보다 겨우 조금 클 뿐이며, 사실상 쓸모가 없다.

왜 이렇게 낮은지는 모집단 KS 거리를 계산하면 명확해진다.

한편 \(n = 200\)에서 5% 임계값은 근사적으로

곧 모집단 거리 \(0.0305\)가 임계값 \(0.096\)의 3분의 1에도 못 미친다. 표집잡음이 없더라도 \(D_n\)이 임계값에 도달하지 못한다. 50% 검정력을 얻으려면 \(n \approx (1.358/0.0305)^2 \approx 1980\)이 필요하다.

\(t_5\)의 분산이 \(5/3 \approx 1.67\)로 1과 크게 다른데도 CDF 거리가 작은 이유는, 두 CDF가 여러 번 교차하면서 차이가 상쇄되기 때문이다. 중앙 근처에서는 \(t_5\)가 더 뾰족해서 \(F_{t_5} > \Phi\)이고, 꼬리에서는 \(t_5\)가 더 두꺼워서 부호가 뒤집힌다. 어느 한 지점의 수직 거리는 작게 유지된다.

같은 자료에 정규성 검정을 적용하면 대비가 뚜렷하다(\(n = 200\), \(\alpha = 0.05\)).

(공정한 비교는 아니다. AD와 SW는 모수를 추정하는 위치-척도 자유 검정인 반면 KS는 완전 지정 대립을 다룬다. 그럼에도 KS가 꼬리 이탈에 얼마나 둔감한지는 분명하다.)

꼬리에 가중을 주는 Anderson-Darling과 순서통계량의 상관구조를 활용하는 Shapiro-Wilk가 이런 상황에서 압도적으로 낫다. \(\square\)

정리 (Glivenko-Cantelli). \(X_1, X_2, \ldots\)가 CDF \(F\)를 갖는 i.i.d.라 하자. 그러면

증명 개요. 임의의 \(x \in \mathbb{R}\)을 고정한다. 강대수의 법칙에 의해 \(F_n(x) = \frac{1}{n}\sum_{i=1}^n \mathbf{1}(X_i \leq x) \to F(x)\)가 거의 확실하게 성립한다. 이는 점별 수렴이다.

균등수렴으로 끌어올리려면, \(\epsilon > 0\)을 고정하고 모든 \(j\)에 대해 \(F(t_j) - F(t_{j-1}) < \epsilon\)이 되도록 점 \(-\infty = t_0 < t_1 < \cdots < t_K = \infty\)를 잡는다. 임의의 \(x \in [t_{j-1}, t_j]\)에 대해 \(F_n\)과 \(F\)의 단조성에서

비슷하게 \(F(x) - F_n(x) < [F(t_j) - F_n(t_{j-1})] + \epsilon\)이다.

\(t_j\)가 유한개이므로 각 \(t_j\)에서의 점별 수렴(강대수의 법칙)이 결국 \(\sup_x |F_n(x) - F(x)| < 2\epsilon\)을 거의 확실하게 보장한다. \(\epsilon\)이 임의였으므로 결과가 따라 나온다. \(\square\)