Anderson-Darling 검정¶

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)

# scipy 의 anderson 은 p-값 대신 유의수준별 기각값을 돌려준다.
# 통계량이 기각값보다 크면 기각이다.
result = stats.anderson(data)
statistic = result.statistic
print(f"Anderson-Darling Test: Statistic={statistic}")

# 기각값은 유의수준이 낮아질수록 커진다. 통계량 하나로 여러 수준의
# 판정을 한꺼번에 읽을 수 있다.
for significance_level, critical_value in zip(result.significance_level, result.critical_values):
    if statistic >= critical_value:
        print(f"At {significance_level}% significance level: Reject H_0. The data is not normally distributed.")
    else:
        print(f"At {significance_level}% significance level: Fail to reject H_0. The data is normally distributed.")

Anderson-Darling Test: Statistic=0.24321791746319832
At 15.0% significance level: Fail to reject H_0. The data is normally distributed.
At 10.0% significance level: Fail to reject H_0. The data is normally distributed.
At 5.0% significance level: Fail to reject H_0. The data is normally distributed.
At 2.5% significance level: Fail to reject H_0. The data is normally distributed.
At 1.0% significance level: Fail to reject H_0. The data is normally distributed.

import numpy as np
from scipy import stats

np.random.seed(0)
a = np.random.normal(1, 10, 1000)   # 보기 1·3 의 자료
np.random.seed(0)
b = np.random.normal(0, 1, 1000)    # 보기 2 의 자료
print(f"a 가 10b + 1 인가: {np.array_equal(a, 10 * b + 1)}")

ra, rb = stats.anderson(a), stats.anderson(b)
print(f"A2(a) = {ra.statistic:.16f}")
print(f"A2(b) = {rb.statistic:.16f}")
print(f"두 값의 차 = {abs(ra.statistic - rb.statistic):.3e}")

# 정의대로 손계산. 표준화에 쓰는 표준편차가 ddof=1 인지 0 인지가 갈린다.
def A2(x, ddof):
    z = np.sort((x - x.mean()) / x.std(ddof=ddof))
    n = len(z)
    F = stats.norm.cdf(z)
    i = np.arange(1, n + 1)
    return -n - np.sum((2 * i - 1) * (np.log(F) + np.log(1 - F[::-1]))) / n

print(f"\n손계산 (ddof=1) = {A2(b, 1):.16f}")
print(f"손계산 (ddof=0) = {A2(b, 0):.16f}")

# 임계값의 출처: Stephens 의 표본크기 보정
base = np.array([0.576, 0.656, 0.787, 0.918, 1.092])
for nn in (20, 50, 100, 1000):
    print(f"n = {nn:>4}: {np.round(base / (1 + 4.0 / nn - 25.0 / nn**2), 3)}")
print(f"scipy 가 준 임계값 (n = 1000): {ra.critical_values}")

a 가 10b + 1 인가: True
A2(a) = 0.2432179174631983
A2(b) = 0.2432179174634257
두 값의 차 = 2.274e-13

손계산 (ddof=1) = 0.2432179174633120
손계산 (ddof=0) = 0.2425011779185979
n =   20: [0.506 0.577 0.692 0.807 0.96 ]
n =   50: [0.538 0.613 0.736 0.858 1.021]
n =  100: [0.555 0.632 0.759 0.885 1.053]
n = 1000: [0.574 0.653 0.784 0.914 1.088]
scipy 가 준 임계값 (n = 1000): [0.574 0.653 0.784 0.914 1.088]

import numpy as np
from statsmodels.stats.diagnostic import normal_ad

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

# statsmodels 의 normal_ad 는 같은 통계량에 p-값까지 붙여 준다.
# 기각값 표를 읽는 대신 p-값으로 바로 판단하고 싶을 때 쓴다.
statistic, p_value = normal_ad(data)
print(f"Anderson-Darling Test: Statistic={statistic}, p-value={p_value}")

Anderson-Darling Test: Statistic=0.2432179174634257, p-value=0.7659878263029309

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

np.random.seed(0)
data = np.random.normal(0, 1, 1000)
n = data.size
ad2, p = normal_ad(data)
print(f"normal_ad:        A2 = {ad2:.16f},  p = {p:.16f}")
print(f"scipy.anderson:   A2 = {stats.anderson(data).statistic:.16f}")

