정규성 검정의 한계와 함정¶

import numpy as np
from scipy.stats import shapiro

np.random.seed(0)

# 두 자료 모두 참으로 정규다. 그런데 표본크기만 다르게 두면 검정의
# 성질이 달라진다. 작은 표본에서는 웬만한 이탈도 잡아내지 못하고,
# 큰 표본에서는 무시할 만한 이탈에도 기각한다.
small_sample = np.random.normal(0, 1, 20)
stat_small, p_value_small = shapiro(small_sample)
print(f"Shapiro-Wilk Test (small sample): p-value={p_value_small}")

# 표본이 커지면 검정력이 올라가는데, 그것이 늘 좋은 일은 아니다.
large_sample = np.random.normal(0, 1, 10000)
stat_large, p_value_large = shapiro(large_sample)
print(f"Shapiro-Wilk Test (large sample): p-value={p_value_large}")

Shapiro-Wilk Test (small sample): p-value=0.7136380887419513
Shapiro-Wilk Test (large sample): p-value=0.7506816458200554

import warnings

import numpy as np
from scipy import stats

warnings.filterwarnings("ignore")   # n > 5000 에서 shapiro 가 내는 경고

# (1) 보기의 두 표본을 다시 뽑아 p 값만이 아니라 통계량까지 함께 본다.
np.random.seed(0)
small = np.random.normal(0, 1, 20)
large = np.random.normal(0, 1, 10000)
for name, x in [("n=20   ", small), ("n=10000", large)]:
    W, p = stats.shapiro(x)
    print(f"{name}  W = {W:.6f}  p = {p:.4f}"
          f"  g1 = {stats.skew(x):+.4f}  g2 = {stats.kurtosis(x):+.4f}")


def pvals(draw, n, R, seed):
    """같은 분포에서 n 개짜리 표본을 R 번 뽑아 샤피로 p 값을 모은다."""
    rng = np.random.default_rng(seed)
    return np.array([stats.shapiro(draw(rng, n)).pvalue for _ in range(R)])


# 귀무가설이 참이면 p 값은 n 과 무관하게 균등분포이고 기각률은 알파와 같다.
normal = lambda rng, n: rng.normal(0, 1, n)
print("\n귀무(정규) 아래 기각률 — n 과 무관하게 0.05 여야 한다")
for n, R, seed in [(20, 4000, 7), (200, 4000, 8), (1000, 2000, 9), (10000, 1500, 23)]:
    p = pvals(normal, n, R, seed)
    print(f"  n = {n:5d}  기각률 {np.mean(p < 0.05):.4f}"
          f"  (MC SE {np.sqrt(0.05 * 0.95 / R):.4f})  p 중앙값 {np.median(p):.3f}")

# 5% 임계값. 같은 p 값이 n 에 따라 얼마나 다른 적합도를 뜻하는지 본다.
print("\nW 의 5% 임계값")
for n, R, seed in [(20, 20000, 120), (10000, 1500, 10100)]:
    rng = np.random.default_rng(seed)
    W = np.array([stats.shapiro(rng.normal(0, 1, n)).statistic for _ in range(R)])
    print(f"  n = {n:5d}  W 5% 분위수 {np.percentile(W, 5):.6f}"
          f"  -> 1-W 의 허용폭 {1 - np.percentile(W, 5):.6f}")

# (2) 자료를 정규에서 살짝 떼어 놓으면 비로소 n 의 효과가 보인다.
t30 = lambda rng, n: rng.standard_t(30, n)
print(f"\nt_30 의 참 초과첨도 = 6/(30-4) = {6 / 26:.4f}  (n 과 무관하다)")
print("t_30 자료에 대한 검정력")
for n, R, seed in [(30, 3000, 11), (200, 3000, 12), (1000, 1500, 13), (5000, 400, 14)]:
    p = pvals(t30, n, R, seed)
    print(f"  n = {n:5d}  검정력 {np.mean(p < 0.05):.3f}  p 중앙값 {np.median(p):.4f}")

n=20     W = 0.968067  p = 0.7136  g1 = +0.1826  g2 = -0.6244
n=10000  W = 0.999844  p = 0.7507  g1 = +0.0274  g2 = -0.0312

귀무(정규) 아래 기각률 — n 과 무관하게 0.05 여야 한다
  n =    20  기각률 0.0437  (MC SE 0.0034)  p 중앙값 0.500
  n =   200  기각률 0.0450  (MC SE 0.0034)  p 중앙값 0.507
  n =  1000  기각률 0.0525  (MC SE 0.0049)  p 중앙값 0.509
  n = 10000  기각률 0.0360  (MC SE 0.0056)  p 중앙값 0.552

