Kolmogorov-Smirnov 검정과 Lilliefors 검정¶

import numpy as np
from scipy import stats

np.random.seed(0)

# 예시 자료를 만든다
# data = np.random.normal(0, 1, 1000)
data = np.random.normal(1, 10, 1000)

# kstest 는 모수를 미리 알고 있다고 전제한다. 그래서 표준화해 N(0,1) 과
# 견주는데, 그 평균과 표준편차를 자료에서 뽑아 썼다는 것이 함정이다.
# 이러면 검정이 지나치게 너그러워져 p-값이 실제보다 크게 나온다.
# 아래 Lilliefors 가 바로 이 문제를 고친 검정이다.
data_ks = (data - data.mean()) / data.std()

stat, p_value = stats.kstest(data_ks, 'norm')
print(f"Kolmogorov-Smirnov Test: Statistic={stat}, p-value={p_value}")

alpha = 0.05
if p_value <= alpha:
    print("Reject H_0: The data is not normally distributed.")
else:
    print("Fail to reject H_0: The data is normally distributed.")

Kolmogorov-Smirnov Test: Statistic=0.01903411267034605, p-value=0.8547733408587939
Fail to reject H_0: The data is normally distributed.

import numpy as np
from scipy import stats

np.random.seed(0)
data = np.random.normal(1, 10, 1000)
n = data.size

# 정의대로 D 를 손계산한다. 계단의 양쪽을 모두 본다.
def D_of(x, ddof):
    z = np.sort((x - x.mean()) / x.std(ddof=ddof))
    F = stats.norm.cdf(z)
    i = np.arange(1, len(z) + 1)
    d_plus = np.max(i / len(z) - F)
    d_minus = np.max(F - (i - 1) / len(z))
    return max(d_plus, d_minus), d_plus, d_minus

D0, dp0, dm0 = D_of(data, 0)
D1, dp1, dm1 = D_of(data, 1)
print(f"ddof=0:  D+ = {dp0:.10f},  D- = {dm0:.10f},  D = {D0:.10f}")
print(f"ddof=1:  D+ = {dp1:.10f},  D- = {dm1:.10f},  D = {D1:.10f}")

data_ks = (data - data.mean()) / data.std()
print(f"\nkstest   D = {stats.kstest(data_ks, 'norm')[0]:.10f}")

# 이 "검정" 의 실제 크기를 재 본다.
rng = np.random.default_rng(7)
R = 4000
ps = np.array([stats.kstest((lambda x: (x - x.mean()) / x.std())(
    rng.standard_normal(n)), 'norm')[1] for _ in range(R)])
for a in (0.05, 0.10, 0.20):
    print(f"명목 {a:.2f}: 실제 기각률 = {np.mean(ps < a):.5f}"
          f"  ({np.sum(ps < a)}/{R})")
print(f"\np 값의 평균 = {ps.mean():.4f}  (균등이면 0.5),  최솟값 = {ps.min():.4f}")

ddof=0:  D+ = 0.0190341127,  D- = 0.0098858448,  D = 0.0190341127
ddof=1:  D+ = 0.0191252945,  D- = 0.0098072163,  D = 0.0191252945

kstest   D = 0.0190341127
명목 0.05: 실제 기각률 = 0.00050  (2/4000)
명목 0.10: 실제 기각률 = 0.00150  (6/4000)
명목 0.20: 실제 기각률 = 0.00725  (29/4000)

p 값의 평균 = 0.7869  (균등이면 0.5),  최솟값 = 0.0408

import numpy as np
from scipy import stats
from statsmodels.stats.diagnostic import lilliefors

np.random.seed(0)

# 예시 자료를 만든다
# data = np.random.normal(0, 1, 1000)
data = np.random.normal(1, 10, 1000)

data_ks = (data - data.mean()) / data.std()

