D'Agostino K제곱 검정¶

(1) 첨도 성분이 \(69.8\%\)다. 분해는 정확히 맞아떨어진다.

\(\chi^2_2\)의 생존함수가 \(e^{-x/2}\)이므로

이고 normaltest의 \(p\)값과 일치한다. 두 성분 모두 \(\lvert Z\rvert > 1.96\)을 넘으므로 둘 다 혼자서도 기각하는데, 첨도 쪽이 \(15.50\) 대 \(6.72\)로 \(2.3\)배 크다. 대수정규 오염이 치우침과 두꺼운 꼬리를 동시에 들여왔고 꼬리 쪽 신호가 더 강하다.

(2) 차이는 전부 첨도 항에서 온다. 두 통계량을 나란히 쪼개면

이다. 왜도 항은 \(1\%\) 안에서 같다. \(JB\)의 왜도 항이 \(\frac{n}{6}g_1^2 = (g_1/\sqrt{6/n})^2\), 곧 \(\mathrm{Var}(g_1) \approx 6/n\)을 쓴 단순 표준화의 제곱인데, \(n = 300\)에서 왜도 쪽 변환은 거의 아무것도 바꾸지 않으므로 \(Z_1^2\)과 사실상 같은 값이 된다.

첨도 항은 \(2.63\)배 벌어진다. 까닭은 앞 쪽에서 본 대로 \(b_2\)의 귀무분포가 크게 치우쳐 있다는 데 있다. \(JB\)는 \(g_2\)를 그대로 제곱하므로 그 치우침을 전혀 고려하지 않고, \(K^2\)은 윌슨–힐퍼티 세제곱근 변환으로 치우침을 먼저 펴 준 다음 제곱한다. 두꺼운 꼬리 쪽 이탈에서는 변환이 값을 크게 깎으므로 \(Z_2^2\)이 작아진다.

그러므로 \(JB\)가 크다는 것은 검정력이 좋다는 뜻이 아니다. 통계량이 커진 만큼 그 귀무분포도 함께 부풀어야 \(\chi^2_2\)가 맞는데, \(JB\)는 중간 크기 표본에서 그 보정 없이 \(\chi^2_2\)를 쓴다. 그래서 \(p\)값이 \(4.70\times10^{-11}\)까지 내려가지만, 이 쪽에서 \(K^2\)이 주는 \(1.50\times10^{-5}\)가 더 믿을 만한 수다. 두 \(p\)값의 자릿수 차이 여섯은 증거의 차이가 아니라 보정의 차이다.

import numpy as np
from scipy import stats

# 정규 240개에 로그정규 60개를 섞었다. 겉보기에는 정규 같지만 오른쪽이
# 조금 늘어난 자료다. 검정이 이것을 잡아내는지 본다.
rng = np.random.default_rng(0)
x = np.concatenate([rng.normal(0, 1, size=240),
                    rng.lognormal(0, 0.6, size=60)])

# K^2 과 그것을 이루는 두 성분을 함께 구한다. 어느 쪽이 큰지를 보면
# 정규에서 벗어난 까닭이 치우침인지 꼬리인지 알 수 있다.
K2, p = stats.normaltest(x)
z1, _ = stats.skewtest(x)
z2, _ = stats.kurtosistest(x)

print(f"Sample size n = {x.size}")
print(f"D'Agostino's K^2 statistic = {K2:.4f}")
print(f"  components: Z1 = {z1:.4f}, Z2 = {z2:.4f}")
print(f"p-value = {p:.4g}")
if p < 0.05:
    print("=> Reject normality at alpha = 0.05.")
else:
    print("=> Fail to reject normality at alpha = 0.05.")

출력:

Sample size n = 300
D'Agostino's K^2 statistic = 22.2123
  components: Z1 = 2.5915, Z2 = 3.9365
p-value = 1.502e-05
=> Reject normality at alpha = 0.05.

두 통계량의 분해를 나란히 확인한다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
x = np.concatenate([rng.normal(0, 1, size=240),
                    rng.lognormal(0, 0.6, size=60)])
