왜도 검정과 첨도 검정¶

import numpy as np
from scipy import stats

np.random.seed(0)

# 지금은 정규자료다. 위 주석을 바꿔 감마분포로 돌리면 검정이 기각되는 것을
# 볼 수 있다. 두 경우를 견주는 것이 이 코드의 쓰임이다.
# data = np.random.gamma(2, 2, 1000)
data = np.random.normal(2, 2, 1000)

skewness_value = stats.skew(data)
print(f"Skewness: {skewness_value:.4f}")

# 귀무가설은 "모집단의 왜도가 0" 이다. 이 검정은 치우침만 보므로,
# 대칭이면서 꼬리만 두꺼운 분포는 잡아내지 못한다.
stat, p_value = stats.skewtest(data)
print(f"Skewness Test: Statistic={stat:.4f}, p-value={p_value:.4f}")

# 결과 해석
alpha = 0.05
if p_value <= alpha:
    print("Reject H_0: The data is not symmetrically distributed (significant skewness).")
else:
    print("Fail to reject H_0: The data is symmetrically distributed (no significant skewness).")

Skewness: 0.0339
Skewness Test: Statistic=0.4402, p-value=0.6598
Fail to reject H_0: The data is symmetrically distributed (no significant skewness).

import numpy as np
from scipy import stats
from scipy.optimize import brentq

n = 1000

# skewtest 의 변환을 g1 의 함수로 쓴다 (단조증가이므로 역이 하나다).
def z1_of(g1, n):
    Y = g1 * np.sqrt((n + 1) * (n + 3) / (6.0 * (n - 2)))
    beta2 = (3.0 * (n**2 + 27 * n - 70) * (n + 1) * (n + 3)
             / ((n - 2.0) * (n + 5) * (n + 7) * (n + 9)))
    W2 = -1 + np.sqrt(2 * (beta2 - 1))
    delta = 1 / np.sqrt(0.5 * np.log(W2))
    alpha = np.sqrt(2.0 / (W2 - 1))
    return delta * np.arcsinh(Y / alpha)

g1_cut = brentq(lambda g: z1_of(g, n) - 1.959964, 0.0, 1.0)
print(f"|Z1| = 1.96 이 되는 g1 = {g1_cut:.4f}")
print(f"어림 1.96*sqrt(6/n)   = {1.959964 * np.sqrt(6 / n):.4f}")

print(f"\nGamma(k=2) 이론 왜도 = 2/sqrt(k) = {2 / np.sqrt(2.0):.4f}")
np.random.seed(0)
gam = np.random.gamma(2, 2, 1000)
zg, pg = stats.skewtest(gam)
print(f"감마 표본 g1 = {stats.skew(gam):.4f},  Z1 = {zg:.4f},  p = {pg:.4g}")

# 정규 자료에서의 실제 기각률이 명목을 지키는지도 본다.
rng = np.random.default_rng(42)
R = 4000
rej = np.mean([stats.skewtest(s)[1] < 0.05 for s in rng.standard_normal((R, n))])
print(f"\n정규 자료 기각률 = {rej:.4f}  (명목 0.05, MC SE {np.sqrt(0.05 * 0.95 / R):.4f})")

|Z1| = 1.96 이 되는 g1 = 0.1515
어림 1.96*sqrt(6/n)   = 0.1518

Gamma(k=2) 이론 왜도 = 2/sqrt(k) = 1.4142
감마 표본 g1 = 1.3583,  Z1 = 13.6656,  p = 1.629e-42

정규 자료 기각률 = 0.0530  (명목 0.05, MC SE 0.0034)

import numpy as np
from scipy import stats

np.random.seed(0)

# data = np.random.gamma(2, 2, 1000)
data = np.random.normal(2, 2, 1000)

kurtosis_value = stats.kurtosis(data)
print(f"Kurtosis: {kurtosis_value:.4f}")