W 의 5% 임계값
  n =    20  W 5% 분위수 0.903472  -> 1-W 의 허용폭 0.096528
  n = 10000  W 5% 분위수 0.999656  -> 1-W 의 허용폭 0.000344

t_30 의 참 초과첨도 = 6/(30-4) = 0.2308  (n 과 무관하다)
t_30 자료에 대한 검정력
  n =    30  검정력 0.072  p 중앙값 0.4752
  n =   200  검정력 0.101  p 중앙값 0.4506
  n =  1000  검정력 0.219  p 중앙값 0.2564
  n =  5000  검정력 0.615  p 중앙값 0.0175

import numpy as np
from scipy.stats import skew, kurtosis, shapiro

np.random.seed(0)

# 감마분포는 오른쪽으로 치우쳐 있다. n=1000 이면 검정력이 충분해
# 이 정도 치우침도 확실하게 기각한다.
skewed_data = np.random.gamma(2, 2, 1000)

# 검정이 기각했을 때 "얼마나" 벗어났는지는 이 값들이 말해 준다.
# 실무에서는 p-값보다 이쪽이 더 쓸모 있다.
print(f"Skewness: {skew(skewed_data)}, Kurtosis: {kurtosis(skewed_data)}")

# Shapiro-Wilk 검정
stat, p_value = shapiro(skewed_data)
print(f"Shapiro-Wilk Test: p-value={p_value}")

Skewness: 1.358332294792333, Kurtosis: 2.4079180806130145
Shapiro-Wilk Test: p-value=7.094914686692259e-25

import warnings

import numpy as np
from scipy import stats

warnings.filterwarnings("ignore")

k, theta, n = 2, 2, 1000
print(f"Gamma({k}, {theta}) 이론값:  왜도 {2 / np.sqrt(k):.4f}   초과첨도 {6 / k:.1f}")

np.random.seed(0)
x = np.random.gamma(k, theta, n)
print(f"표본 왜도      g1 = {stats.skew(x):.4f}   G1 = {stats.skew(x, bias=False):.4f}")
print(f"표본 초과첨도  g2 = {stats.kurtosis(x):.4f}   G2 = {stats.kurtosis(x, bias=False):.4f}")

# 벗어난 폭이 우연인가. 정규 가정의 SE 를 적어 보고, 같은 감마분포에서
# 표집분포를 직접 만들어 그 SE 를 쓸 수 있는지 본다.
print(f"\n정규 가정 SE:  왜도 sqrt(6/n) = {np.sqrt(6 / n):.4f}"
      f"   초과첨도 sqrt(24/n) = {np.sqrt(24 / n):.4f}")

rng = np.random.default_rng(31)
R = 20000
Y = rng.gamma(k, theta, size=(R, n))
g1s, g2s = stats.skew(Y, axis=1), stats.kurtosis(Y, axis=1)
print(f"Gamma(2,2) 실제 표집분포 (R = {R}):")
print(f"  g1  평균 {g1s.mean():.4f}  SD {g1s.std():.4f}"
      f"  95% 범위 [{np.percentile(g1s, 2.5):.3f}, {np.percentile(g1s, 97.5):.3f}]")
print(f"  g2  평균 {g2s.mean():.4f}  SD {g2s.std():.4f}"
      f"  95% 범위 [{np.percentile(g2s, 2.5):.3f}, {np.percentile(g2s, 97.5):.3f}]")
print(f"  관측 g1 은 평균에서 {(stats.skew(x) - g1s.mean()) / g1s.std():+.2f} SD")
print(f"  관측 g2 는 평균에서 {(stats.kurtosis(x) - g2s.mean()) / g2s.std():+.2f} SD")

# 교과서의 왜도 표준오차 공식은 어느 판본의 것인가. 정규 자료에서 확인한다.
se_G1 = lambda m: np.sqrt(6 * m * (m - 1) / ((m - 2) * (m + 1) * (m + 3)))
print("\n왜도 표준오차 공식 sqrt(6n(n-1)/((n-2)(n+1)(n+3))) 은 누구의 것인가"
      "  (정규 자료, R = 60000)")