n = x.size
K2, p = stats.normaltest(x)
z1, _ = stats.skewtest(x)
z2, _ = stats.kurtosistest(x)
print(f"Z1^2 = {z1**2:.4f},  Z2^2 = {z2**2:.4f},  합 = {z1**2 + z2**2:.4f}  (K2 = {K2:.4f})")
print(f"첨도 성분의 비중 = {100 * z2**2 / K2:.1f}%")
print(f"p = exp(-K2/2) = {np.exp(-K2 / 2):.4g},  normaltest p = {p:.4g}")

g1, g2 = stats.skew(x), stats.kurtosis(x)
jb, pjb = stats.jarque_bera(x)
jb_s, jb_k = n / 6 * g1**2, n / 24 * g2**2
print(f"\nJB = {jb:.4f} = {jb_s:.4f} + {jb_k:.4f},  첨도 비중 {100 * jb_k / jb:.1f}%")
print(f"  p(JB) = {pjb:.4g}")
print(f"\n왜도 쪽:  JB 항 {jb_s:.4f}  대  Z1^2 {z1**2:.4f}   (비 {jb_s / z1**2:.2f})")
print(f"첨도 쪽:  JB 항 {jb_k:.4f}  대  Z2^2 {z2**2:.4f}   (비 {jb_k / z2**2:.2f})")

출력:

Z1^2 = 6.7159,  Z2^2 = 15.4964,  합 = 22.2123  (K2 = 22.2123)
첨도 성분의 비중 = 69.8%
p = exp(-K2/2) = 1.502e-05,  normaltest p = 1.502e-05

JB = 47.5606 = 6.8005 + 40.7601,  첨도 비중 85.7%
  p(JB) = 4.703e-11

왜도 쪽:  JB 항 6.8005  대  Z1^2 6.7159   (비 1.01)
첨도 쪽:  JB 항 40.7601  대  Z2^2 15.4964   (비 2.63)

분해가 양쪽 모두 끝까지 맞고, \(e^{-K^2/2}\)가 normaltest의 \(p\)값과 일치한다. 마지막 두 줄이 (2)의 답이다. 같은 자료에서 두 옴니버스 검정이 왜도에는 똑같이 반응하고 첨도에만 다르게 반응한다. \(\square\)

import numpy as np
from scipy import stats

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

K2, p = stats.normaltest(x)
print(f"K^2 = {K2:.4f}, p = {p:.4f}")

출력:

K^2 = 5.9058, p = 0.0522

\(p = 0.0522\)로 0.05를 간신히 넘었다. 형식적으로는 기각하지 못하지만 문턱에 아슬아슬하게 걸쳐 있다.

자료를 정확히 표준정규에서 생성했으므로 이것이 우연임을 우리는 안다. 귀무가설 아래에서 \(p\)값은 \(\text{Uniform}(0,1)\)을 따르므로 0.05 근처의 값이 나올 확률도 다른 구간과 똑같다.

이 예는 \(\alpha = 0.05\)라는 문턱의 자의성을 잘 보여준다. \(p = 0.0522\)와 \(p = 0.0478\)은 증거로서 사실상 구별되지 않는데도 이분법적 판정은 정반대가 된다. \(p\)값을 그대로 보고하고 문턱 통과 여부만 말하지 않는 편이 낫다. \(\square\)

import numpy as np
from scipy import stats

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

K2, p = stats.normaltest(x)
z1, p1 = stats.skewtest(x)
z2, p2 = stats.kurtosistest(x)

print(f"K^2 = {K2:.4f}, p = {p:.4g}")
print(f"Skewness:  Z1 = {z1:.4f}, p = {p1:.4g}")
print(f"Kurtosis:  Z2 = {z2:.4f}, p = {p2:.4g}")

출력:

K^2 = 194.4187, p = 6.06e-43
Skewness:  Z1 = 0.2788, p = 0.7804
Kurtosis:  Z2 = -13.9406, p = 3.588e-44

균등분포는 대칭이고(\(\gamma_1 = 0\)) 저첨이다(\(\gamma_2 = -1.2\)). 따라서 \(Z_1 = 0.279\)는 작고 유의하지 않은 반면 \(Z_2 = -13.94\)는 크고 매우 유의하다.