# 귀무가설은 "모집단의 초과첨도가 0" 이다. 앞의 왜도 검정과 짝을 이루며,
# 둘을 합친 것이 D'Agostino 의 K^2 이다.
stat, p_value = stats.kurtosistest(data)
print(f"Kurtosis Test: Statistic={stat:.4f}, p-value={p_value:.4f}")

# 결과 해석
alpha = 0.05
if p_value <= alpha:
    print("Reject H_0: The data does not have normal kurtosis.")
else:
    print("Fail to reject H_0: The data has normal kurtosis.")

Kurtosis: -0.0468
Kurtosis Test: Statistic=-0.1980, p-value=0.8431
Fail to reject H_0: The data has normal kurtosis.

import numpy as np
from scipy import stats
from scipy.optimize import brentq

n = 1000

# kurtosistest 의 변환을 b2 의 함수로 쓴다 (단조증가이므로 역이 하나다).
def z2_of(b2, n):
    E = 3.0 * (n - 1) / (n + 1)
    varb2 = 24.0 * n * (n - 2) * (n - 3) / ((n + 1) ** 2 * (n + 3) * (n + 5))
    u = (b2 - E) / np.sqrt(varb2)
    sqrtbeta1 = (6.0 * (n * n - 5 * n + 2) / ((n + 7) * (n + 9))
                 * np.sqrt(6.0 * (n + 3) * (n + 5) / (n * (n - 2) * (n - 3))))
    A = 6.0 + 8.0 / sqrtbeta1 * (2.0 / sqrtbeta1 + np.sqrt(1 + 4.0 / sqrtbeta1**2))
    term1 = 1 - 2 / (9.0 * A)
    term2 = ((1 - 2.0 / A) / (1 + u * np.sqrt(2 / (A - 4.0)))) ** (1 / 3)
    return (term1 - term2) / np.sqrt(2 / (9.0 * A))

lo = brentq(lambda g: z2_of(g + 3.0, n) + 1.959964, -0.5, 0.0)
hi = brentq(lambda g: z2_of(g + 3.0, n) - 1.959964, 0.0, 2.0)
print(f"기각하지 않는 g2 구간 = [{lo:.4f}, {hi:.4f}]")
print(f"폭의 비 (위/아래) = {hi / abs(lo):.3f}")
print(f"대칭 어림 1.96*sqrt(24/n) = +-{1.959964 * np.sqrt(24 / n):.4f}")

print(f"\nGamma(k=2) 이론 초과첨도 = 3(k+2)/k - 3 = 6/k = {6 / 2.0:.4f}")
np.random.seed(0)
gam = np.random.gamma(2, 2, 1000)
z1, p1 = stats.skewtest(gam)
z2, p2 = stats.kurtosistest(gam)
print(f"감마 표본 g2 = {stats.kurtosis(gam):.4f},  Z2 = {z2:.4f},  p = {p2:.4g}")
print(f"감마 표본 Z1 = {z1:.4f} (왜도)  대  Z2 = {z2:.4f} (첨도)")

# 정규 자료에서의 실제 기각률
rng = np.random.default_rng(43)
R = 4000
rej = np.mean([stats.kurtosistest(s)[1] < 0.05 for s in rng.standard_normal((R, n))])
print(f"\n정규 자료 기각률 = {rej:.4f}  (명목 0.05, MC SE {np.sqrt(0.05 * 0.95 / R):.4f})")

기각하지 않는 g2 구간 = [-0.2744, 0.3272]
폭의 비 (위/아래) = 1.192
대칭 어림 1.96*sqrt(24/n) = +-0.3036

Gamma(k=2) 이론 초과첨도 = 3(k+2)/k - 3 = 6/k = 3.0000
감마 표본 g2 = 2.4079,  Z2 = 7.7918,  p = 6.607e-15
감마 표본 Z1 = 13.6656 (왜도)  대  Z2 = 7.7918 (첨도)

정규 자료 기각률 = 0.0470  (명목 0.05, MC SE 0.0034)