for m in (20, 100, 1000):
    rng = np.random.default_rng(900 + m)
    Z = rng.normal(0, 1, size=(60000, m))
    print(f"  n = {m:5d}  공식 {se_G1(m):.5f}"
          f"   SD(G1) {stats.skew(Z, axis=1, bias=False).std():.5f}"
          f"   SD(g1) {stats.skew(Z, axis=1).std():.5f}")

Gamma(2, 2) 이론값:  왜도 1.4142   초과첨도 3.0
표본 왜도      g1 = 1.3583   G1 = 1.3604
표본 초과첨도  g2 = 2.4079   G2 = 2.4260

정규 가정 SE:  왜도 sqrt(6/n) = 0.0775   초과첨도 sqrt(24/n) = 0.1549
Gamma(2,2) 실제 표집분포 (R = 20000):
  g1  평균 1.4010  SD 0.1706  95% 범위 [1.124, 1.784]
  g2  평균 2.8993  SD 1.2936  95% 범위 [1.289, 5.971]
  관측 g1 은 평균에서 -0.25 SD
  관측 g2 는 평균에서 -0.38 SD

왜도 표준오차 공식 sqrt(6n(n-1)/((n-2)(n+1)(n+3))) 은 누구의 것인가  (정규 자료, R = 60000)
  n =    20  공식 0.51210   SD(G1) 0.51384   SD(g1) 0.47447
  n =   100  공식 0.24138   SD(G1) 0.24151   SD(g1) 0.23788
  n =  1000  공식 0.07734   SD(G1) 0.07785   SD(g1) 0.07774

import numpy as np
from scipy.stats import shapiro

np.random.seed(0)

# 소득처럼 오른쪽으로 크게 치우친 자료는 정규성 검정이 언제나 기각한다.
# 검정을 돌리기 전에 이미 알 수 있는 일이므로, 이럴 때는 검정이 아니라
# 변환이나 비모수 방법으로 바로 넘어가는 편이 낫다.
income_data = np.random.exponential(scale=50000, size=1000)

# Shapiro-Wilk 검정
stat, p_value = shapiro(income_data)
print(f"Shapiro-Wilk Test on Skewed Data: p-value={p_value}")

Shapiro-Wilk Test on Skewed Data: p-value=4.8578866080562626e-33

import numpy as np
from scipy import stats

# 1) 자료의 단위는 아무 몫도 하지 않는다. W 가 양의 아핀변환에 불변이기 때문이다.
np.random.seed(0)
income = np.random.exponential(scale=50000, size=1000)
for name, y in [("원 단위 (scale=50000)", income),
                ("만원 단위 (/10000)", income / 10000),
                ("표준화 (-mean)/sd", (income - income.mean()) / income.std())]:
    W, p = stats.shapiro(y)
    print(f"{name:24s}  W = {W:.10f}  p = {p:.6e}")

x = np.arange(1.0, 11.0)
print(f"\n자료에 3 을 곱하고 7 을 더해도 W 가 그대로인가: "
      f"{stats.shapiro(x).statistic:.10f} vs {stats.shapiro(3 * x + 7).statistic:.10f}")

# 2) 치우침의 크기는 척도와 무관하게 정해져 있다.
print("\nExponential 이론값: 왜도 2, 초과첨도 6")
print(f"표본 g1 = {stats.skew(income):.4f}   g2 = {stats.kurtosis(income):.4f}")

# 3) 그러므로 결과는 검정을 돌리기 전에 이미 정해져 있다. n 이 얼마부터 그런가.
print("\nExponential 자료에 대한 샤피로 검정력 (alpha = 0.05, R = 20000)")
R = 20000
for n in (5, 10, 20, 30, 50, 100):
    rng = np.random.default_rng(500 + n)
    p = np.array([stats.shapiro(rng.exponential(1.0, n)).pvalue for _ in range(R)])
    print(f"  n = {n:4d}  검정력 {np.mean(p < 0.05):.4f}   놓친 횟수 {np.sum(p >= 0.05):5d}")

원 단위 (scale=50000)        W = 0.8050169741  p = 4.857887e-33
만원 단위 (/10000)            W = 0.8050169741  p = 4.857887e-33
표준화 (-mean)/sd            W = 0.8050169741  p = 4.857887e-33

자료에 3 을 곱하고 7 을 더해도 W 가 그대로인가: 0.9701646111 vs 0.9701646111

Exponential 이론값: 왜도 2, 초과첨도 6
표본 g1 = 2.0526   g2 = 6.4761