분해를 보면 기여가 명확하다.

왜도 성분의 기여가 전체의 \(0.04\%\)에 불과하다. \(K^2\)의 기각은 전적으로 첨도 성분이 이끈다. \(\square\)

\(H_0\) 아래에서 D'Agostino의 변환은 \(g_1\)과 \(g_2\)를 각각 근사적으로 \(\mathcal{N}(0,1)\)인 \(Z_1\)과 \(Z_2\)로 보낸다. 이 변환들은 \(Z_1\)과 \(Z_2\)가 근사적으로 독립이 되도록 설계되어 있다(하나는 홀수 차수 중심적률에, 다른 하나는 짝수 차수에 의존한다). 독립인 두 표준정규의 제곱합은 정의상 \(\chi^2_2\)이다.

근사가 성립하려면 다음이 필요하다.

1. \(Z_1\)과 \(Z_2\)의 정규근사가 성립할 만큼 \(n\)이 커야 한다(SciPy는 \(Z_1\)에 \(n \geq 8\), \(Z_2\)에 \(n \geq 20\)을 요구한다).
2. 자료가 i.i.d.여야 한다. 자료에 계열상관이 있으면 \(g_1\)과 \(g_2\)의 분산이 달라져 \(\chi^2_2\) 근사가 무너진다.

독립성 가정도 근사에 불과하다는 점을 유념하라. 유한표본에서 \(Z_1\)과 \(Z_2\) 사이에는 약한 상관이 남아 있으며, 이것이 연습문제 4에서 관찰할 크기 왜곡의 한 원인이다. \(\square\)

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.normaltest(x)
    if p < alpha:
        rejections += 1

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

출력:

Empirical size: 0.0571

추정값의 표준오차는 \(\sqrt{0.05 \times 0.95 / 10000} = 0.0022\)이다. 관측된 \(0.0571\)은 명목값 \(0.05\)보다 \(3.2\) 표준오차 위에 있다. 우연으로 설명하기 어렵다.

반복수를 20,000회로 늘려 여러 \(n\)에서 확인하면 체계적인 경향이 드러난다.

(SE \(= 0.0015\).)

곧 \(K^2\) 검정은 작은 표본에서 약간 자유주의적(liberal)이다. 실제 제1종 오류율이 명목값보다 높다. \(n\)이 커지면서 차이가 줄어들어 \(n = 500\)에서는 무시할 만해진다.

원인은 두 가지이다. 하나는 \(Z_1\), \(Z_2\) 각각의 정규근사가 작은 \(n\)에서 불완전하다는 것이고, 다른 하나는 연습문제 3에서 언급한 대로 유한표본에서 두 통계량이 완전히 독립이 아니라는 것이다.

실용적 함의: \(n = 50\)에서 \(\alpha = 0.05\)로 \(K^2\) 검정을 쓰면 실제 크기는 약 0.058이다. 크기 조정이 결정적으로 중요하다면 모의실험으로 임계값을 얻거나, 작은 표본에서 크기가 더 정확한 Shapiro-Wilk를 쓰는 편이 낫다. \(\square\)

\(Z_1, Z_2 \overset{\text{iid}}{\sim} \mathcal{N}(0,1)\)이면 정의상 \(Z_1^2 \sim \chi^2_1\)이고 \(Z_2^2 \sim \chi^2_1\)이다. 두 변수가 독립이므로 \(K^2 = Z_1^2 + Z_2^2 \sim \chi^2_2\)이다. \(\chi^2_2\) 분포의 밀도는

이며 이는 \(\text{Exponential}(1/2)\)의 밀도이다. 적분하면

주장이 확인되었다. 자유도 2인 카이제곱분포가 평균 2인 지수분포와 같다는 사실은 유용한 특수 성질이다.

한 가지 따름 결과로, \(K^2\) 검정의 \(p\)값을 닫힌 형태로 계산할 수 있다.

본문 보기로 확인해 보자. \(K^2 = 22.2123\)이므로 \(p = e^{-11.106} = 1.50 \times 10^{-5}\)이고, SciPy가 보고한 \(1.502 \times 10^{-5}\)와 일치한다. \(\square\)