형식적 정규성 검정 모음¶

import numpy as np
from scipy import stats

# 참으로 정규인 자료에 다섯 검정을 모두 돌려 견준다. 모두 기각하지
# 않아야 정상이다. 검정마다 무엇에 민감한지가 달라, 정규에서 벗어난
# 자료로 바꿔 돌려 보면 차이가 드러난다.
rng = np.random.default_rng(0)
data = rng.normal(0, 1, size=100)

# Shapiro-Wilk — 작은 표본에서 검정력이 가장 좋다. 두루 쓰기 좋은 기본값이다.
W, p_sw = stats.shapiro(data)
print(f"Shapiro-Wilk:     W = {W:.4f}, p = {p_sw:.4g}")

# D'Agostino K^2 — 왜도와 첨도만 본다. 왜 기각되었는지 알기 쉽다.
K2, p_k2 = stats.normaltest(data)
print(f"D'Agostino K^2:   K2 = {K2:.4f}, p = {p_k2:.4g}")

# Jarque-Bera — K^2 과 비슷하지만 표본이 아주 클 때라야 근사가 맞는다.
JB, p_jb = stats.jarque_bera(data)
print(f"Jarque-Bera:      JB = {JB:.4f}, p = {p_jb:.4g}")

# KS — 모수를 못박은 경우에만 이 p-값이 맞다. 자료에서 추정했다면
# Lilliefors 로 가야 한다.
D, p_ks = stats.kstest(data, 'norm', args=(0, 1))
print(f"KS (vs N(0,1)):   D = {D:.4f}, p = {p_ks:.4g}")

# Anderson-Darling — 꼬리 쪽 이탈에 민감하다. p-값 대신 기각값 표를 준다.
ad = stats.anderson(data, dist="norm")
print(f"Anderson-Darling: A^2 = {ad.statistic:.4f}")
for cv, sl in zip(ad.critical_values, ad.significance_level):
    print(f"  {sl:.0f}% critical value: {cv:.4f}")

Shapiro-Wilk:     W = 0.9819, p = 0.1873
D'Agostino K^2:   K2 = 1.9679, p = 0.3738
Jarque-Bera:      JB = 1.4236, p = 0.4908
KS (vs N(0,1)):   D = 0.0723, p = 0.6465
Anderson-Darling: A^2 = 0.4494
  15% critical value: 0.5550
  10% critical value: 0.6320
  5% critical value: 0.7590
  2% critical value: 0.8850
  1% critical value: 1.0530

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
data = rng.normal(0, 1, size=100)
n = data.size

# (1) 다섯 검정 가운데 KS 만 귀무가설이 다르다. 자료를 0.5 만큼 옮겨 본다.
shift = data + 0.5
print("원자료와 0.5 옮긴 자료의 p 값")
print(f"{'검정':<18}{'원자료':>12}{'+0.5':>12}")
print(f"{'Shapiro-Wilk':<18}{stats.shapiro(data)[1]:>12.4g}{stats.shapiro(shift)[1]:>12.4g}")
print(f"{'DAgostino K^2':<18}{stats.normaltest(data)[1]:>12.4g}{stats.normaltest(shift)[1]:>12.4g}")
print(f"{'Jarque-Bera':<18}{stats.jarque_bera(data)[1]:>12.4g}{stats.jarque_bera(shift)[1]:>12.4g}")
print(f"{'KS vs N(0,1)':<18}{stats.kstest(data, 'norm', args=(0, 1))[1]:>12.4g}"
      f"{stats.kstest(shift, 'norm', args=(0, 1))[1]:>12.4g}")
print(f"{'Anderson A^2':<18}{stats.anderson(data).statistic:>12.4f}"
      f"{stats.anderson(shift).statistic:>12.4f}")

# (2) 같은 n 에서 두 대립에 대한 검정력
def power(draw, R, alpha=0.05):
    out = {k: 0 for k in ("SW", "K2", "JB", "KS", "AD")}
    for _ in range(R):
        x = draw()
        out["SW"] += stats.shapiro(x)[1] < alpha
        out["K2"] += stats.normaltest(x)[1] < alpha
        out["JB"] += stats.jarque_bera(x)[1] < alpha
        out["KS"] += stats.kstest(x, 'norm', args=(0, 1))[1] < alpha
        r = stats.anderson(x)
        out["AD"] += r.statistic > r.critical_values[2]
    return {k: v / R for k, v in out.items()}

rg = np.random.default_rng(3)
R = 4000
print(f"\nn = {n}, alpha = 0.05, R = {R}  (MC SE 는 0.05 에서 0.0034)")
rows = [("정규 N(0,1) — 크기", lambda: rg.standard_normal(n)),
        ("t_5 — 두꺼운 꼬리", lambda: rg.standard_t(5, n)),
        ("Lognormal(0,0.5) 표준화 — 치우침",
         lambda: (lambda z: (z - z.mean()) / z.std(ddof=1))(rg.lognormal(0, 0.5, n)))]
print(f"{'대립':<34}{'SW':>8}{'K2':>8}{'JB':>8}{'KS':>8}{'AD':>8}")
for name, d in rows:
    pw = power(d, R)
    print(f"{name:<34}{pw['SW']:>8.4f}{pw['K2']:>8.4f}{pw['JB']:>8.4f}"
          f"{pw['KS']:>8.4f}{pw['AD']:>8.4f}")

원자료와 0.5 옮긴 자료의 p 값
검정                         원자료        +0.5
Shapiro-Wilk            0.1873      0.1873
DAgostino K^2           0.3738      0.3738
Jarque-Bera             0.4908      0.4908
KS vs N(0,1)            0.6465   1.562e-05
Anderson A^2            0.4494      0.4494

n = 100, alpha = 0.05, R = 4000  (MC SE 는 0.05 에서 0.0034)
대립                                      SW      K2      JB      KS      AD
정규 N(0,1) — 크기                      0.0558  0.0580  0.0440  0.0475  0.0475
t_5 — 두꺼운 꼬리                        0.5503  0.5880  0.6122  0.0630  0.4710
Lognormal(0,0.5) 표준화 — 치우침          0.9990  0.9925  0.9905  0.4095  0.9932

import numpy as np
from scipy import stats

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

W, p_sw = stats.shapiro(data)
K2, p_k2 = stats.normaltest(data)
JB, p_jb = stats.jarque_bera(data)

print(f"Shapiro-Wilk:   p = {p_sw:.4g}")
print(f"D'Agostino K^2: p = {p_k2:.4g}")
print(f"Jarque-Bera:    p = {p_jb:.4g}")

Shapiro-Wilk:   p = 0.4893
D'Agostino K^2: p = 0.301
Jarque-Bera:    p = 0.3226

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
data = rng.lognormal(0, 0.5, size=200)

W, p_sw = stats.shapiro(data)
K2, p_k2 = stats.normaltest(data)
JB, p_jb = stats.jarque_bera(data)

print(f"Shapiro-Wilk:   p = {p_sw:.4g}")
print(f"D'Agostino K^2: p = {p_k2:.4g}")
print(f"Jarque-Bera:    p = {p_jb:.4g}")

Shapiro-Wilk:   p = 1.153e-12
D'Agostino K^2: p = 1.664e-22
Jarque-Bera:    p = 3.108e-107

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n, reps, alpha = 50, 10000, 0.05
rejections = 0

for _ in range(reps):
    x = rng.normal(0, 1, size=n)
    _, p = stats.shapiro(x)
    if p < alpha:
        rejections += 1

empirical_size = rejections / reps
print(f"Empirical size: {empirical_size:.4f}")

Empirical size: 0.0487