Exponential 자료에 대한 샤피로 검정력 (alpha = 0.05, R = 20000)
  n =    5  검정력 0.1601   놓친 횟수 16798
  n =   10  검정력 0.4373   놓친 횟수 11254
  n =   20  검정력 0.8383   놓친 횟수  3235
  n =   30  검정력 0.9677   놓친 횟수   645
  n =   50  검정력 0.9995   놓친 횟수    10
  n =  100  검정력 1.0000   놓친 횟수     0

import numpy as np

np.random.seed(0)

# 자료 자체는 크게 치우쳐 있다.
skewed_data = np.random.gamma(2, 2, 100)

mean_sample = np.mean(skewed_data)
std_error = np.std(skewed_data) / np.sqrt(len(skewed_data))

# 그런데 t 검정이 정규성을 요구하는 대상은 자료가 아니라 표본평균의
# 표집분포다. 중심극한정리 덕분에 n=100 이면 그쪽은 이미 정규에 가깝다.
# 자료의 정규성 검정이 기각되었다고 t 검정을 못 쓰는 것이 아닌 까닭이다.
print(f"Sample mean: {mean_sample}, Standard error: {std_error}")

Sample mean: 4.097224928182177, Standard error: 0.2886128975105421

import numpy as np
from scipy import stats

k, theta, n = 2, 2, 100
mu, sigma = k * theta, np.sqrt(k) * theta
g1_pop, g2_pop = 2 / np.sqrt(k), 6 / k
print(f"Gamma({k},{theta}) 모집단:  mu = {mu}  sigma = {sigma:.4f}"
      f"  왜도 {g1_pop:.4f}  초과첨도 {g2_pop:.1f}")
print(f"표본평균의 이론값:  SE = sigma/sqrt(n) = {sigma / np.sqrt(n):.4f}"
      f"  왜도 = {g1_pop / np.sqrt(n):.4f}  초과첨도 = {g2_pop / n:.4f}")

np.random.seed(0)
x = np.random.gamma(k, theta, n)
print(f"\n보기의 표본:  평균 {x.mean():.4f}"
      f"  SE(ddof=0) {x.std() / np.sqrt(n):.4f}"
      f"  SE(ddof=1) {x.std(ddof=1) / np.sqrt(n):.4f}")

# 표본평균의 표집분포를 직접 만들어 세 이론값을 확인한다.
rng = np.random.default_rng(77)
R = 200000
xbar = rng.gamma(k, theta, size=(R, n)).mean(axis=1)
print(f"\n표집분포 모의 (R = {R}):  SD {xbar.std():.4f}"
      f"  왜도 {stats.skew(xbar):.4f}  초과첨도 {stats.kurtosis(xbar):.4f}")

# t 구간의 포함률. t 검정이 요구하는 것은 자료가 아니라 xbar 의 정규성이다.
print("\nmu 에 대한 95% t 구간의 포함률 (명목 0.95, R = 40000)")
R2 = 40000
for m in (5, 10, 30, 100):
    rng = np.random.default_rng(1000 + m)
    Y = rng.gamma(k, theta, size=(R2, m))
    half = stats.t.ppf(0.975, m - 1) * Y.std(axis=1, ddof=1) / np.sqrt(m)
    cov = np.mean(np.abs(Y.mean(axis=1) - mu) <= half)
    print(f"  n = {m:4d}  포함률 {cov:.4f}  (MC SE {np.sqrt(cov * (1 - cov) / R2):.4f})"
          f"   xbar 의 이론 왜도 {g1_pop / np.sqrt(m):.3f}")

Gamma(2,2) 모집단:  mu = 4  sigma = 2.8284  왜도 1.4142  초과첨도 3.0
표본평균의 이론값:  SE = sigma/sqrt(n) = 0.2828  왜도 = 0.1414  초과첨도 = 0.0300

보기의 표본:  평균 4.0972  SE(ddof=0) 0.2886  SE(ddof=1) 0.2901

표집분포 모의 (R = 200000):  SD 0.2826  왜도 0.1354  초과첨도 0.0210

mu 에 대한 95% t 구간의 포함률 (명목 0.95, R = 40000)
  n =    5  포함률 0.9136  (MC SE 0.0014)   xbar 의 이론 왜도 0.632
  n =   10  포함률 0.9234  (MC SE 0.0013)   xbar 의 이론 왜도 0.447
  n =   30  포함률 0.9386  (MC SE 0.0012)   xbar 의 이론 왜도 0.258
  n =  100  포함률 0.9472  (MC SE 0.0011)   xbar 의 이론 왜도 0.141