# 정의대로 손계산 — 표준화에 ddof=1 을 쓴다.
z = np.sort((data - data.mean()) / data.std(ddof=1))
F = stats.norm.cdf(z)
i = np.arange(1, n + 1)
A2_hand = -n - np.sum((2 * i - 1) * (np.log(F) + np.log(1 - F[::-1]))) / n
print(f"손계산 (ddof=1):  A2 = {A2_hand:.16f}")

# statsmodels 의 p 값: 표본크기 보정 뒤 구간별 근사식
ad2a = ad2 * (1 + 0.75 / n + 2.25 / n**2)
print(f"\n보정된 A2* = {ad2a:.10f}  ->  구간 [0.200, 0.340)")
p_hand = 1 - np.exp(-8.318 + 42.796 * ad2a - 59.938 * ad2a**2)
print(f"손계산 p = 1 - exp(-8.318 + 42.796 A2* - 59.938 A2*^2) = {p_hand:.16f}")
print(f"normal_ad p                                            = {p:.16f}")

normal_ad:        A2 = 0.2432179174634257,  p = 0.7659878263029309
scipy.anderson:   A2 = 0.2432179174634257
손계산 (ddof=1):  A2 = 0.2432179174633120

보정된 A2* = 0.2434008781  ->  구간 [0.200, 0.340)
손계산 p = 1 - exp(-8.318 + 42.796 A2* - 59.938 A2*^2) = 0.7659878263029309
normal_ad p                                            = 0.7659878263029309

import numpy as np
from statsmodels.stats.diagnostic import normal_ad

np.random.seed(0)

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

# Anderson-Darling 은 꼬리 쪽 이탈에 특히 민감하다. 꼬리가 문제가 되는
# 금융 자료에서 이 검정을 즐겨 쓰는 까닭이다.
statistic, p_value = normal_ad(data)
print(f"Anderson-Darling Test: Statistic={statistic}, 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.")

Anderson-Darling Test: Statistic=0.243217917463312, p-value=0.7659878263032931
Fail to reject H_0: The data is normally distributed.

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

np.random.seed(0)
data = np.random.normal(1, 10, 1000)   # 이 보기가 실제로 쓰는 자료
n = data.size
print("이 보기의 자료 — 치우쳐 있는가?")
print(f"  g1 = {stats.skew(data):+.4f},  G1 = {stats.skew(data, bias=False):+.4f}")
print(f"  g2 = {stats.kurtosis(data):+.4f}")
print(f"  skewtest p = {stats.skewtest(data)[1]:.4f}")
ad2, p = normal_ad(data)
print(f"  A2 = {ad2:.10f},  p = {p:.10f}")

# 정말로 치우친 자료를 넣으면
np.random.seed(0)
sk = np.random.lognormal(0, 0.6, 1000)
ad2s, ps = normal_ad(sk)
ad2a = ad2s * (1 + 0.75 / n + 2.25 / n**2)
print(f"\n대수정규 자료 (sigma = 0.6)")
print(f"  g1 = {stats.skew(sk):+.4f}   (이론 {(np.exp(0.36) + 2) * np.sqrt(np.exp(0.36) - 1):+.4f})")
print(f"  A2 = {ad2s:.4f},  보정된 A2* = {ad2a:.4f}  ->  네 번째 구간 [0.600, 13]")
print(f"  근사식 p = {np.exp(1.2937 - 5.709 * ad2a + 0.0186 * ad2a**2):.6g}")
print(f"  normal_ad p = {ps:.6g}")
print(f"  scipy 임계값 5% = {stats.anderson(sk).critical_values[2]:.3f}"
      f"  (A2 가 {ad2s / stats.anderson(sk).critical_values[2]:.1f} 배)")

이 보기의 자료 — 치우쳐 있는가?
  g1 = +0.0339,  G1 = +0.0339
  g2 = -0.0468
  skewtest p = 0.6598
  A2 = 0.2432179175,  p = 0.7659878263

대수정규 자료 (sigma = 0.6)
  g1 = +1.8107   (이론 +2.2601)
  A2 = 38.0734,  보정된 A2* = 38.1020  ->  네 번째 구간 [0.600, 13]
  근사식 p = 6.59708e-83
  normal_ad p = 0
  scipy 임계값 5% = 0.784  (A2 가 48.6 배)