첨도 검정¶

import numpy as np
from scipy import stats

# 정규에 자유도 4 짜리 t 를 섞었다. 중심은 그대로이고 꼬리만 두꺼워진다.
# 왜도 검정으로는 잡히지 않고 첨도 검정에만 걸리는 자료다.
rng = np.random.default_rng(0)
x = np.concatenate([rng.normal(0, 1, size=220),
                    rng.standard_t(df=4, size=80)])

# 표본첨도를 바로 쓰지 않고 Z 로 바꾸는 까닭은, 첨도의 표집분포가 표본이
# 상당히 클 때까지도 크게 치우쳐 있기 때문이다. 이 검정은 n>=20 부터 쓸 만하다.
g2 = stats.kurtosis(x, fisher=True, bias=False)
z, p = stats.kurtosistest(x)

print(f"Sample size n = {x.size}")
print(f"Sample excess kurtosis (Fisher) g2 = {g2:.4f}")
print(f"D'Agostino kurtosis test: Z = {z:.4f}, p-value = {p:.4g}")
if p < 0.05:
    print("=> Evidence of non-normal kurtosis (departing from normality).")
else:
    print("=> No strong evidence of non-normal kurtosis.")

Sample size n = 300
Sample excess kurtosis (Fisher) g2 = 1.0493
D'Agostino kurtosis test: Z = 2.7603, p-value = 0.005774
=> Evidence of non-normal kurtosis (departing from normality).

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
x = np.concatenate([rng.normal(0, 1, size=220),
                    rng.standard_t(df=4, size=80)])
n = x.size

b2 = stats.kurtosis(x, fisher=False)              # 비보정, 3 을 빼지 않은 것
g2 = stats.kurtosis(x)                            # 비보정 초과첨도 = b2 - 3
G2 = stats.kurtosis(x, fisher=True, bias=False)   # 보정 초과첨도
print(f"b2 (비보정)      = {b2:.6f}")
print(f"g2 = b2 - 3      = {g2:.6f}")
print(f"G2 (보정)        = {G2:.6f}")

# 1단계: 정규성 아래의 평균과 분산으로 표준화.
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)
print(f"\nE[b2] = 3(n-1)/(n+1) = {E:.6f},  SD(b2) = {np.sqrt(varb2):.6f}")
print(f"단순 표준화 u = {u:.6f}")

# 2단계: 귀무분포의 왜도를 윌슨–힐퍼티 세제곱근으로 펴 준다.
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)
Z = (term1 - term2) / np.sqrt(2 / (9.0 * A))
print(f"sqrt(beta1) = {sqrtbeta1:.6f},  A = {A:.6f}")
print(f"손계산 Z2 = {Z:.6f},  p = {2 * stats.norm.sf(abs(Z)):.6f}")

z_scipy, p_scipy = stats.kurtosistest(x)
print(f"scipy  Z2 = {z_scipy:.6f},  p = {p_scipy:.6f}")
print(f"\n단순 표준화만 했다면 p = {2 * stats.norm.sf(abs(u)):.6f}")

b2 (비보정)      = 4.011938
g2 = b2 - 3      = 1.011938
G2 (보정)        = 1.049279

E[b2] = 3(n-1)/(n+1) = 2.980066,  SD(b2) = 0.275877
단순 표준화 u = 3.740339
sqrt(beta1) = 0.808923,  A = 56.827255
손계산 Z2 = 2.760321,  p = 0.005774
scipy  Z2 = 2.760321,  p = 0.005774

단순 표준화만 했다면 p = 0.000184

import numpy as np
from scipy import stats

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

g2 = stats.kurtosis(x, fisher=True, bias=False)
z, p = stats.kurtosistest(x)

print(f"g2 = {g2:.4f}")
print(f"Z = {z:.4f}, p = {p:.4f}")

g2 = 0.1380
Z = 0.6014, p = 0.5476

import numpy as np
from scipy import stats

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

g2 = stats.kurtosis(x, fisher=True, bias=False)
z, p = stats.kurtosistest(x)

print(f"g2 = {g2:.4f}, Z = {z:.4f}, p = {p:.4g}")

g2 = -1.1951, Z = -12.5358, p = 4.755e-36

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n, reps, alpha = 100, 5000, 0.05

for df in [5, 10]:
    rejections = 0
    for _ in range(reps):
        x = rng.standard_t(df=df, size=n)
        _, p = stats.kurtosistest(x)
        if p < alpha:
            rejections += 1
    print(f"t({df}): power = {rejections / reps:.4f}")

t(5): power = 0.6058
t(10): power = 0.2448