민감도와 실질적 유의성¶

import numpy as np
from scipy import stats

# ===================================================================
# 정규성이 깨졌을 때 두 검정이 얼마나 다르게 버티는가
#
# 같은 자료에 평균에 대한 t 검정과 분산에 대한 카이제곱 검정을 함께
# 돌린다. 둘 다 귀무가설이 참인데, t 검정의 오류율은 0.05 근처에
# 머물고 분산 검정의 오류율은 크게 벗어난다.
#
# 평균은 중심극한정리의 보호를 받지만 분산은 그렇지 않기 때문이다.
# "정규성이 필요하다"는 말의 무게가 검정마다 다르다는 뜻이다.
# ===================================================================

np.random.seed(42)
n = 30
alpha = 0.05
n_simulations = 10_000
df_true = 5  # 꼬리가 두꺼운 t 분포

rejections_t = 0
rejections_chi2 = 0

for _ in range(n_simulations):
    sample = stats.t.rvs(df=df_true, size=n)

    # 평균이 0 인지 검정한다. t 분포의 평균은 참으로 0 이므로 귀무가설이 참이다.
    _, p_t = stats.ttest_1samp(sample, popmean=0)
    if p_t < alpha:
        rejections_t += 1

    # 분산이 df/(df-2) 인지 검정한다. 이것이 t 분포의 참 분산이므로
    # 이쪽도 귀무가설이 참이다. 조건을 똑같이 맞춘 것이다.
    true_var = df_true / (df_true - 2)
    chi2_stat = (n - 1) * np.var(sample, ddof=1) / true_var
    p_chi2 = 2 * min(
        stats.chi2.cdf(chi2_stat, df=n - 1),
        1 - stats.chi2.cdf(chi2_stat, df=n - 1)
    )
    if p_chi2 < alpha:
        rejections_chi2 += 1

if __name__ == "__main__":
    print(f"Data: t-distribution with df = {df_true}, n = {n}")
    print(f"Nominal alpha: {alpha}")
    print(f"Actual Type I error (t-test):          "
          f"{rejections_t / n_simulations:.4f}")
    print(f"Actual Type I error (chi-squared test): "
          f"{rejections_chi2 / n_simulations:.4f}")

Data: t-distribution with df = 5, n = 30
Nominal alpha: 0.05
Actual Type I error (t-test):          0.0487
Actual Type I error (chi-squared test): 0.1713

import numpy as np
from scipy import stats

nu, n, alpha = 5, 30, 0.05
sigma2 = nu / (nu - 2)                        # t_5 의 참 분산
mu4 = 3 * nu ** 2 / ((nu - 2) * (nu - 4))     # 4차 적률 (nu > 4 에서만 유한)
g2 = mu4 / sigma2 ** 2 - 3
print(f"t_{nu}:  sigma^2 = {sigma2:.4f}   mu4 = {mu4:.4f}   초과첨도 = {g2:.1f}")

# S^2 의 분산. 카이제곱이 가정하는 값과 실제 값을 나란히 적는다.
var_chi2 = 2 * sigma2 ** 2 / (n - 1)
var_true = mu4 / n - sigma2 ** 2 * (n - 3) / (n * (n - 1))
print("\nVar(S^2):")
print(f"  카이제곱이 가정하는 값  2 sigma^4/(n-1)   = {var_chi2:.6f}")
print(f"  실제값  mu4/n - sigma^4 (n-3)/(n(n-1))    = {var_true:.6f}")
print(f"  근사식  sigma^4 (g2 + 2)/n                = "
      f"{sigma2 ** 2 * (g2 + 2) / n:.6f}")
print(f"  비  {var_true / var_chi2:.3f} 배   (표준편차로는 "
      f"{np.sqrt(var_true / var_chi2):.3f} 배)")

# 그 비로 기각률을 정규근사해 보고, 모의값과 맞는지 본다.
pred = 2 * stats.norm.cdf(-stats.norm.ppf(1 - alpha / 2)
                          * np.sqrt(var_chi2 / var_true))
print(f"\n정규근사 예측 기각률 = {pred:.4f}")

rng = np.random.default_rng(55)
R = 400000
S2 = rng.standard_t(nu, size=(R, n)).var(axis=1, ddof=1)
lo = stats.chi2.ppf(alpha / 2, n - 1) * sigma2 / (n - 1)
hi = stats.chi2.ppf(1 - alpha / 2, n - 1) * sigma2 / (n - 1)
print(f"카이제곱 기각역: S^2 < {lo:.4f} 또는 S^2 > {hi:.4f}")
print(f"  모의 (R = {R}):  아래쪽 {np.mean(S2 < lo):.4f}"
      f"   위쪽 {np.mean(S2 > hi):.4f}   합 {np.mean((S2 < lo) | (S2 > hi)):.4f}")
print(f"  정규 자료라면 각각 {alpha / 2}, {alpha / 2}, 합 {alpha}")

# 예측이 어긋난 까닭: S^2 의 분포가 정규와 거리가 멀다.
# t_5 는 6차 적률이 없으므로 S^2 의 왜도가 유한하지 않고, 모의 왜도는
# 씨앗마다 들쭉날쭉하다. Var(S^2) 의 모의 추정값도 함께 흔들린다.
print("\n씨앗을 바꾸어 본다 (각 R = 400000)")
for seed in (55, 56, 57, 58, 59):
    rng = np.random.default_rng(seed)
    s2 = rng.standard_t(nu, size=(R, n)).var(axis=1, ddof=1)
    print(f"  seed {seed}:  E[S^2] {s2.mean():.4f}"
          f"   모의 Var(S^2) {s2.var():.4f}"
          f"   모의 왜도 {stats.skew(s2):5.1f}   최대 S^2 {s2.max():4.0f}")
print(f"  이론:      E[S^2] {sigma2:.4f}   Var(S^2) {var_true:.4f}   왜도 유한하지 않다")

t_5:  sigma^2 = 1.6667   mu4 = 25.0000   초과첨도 = 6.0

Var(S^2):
  카이제곱이 가정하는 값  2 sigma^4/(n-1)   = 0.191571
  실제값  mu4/n - sigma^4 (n-3)/(n(n-1))    = 0.747126
  근사식  sigma^4 (g2 + 2)/n                = 0.740741
  비  3.900 배   (표준편차로는 1.975 배)

정규근사 예측 기각률 = 0.3210
카이제곱 기각역: S^2 < 0.9222 또는 S^2 > 2.6277
  모의 (R = 400000):  아래쪽 0.0907   위쪽 0.0825   합 0.1732
  정규 자료라면 각각 0.025, 0.025, 합 0.05

씨앗을 바꾸어 본다 (각 R = 400000)
  seed 55:  E[S^2] 1.6647   모의 Var(S^2) 0.7191   모의 왜도  17.7   최대 S^2  145
  seed 56:  E[S^2] 1.6641   모의 Var(S^2) 0.6853   모의 왜도   9.5   최대 S^2   80
  seed 57:  E[S^2] 1.6671   모의 Var(S^2) 0.7156   모의 왜도  11.8   최대 S^2   92
  seed 58:  E[S^2] 1.6689   모의 Var(S^2) 0.7878   모의 왜도  24.1   최대 S^2  160
  seed 59:  E[S^2] 1.6667   모의 Var(S^2) 0.6765   모의 왜도   6.9   최대 S^2   55
  이론:      E[S^2] 1.6667   Var(S^2) 0.7471   왜도 유한하지 않다