# 같은 자료에 두 검정을 나란히 돌린다. 통계량은 거의 같은데 p-값이 갈린다.
stat, p_value = stats.kstest(data_ks, 'norm')
print(f"Kolmogorov-Smirnov Test: Statistic={stat}, p-value={p_value}")

# Lilliefors 는 모수를 자료에서 추정했다는 사실을 셈에 넣은 기각값을 쓴다.
# 그래서 같은 통계량에도 더 작은 p-값을 준다. 모수를 추정했다면 이쪽이 맞다.
stat, p_value = lilliefors(data)
print(f"Lilliefors Test: Statistic={stat}, p-value={p_value}")

# 결과 해석
alpha = 0.05
if p_value <= alpha:
    print("Reject H_0: The data is not normally distributed.")
else:
    print("Fail to reject H_0: The data is normally distributed.")

Kolmogorov-Smirnov Test: Statistic=0.01903411267034605, p-value=0.8547733408587939
Lilliefors Test: Statistic=0.019125294462402076, p-value=0.5818164701330186
Fail to reject H_0: The data is normally distributed.

import numpy as np
from scipy import stats
from statsmodels.stats.diagnostic import lilliefors

np.random.seed(0)
data = np.random.normal(1, 10, 1000)

z0 = (data - data.mean()) / data.std(ddof=0)
z1 = (data - data.mean()) / data.std(ddof=1)
print(f"ddof=0 으로 표준화한 KS 통계량 = {stats.kstest(z0, 'norm')[0]:.12f}")
print(f"ddof=1 으로 표준화한 KS 통계량 = {stats.kstest(z1, 'norm')[0]:.12f}")
print(f"lilliefors 의 통계량           = {lilliefors(data)[0]:.12f}")
print(f"\n같은 통계량에 붙는 두 p 값")
print(f"  kstest     p = {stats.kstest(z0, 'norm')[1]:.6f}")
print(f"  lilliefors p = {lilliefors(data)[1]:.6f}")
print(f"  ddof=1 통계량을 kstest 의 귀무분포로 읽으면 p = {stats.kstest(z1, 'norm')[1]:.6f}")

# 릴리에포르가 실제 크기를 얼마나 바로잡는가
rng = np.random.default_rng(1234)
R = 2000
print(f"\n정규 자료에서의 실제 기각률 (R = {R})")
for n in (50, 200, 1000):
    pl = np.array([lilliefors(rng.standard_normal(n))[1] for _ in range(R)])
    pk = np.array([stats.kstest((lambda x: (x - x.mean()) / x.std())(
        rng.standard_normal(n)), 'norm')[1] for _ in range(R)])
    print(f"{n:>6}  릴리에포르 {np.mean(pl < 0.05):.4f}  KS {np.mean(pk < 0.05):.4f}"
          f"   |  릴리에포르 {np.mean(pl < 0.10):.4f}  KS {np.mean(pk < 0.10):.4f}")
print(f"\n몬테카를로 표준오차: 0.05 에서 {np.sqrt(0.05 * 0.95 / R):.4f},"
      f"  0.10 에서 {np.sqrt(0.10 * 0.90 / R):.4f}")

ddof=0 으로 표준화한 KS 통계량 = 0.019034112670
ddof=1 으로 표준화한 KS 통계량 = 0.019125294462
lilliefors 의 통계량           = 0.019125294462

같은 통계량에 붙는 두 p 값
  kstest     p = 0.854773
  lilliefors p = 0.581816
  ddof=1 통계량을 kstest 의 귀무분포로 읽으면 p = 0.850822

정규 자료에서의 실제 기각률 (R = 2000)
    50  릴리에포르 0.0540  KS 0.0000   |  릴리에포르 0.1080  KS 0.0005
   200  릴리에포르 0.0470  KS 0.0000   |  릴리에포르 0.0970  KS 0.0010
  1000  릴리에포르 0.0365  KS 0.0000   |  릴리에포르 0.0745  KS 0.0020

몬테카를로 표준오차: 0.05 에서 0.0049,  0.10 에서 0.0067