올바른 검정 고르기¶

import numpy as np
from scipy import stats

# ===================================================================
# 같은 자료에 네 검정을 모두 돌려 결론을 견준다
#
# 검정마다 무엇에 민감한지가 달라, 같은 자료에서도 p-값이 꽤 갈린다.
# 하나만 골라 돌리고 끝낼 일이 아니라는 것이 요점이다.
# ===================================================================

np.random.seed(42)
n = 100

# 자유도 5 인 t. 대칭이지만 꼬리가 두꺼운 자료다.
data = stats.t.rvs(df=5, size=n)

if __name__ == "__main__":
    # Shapiro-Wilk — 두루 쓰기 좋고 작은 표본에서 검정력이 높다.
    sw_stat, sw_p = stats.shapiro(data)
    print(f"Shapiro-Wilk:      W = {sw_stat:.4f}, p = {sw_p:.4f}")

    # Anderson-Darling — 꼬리에 무게를 싣는다. 이 자료에 가장 민감할 것이다.
    ad_result = stats.anderson(data, dist="norm")
    print(f"Anderson-Darling:  A2 = {ad_result.statistic:.4f}, "
          f"critical (5%) = {ad_result.critical_values[2]:.4f}")

    # KS — 모수를 자료에서 추정해 넘겼으므로 p-값이 실제보다 크게 나온다.
    # 제대로 하려면 Lilliefors 기각값을 써야 한다.
    ks_stat, ks_p = stats.kstest(data, "norm", args=(np.mean(data), np.std(data)))
    print(f"KS (estimated):    D = {ks_stat:.4f}, p = {ks_p:.4f}")

    # Jarque-Bera — 왜도와 첨도만 본다. 큰 표본에서 쓸 만하다.
    jb_stat, jb_p = stats.jarque_bera(data)
    print(f"Jarque-Bera:       JB = {jb_stat:.4f}, p = {jb_p:.4f}")

Shapiro-Wilk:      W = 0.6287, p = 0.0000
Anderson-Darling:  A2 = 6.0914, critical (5%) = 0.7590
KS (estimated):    D = 0.1783, p = 0.0030
Jarque-Bera:       JB = 5365.9041, p = 0.0000

import numpy as np
from scipy import stats

np.random.seed(42)
n = 100
d = stats.t.rvs(df=5, size=n)

print("t_5 의 이론값:  왜도 0,  초과첨도 6/(5-4) = 6")
print(f"이 표본:  g1 = {stats.skew(d):.4f}   g2 = {stats.kurtosis(d):.4f}")
print(f"  (보정판 G1 = {stats.skew(d, bias=False):.4f},"
      f" G2 = {stats.kurtosis(d, bias=False):.4f})")

# 자르크-베라의 닫힌 꼴을 직접 확인한다.
g1, g2 = stats.skew(d), stats.kurtosis(d)
print(f"\nJB = n/6*(g1^2 + g2^2/4) = {n / 6 * (g1 ** 2 + g2 ** 2 / 4):.4f}"
      f"   scipy = {stats.jarque_bera(d).statistic:.4f}")
print(f"t_5 의 이론값을 넣으면 JB = {n / 6 * (0 ** 2 + 6 ** 2 / 4):.1f}")

# 이 표본의 극단값 하나가 모든 것을 끌고 간다.
z = (d - d.mean()) / d.std(ddof=1)
print(f"\n표준화 최대/최소:  {z.max():+.2f}  {z.min():+.2f}")
drop = np.sort(d)[:-1]
print("가장 큰 관측 하나를 빼면 (n = 99):")
print(f"  g1 = {stats.skew(drop):+.4f}   g2 = {stats.kurtosis(drop):+.4f}")
print(f"  W = {stats.shapiro(drop).statistic:.4f}  p = {stats.shapiro(drop).pvalue:.4f}")
print(f"  A2 = {stats.anderson(drop, dist='norm').statistic:.4f}")
print(f"  JB = {stats.jarque_bera(drop).statistic:.4f}"
      f"  p = {stats.jarque_bera(drop).pvalue:.4f}")

# KS 를 추정 모수로 쓰면 p 값이 얼마나 과대평가되는가.
D_obs = stats.kstest(d, "norm", args=(np.mean(d), np.std(d))).statistic


def ks_est(x):
    return stats.kstest(x, "norm", args=(np.mean(x), np.std(x))).statistic


rng = np.random.default_rng(404)
R = 40000
Dnull = np.array([ks_est(rng.normal(0, 1, n)) for _ in range(R)])
p_lillie = np.mean(Dnull >= D_obs)
print(f"\nKS:  D = {D_obs:.4f}")
print(f"  표준 KS p값 (모수를 안다고 가정)        = "
      f"{stats.kstest(d, 'norm', args=(np.mean(d), np.std(d))).pvalue:.4f}")
print(f"  릴리에포르 p값 (모수 추정을 반영, R={R}) = {p_lillie:.5f}")
print(f"  5% 임계값:  표준 KS {stats.ksone.ppf(0.975, n):.4f}"
      f"   릴리에포르 {np.percentile(Dnull, 95):.4f}")

# 그 임계값을 잘못 쓰면 검정의 실제 크기는 얼마가 되는가.
crit_std = stats.ksone.ppf(0.975, n)
print(f"\n표준 KS 임계값 {crit_std:.4f} 를 추정 모수에 그대로 쓰면")
print(f"  명목 0.05 인 검정의 실제 크기 = {np.mean(Dnull > crit_std):.5f}")

t_5 의 이론값:  왜도 0,  초과첨도 6/(5-4) = 6
이 표본:  g1 = 4.8256   g2 = 34.5640
  (보정판 G1 = 4.8994, G2 = 36.4190)

JB = n/6*(g1^2 + g2^2/4) = 5365.9041   scipy = 5365.9041
t_5 의 이론값을 넣으면 JB = 150.0

표준화 최대/최소:  +7.75  -1.45
가장 큰 관측 하나를 빼면 (n = 99):
  g1 = +0.8195   g2 = +2.7637
  W = 0.9548  p = 0.0019
  A2 = 0.6086
  JB = 42.5870  p = 0.0000

KS:  D = 0.1783
  표준 KS p값 (모수를 안다고 가정)        = 0.0030
  릴리에포르 p값 (모수 추정을 반영, R=40000) = 0.00000
  5% 임계값:  표준 KS 0.1340   릴리에포르 0.0892

표준 KS 임계값 0.1340 를 추정 모수에 그대로 쓰면
  명목 0.05 인 검정의 실제 크기 = 0.00015