잔차의 정규성 확인¶

import numpy as np
import pandas as pd
from scipy import special
from statsmodels.formula.api import ols

# 이 페이지의 모든 진단은 아래 모형 하나를 놓고 수행한다.
# 집단마다 표준편차를 1.0, 1.3, 1.6으로 다르게 주어 진단이 무엇을 잡아내는지
# (그리고 무엇을 못 잡아내는지) 볼 수 있게 했다.
rng = np.random.default_rng(42)
n = 20
data = pd.DataFrame({
    "group": np.repeat(["A", "B", "C"], n),
    "response": np.concatenate([
        rng.normal(10.0, 1.0, n),
        rng.normal(10.8, 1.3, n),
        rng.normal(12.0, 1.6, n),
    ]),
})
model = ols("response ~ C(group)", data=data).fit()

print(data.groupby("group").response.agg(["count", "mean", "std"]).round(3))
print(f"\nF = {model.fvalue:.4f}, p = {model.f_pvalue:.4f}")

# (1) 에서 유도한 c4 와 SD(S)/sigma 를 n = 20 에서 계산한다.
c4 = np.sqrt(2 / (n - 1)) * special.gamma(n / 2) / special.gamma((n - 1) / 2)
rel_sd = np.sqrt(1 - c4 ** 2)
print(f"\nc4 = {c4:.5f}   (치우침 {100 * (1 - c4):.2f}%)")
print(f"SD(S)/sigma = {rel_sd:.5f}   (표집 변동 {100 * rel_sd:.2f}%)")

print(f"\n{'sigma':>7}{'E[S]':>9}{'SD(S)':>9}{'관측 s':>9}{'z':>8}")
obs = data.groupby("group").response.std().values
for sigma, s in zip([1.0, 1.3, 1.6], obs):
    print(f"{sigma:>7.1f}{c4 * sigma:>9.4f}{rel_sd * sigma:>9.4f}"
          f"{s:>9.4f}{(s - c4 * sigma) / (rel_sd * sigma):>+8.3f}")

       count    mean    std
group                      
A         20   9.967  0.870
B         20  10.942  1.034
C         20  12.191  1.145

F = 23.7708, p = 0.0000

c4 = 0.98693   (치우침 1.31%)
SD(S)/sigma = 0.16112   (표집 변동 16.11%)

  sigma     E[S]    SD(S)     관측 s       z
    1.0   0.9869   0.1611   0.8702  -0.725
    1.3   1.2830   0.2095   1.0335  -1.191
    1.6   1.5791   0.2578   1.1450  -1.684

import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from scipy import stats

# 그림에서 읽으려는 것을 먼저 수로 적어 둔다.
e = np.sort(model.resid.values)
N = len(e)
q = stats.norm.ppf(np.arange(1, N + 1) / (N + 1))   # sm.qqplot 의 기본 위치
m, s = e.mean(), e.std(ddof=1)
line = m + s * q                                     # line='s' 가 그리는 기준선

print(f"n = {N},  잔차 표준편차 s = {s:.4f}")
print(f"정렬 잔차와 이론 분위수의 상관 r = {np.corrcoef(e, q)[0, 1]:.5f}   (r^2 = {np.corrcoef(e, q)[0, 1] ** 2:.5f})")
print(f"Shapiro-Wilk W                  = {stats.shapiro(model.resid).statistic:.5f}")
k = np.argmax(np.abs(e - line))
print(f"기준선에서 가장 먼 점: 순위 {k + 1},  관측 {e[k]:+.4f},  기준선 {line[k]:+.4f}"
      f"  (차이 {e[k] - line[k]:+.4f} = {abs(e[k] - line[k]) / s:.2f}s)")
print(f"왜도 {stats.skew(e):+.4f},  초과첨도 {stats.kurtosis(e):+.4f}")

# Q-Q 그림이 가리는 것: 집단마다 잔차의 흩어짐이 다른지는 보이지 않는다.
g = model.model.data.frame["group"]
print("\n집단별 잔차 표준편차 (Q-Q 그림에는 드러나지 않는다)")
print(model.resid.groupby(g).std().round(4).to_string())

# line='s'는 표본의 평균과 표준편차로 정한 기준선이다.
# line='45'(y=x)는 잔차가 표준화되어 있을 때만 맞으므로 여기서는 쓰지 않는다.
sm.qqplot(model.resid, line='s')
plt.title("Q-Q Plot of Residuals")
plt.show()

n = 60,  잔차 표준편차 s = 1.0050
정렬 잔차와 이론 분위수의 상관 r = 0.99115   (r^2 = 0.98238)
Shapiro-Wilk W                  = 0.98696
기준선에서 가장 먼 점: 순위 60,  관측 +2.6418,  기준선 +2.1454  (차이 +0.4965 = 0.49s)
왜도 +0.0390,  초과첨도 -0.0756

집단별 잔차 표준편차 (Q-Q 그림에는 드러나지 않는다)
group
A    0.8702
B    1.0335
C    1.1450

import numpy as np
from scipy.stats import shapiro
from statsmodels.stats.stattools import durbin_watson

# 잔차 전체를 한 번에 넣는다. 집단별로 따로 검정하면 다중검정 문제가 생기고,
# 분산분석이 요구하는 것도 "각 집단의 잔차"가 아니라 하나의 오차 분포다.
stat, p_value = shapiro(model.resid)
print(f"Shapiro-Wilk Test: W = {stat:.4f}, p-value = {p_value:.4f}")

# 같은 잔차를 섞어서 1000번 다시 돌린다. 순서를 바꾸면 W 가 바뀌는가.
e = model.resid.values
rng_sh = np.random.default_rng(7)
perms = [rng_sh.permutation(e) for _ in range(1000)]
W = np.array([shapiro(v).statistic for v in perms])
d = np.array([durbin_watson(v) for v in perms])

print(f"\n원래 순서의 W = {stat!r}")
print(f"섞은 1000개의 W 최소 = {W.min()!r}")
print(f"               W 최대 = {W.max()!r}")
print(f"               폭 = {W.max() - W.min():.3e}  (상대 {(W.max() - W.min()) / stat:.2e})")
print(f"정렬한 잔차의 W = {shapiro(np.sort(e)).statistic!r}")
print(f"뒤집은 잔차의 W = {shapiro(e[::-1]).statistic!r}")

print(f"\n비교 — 순서를 보는 통계량은 이렇게 움직인다")
print(f"원래 순서의 d = {durbin_watson(e):.4f}")
print(f"섞은 1000개의 d 최소 {d.min():.4f}, 최대 {d.max():.4f}, 폭 {d.max() - d.min():.4f}")

Shapiro-Wilk Test: W = 0.9870, p-value = 0.7711

원래 순서의 W = 0.9869641022270246
섞은 1000개의 W 최소 = 0.9869641022270237
               W 최대 = 0.9869641022270246
               폭 = 8.882e-16  (상대 9.00e-16)
정렬한 잔차의 W = 0.9869641022270246
뒤집은 잔차의 W = 0.9869641022270246

비교 — 순서를 보는 통계량은 이렇게 움직인다
원래 순서의 d = 2.1101
섞은 1000개의 d 최소 1.1633, 최대 2.6866, 폭 1.5233

import numpy as np
from scipy import stats

# 표본크기에 따라 검정력이 어떻게 달라지는가. 이론값은 정규일 때 0.05 다.
rng_n = np.random.default_rng(2026)
B = 4000
print(f"반복 {B}회, 명목수준 0.05, 몬테카를로 오차 "
      f"{np.sqrt(0.05 * 0.95 / B):.4f}")

cases = [
    ("정규 (가정 성립)", lambda m: rng_n.normal(0, 1, m)),
    ("t_10 (완만한 이탈)", lambda m: rng_n.standard_t(10, m)),
    ("지수 (심한 이탈)", lambda m: rng_n.exponential(1.0, m)),
]
ns = [10, 25, 60, 200, 800]
print("\n" + " ".join(f"{'n=' + str(m):>8s}" for m in ns) + "  자료 분포")
for name, draw in cases:
    row = []
    for m in ns:
        hit = sum(stats.shapiro(draw(m)).pvalue < 0.05 for _ in range(B))
        row.append(hit / B)
    print(" ".join(f"{v:8.4f}" for v in row) + f"  {name}")

반복 4000회, 명목수준 0.05, 몬테카를로 오차 0.0034

    n=10     n=25     n=60    n=200    n=800  자료 분포
  0.0503   0.0493   0.0568   0.0575   0.0535  정규 (가정 성립)
  0.0747   0.1108   0.1605   0.3510   0.8285  t_10 (완만한 이탈)
  0.4457   0.9260   1.0000   1.0000   1.0000  지수 (심한 이탈)

import matplotlib.pyplot as plt
import numpy as np

# 그림에서 읽으려는 것을 먼저 수로 적어 둔다.
e = model.resid.values
print(f"n = {len(e)},  표준편차 {e.std(ddof=1):.4f},  "
      f"범위 [{e.min():.4f}, {e.max():.4f}]")

print(f"\n{'구간 수':>7s} {'구간 폭':>8s} {'구간당 평균':>10s} {'빈 구간':>7s} "
      f"{'봉우리':>6s}  구간별 도수")
for nb in [8, 12, 20, 30]:
    cnt, edges = np.histogram(e, bins=nb)
    peaks = sum(1 for j in range(nb)
                if cnt[j] > 0
                and (j == 0 or cnt[j] > cnt[j - 1])
                and (j == nb - 1 or cnt[j] >= cnt[j + 1]))
    print(f"{nb:7d} {edges[1] - edges[0]:8.3f} {len(e) / nb:10.2f} "
          f"{int((cnt == 0).sum()):7d} {peaks:6d}  {cnt}")

# 히스토그램은 Q-Q 그림보다 거칠지만 치우침과 봉우리 수를 한눈에 보여 준다.
# 계급 수에 따라 모양이 달라지므로 단독으로 판단하지 말고 Q-Q 그림과 함께 본다.
plt.hist(model.resid, bins=20, density=True, alpha=0.7, edgecolor='black')
plt.xlabel("Residuals")
plt.ylabel("Density")
plt.title("Histogram of Residuals")
plt.show()

n = 60,  표준편차 1.0050,  범위 [-2.5224, 2.6418]

   구간 수     구간 폭     구간당 평균    빈 구간    봉우리  구간별 도수
      8    0.646       7.50       0      2  [ 2  4 15  8 16 11  2  2]
     12    0.430       5.00       0      3  [ 1  1  4 10  6  7 11 11  5  2  1  1]
     20    0.258       3.00       3      5  [1 0 1 1 3 5 5 5 4 4 6 7 6 5 3 2 0 0 1 1]
     30    0.172       2.00       8     10  [1 0 0 1 0 1 0 6 2 5 4 1 3 1 4 5 5 3 3 6 2 2 2 1 0 0 0 1 0 1]

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(4141)
B = 6_000
dists = {
    "정규 (H0 참)": lambda n: rng.normal(0, 1, n),
    "t(10) 약한 두꺼움": lambda n: rng.standard_t(10, n),
    "t(5)": lambda n: rng.standard_t(5, n),
    "균등": lambda n: rng.uniform(-1, 1, n),
    "지수 (치우침 2)": lambda n: rng.exponential(1, n),
    "로그정규 (치우침 6)": lambda n: np.exp(rng.normal(0, 1, n)),
}
ns = [8, 15, 30, 50, 100]
print("샤피로-윌크의 기각률 (명목 0.05)")
print(f"{'분포':>20s} " + " ".join(f"{'n=' + str(n):>8s}" for n in ns))
for lab, f in dists.items():
    row = [sum(stats.shapiro(f(n)).pvalue < 0.05 for _ in range(B)) / B
           for n in ns]
    print(f"{lab:>20s} " + " ".join(f"{r:8.4f}" for r in row))

샤피로-윌크의 기각률 (명목 0.05)
                  분포      n=8     n=15     n=30     n=50    n=100
           정규 (H0 참)   0.0497   0.0522   0.0515   0.0493   0.0507
        t(10) 약한 두꺼움   0.0717   0.0893   0.1188   0.1567   0.2363
                t(5)   0.0993   0.1555   0.2508   0.3558   0.5720
                  균등   0.0690   0.1247   0.3840   0.7545   0.9952
          지수 (치우침 2)   0.3400   0.6715   0.9668   0.9995   1.0000
        로그정규 (치우침 6)   0.4793   0.8273   0.9925   1.0000   1.0000

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(4242)
B = 8_000
print("k=3, n=15, 각 집단은 완벽한 정규. 집단 평균만 δ 씩 떨어져 있다.")
print(f"{'집단 평균 차이 δ':>16s} {'전체 원자료':>11s} {'잔차':>9s}")
for d in [0.0, 1.0, 2.0, 4.0]:
    a = b = 0
    for _ in range(B):
        gs = [rng.normal(i * d, 1, 15) for i in range(3)]
        y = np.concatenate(gs)
        r = np.concatenate([g - g.mean() for g in gs])
        a += stats.shapiro(y).pvalue < 0.05
        b += stats.shapiro(r).pvalue < 0.05
    print(f"{d:16.1f} {a / B:11.4f} {b / B:9.4f}")

k=3, n=15, 각 집단은 완벽한 정규. 집단 평균만 δ 씩 떨어져 있다.
      집단 평균 차이 δ      전체 원자료        잔차
             0.0      0.0471    0.0491
             1.0      0.0356    0.0522
             2.0      0.1042    0.0503
             4.0      0.8828    0.0467

import numpy as np
from scipy import stats

rng = np.random.default_rng(4747)
dists = {
    "정규": lambda n: rng.normal(0, 1, n),
    "t(5) 두꺼운 꼬리": lambda n: rng.standard_t(5, n),
    "균등 (가벼운 꼬리)": lambda n: rng.uniform(-1, 1, n),
    "지수 (오른쪽 치우침)": lambda n: rng.exponential(1, n),
    "좌우 뒤집은 지수": lambda n: -rng.exponential(1, n),
    "이산 (0/1/2)": lambda n: rng.choice([0.0, 1, 2], n, p=[.3, .4, .3]),
}
qs = [0.01, 0.05, 0.50, 0.95, 0.99]
print("표본 20만, 표준화 후 분위수 (Q-Q 그림의 세로축에 해당)")
print(f"{'분포':>20s} {'치우침':>7s} {'초과첨도':>8s} "
      + " ".join(f"{'q' + str(int(q * 100)):>8s}" for q in qs))
print(f"{'표준정규 기준':>20s} {0.0:7.2f} {0.0:8.2f} "
      + " ".join(f"{stats.norm.ppf(q):8.3f}" for q in qs))
for lab, f in dists.items():
    x = f(200_000)
    z = (x - x.mean()) / x.std(ddof=1)
    print(f"{lab:>20s} {stats.skew(x):7.2f} {stats.kurtosis(x):8.2f} "
          + " ".join(f"{np.quantile(z, q):8.3f}" for q in qs))

표본 20만, 표준화 후 분위수 (Q-Q 그림의 세로축에 해당)
                  분포     치우침     초과첨도       q1       q5      q50      q95      q99
             표준정규 기준    0.00     0.00   -2.326   -1.645    0.000    1.645    2.326
                  정규   -0.01     0.01   -2.334   -1.649    0.001    1.648    2.329
         t(5) 두꺼운 꼬리    0.02     4.23   -2.608   -1.560    0.001    1.562    2.629
         균등 (가벼운 꼬리)    0.00    -1.20   -1.693   -1.557    0.002    1.557    1.694
        지수 (오른쪽 치우침)    2.00     5.89   -0.989   -0.948   -0.308    2.000    3.622
           좌우 뒤집은 지수   -1.98     5.77   -3.594   -1.999    0.307    0.949    0.990
          이산 (0/1/2)    0.00    -1.33   -1.288   -1.288    0.003    1.293    1.293

import numpy as np
from scipy import stats

rng = np.random.default_rng(4343)
B = 6_000
print("k=3, n=15, 균형. 모든 분포를 분산 1 로 표준화. 이동폭 0.9")
print(f"{'분포':>20s} {'F 크기':>8s} {'KW 크기':>8s} {'F 검정력':>9s} {'KW 검정력':>9s}")
for lab, f in [
    ("정규", lambda n: rng.normal(0, 1, n)),
    ("t(5)", lambda n: rng.standard_t(5, n) / np.sqrt(5 / 3)),
    ("지수", lambda n: rng.exponential(1, n)),
    ("로그정규", lambda n: (np.exp(rng.normal(0, 1, n)) - np.exp(.5)) / 2.16),
    ("균등", lambda n: rng.uniform(-np.sqrt(3), np.sqrt(3), n)),
]:
    a = b = c = e = 0
    for _ in range(B):
        gs = [f(15) for _ in range(3)]
        a += stats.f_oneway(*gs).pvalue < 0.05
        b += stats.kruskal(*gs).pvalue < 0.05
        gs = [f(15), f(15), f(15) + 0.9]
        c += stats.f_oneway(*gs).pvalue < 0.05
        e += stats.kruskal(*gs).pvalue < 0.05
    print(f"{lab:>20s} {a / B:8.4f} {b / B:8.4f} {c / B:9.4f} {e / B:9.4f}")

k=3, n=15, 균형. 모든 분포를 분산 1 로 표준화. 이동폭 0.9
                  분포     F 크기    KW 크기     F 검정력    KW 검정력
                  정규   0.0492   0.0455    0.6955    0.6628
                t(5)   0.0443   0.0452    0.7088    0.7577
                  지수   0.0395   0.0425    0.7142    0.9103
                로그정규   0.0367   0.0455    0.7807    0.9935
                  균등   0.0490   0.0465    0.6773    0.6020

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

def perm_p(gs, rng, R=400):
    obs = stats.f_oneway(*gs).statistic
    z = np.concatenate(gs)
    ns = [len(g) for g in gs]
    c = 0
    for _ in range(R):
        p = rng.permutation(z)
        i, parts = 0, []
        for n in ns:
            parts.append(p[i:i + n])
            i += n
        c += stats.f_oneway(*parts).statistic >= obs
    return (c + 1) / (R + 1)

rng = np.random.default_rng(4545)
B = 2_000
print("k=3, n=12, 모든 집단이 같은 분포, 명목 0.05")
print(f"{'분포':>16s} {'F 크기':>8s} {'순열 크기':>9s} {'KW 크기':>8s}")
for lab, f in [("정규", lambda n: rng.normal(0, 1, n)),
               ("t(3) 매우 두꺼움", lambda n: rng.standard_t(3, n)),
               ("로그정규", lambda n: np.exp(rng.normal(0, 1, n))),
               ("이항(0/1)", lambda n: rng.binomial(1, 0.2, n).astype(float))]:
    a = b = c = 0
    for _ in range(B):
        gs = [f(12) for _ in range(3)]
        a += stats.f_oneway(*gs).pvalue < 0.05
        b += perm_p(gs, rng) < 0.05
        c += stats.kruskal(*gs).pvalue < 0.05
    print(f"{lab:>16s} {a / B:8.4f} {b / B:9.4f} {c / B:8.4f}")

k=3, n=12, 모든 집단이 같은 분포, 명목 0.05
              분포     F 크기     순열 크기    KW 크기
              정규   0.0440    0.0460   0.0420
     t(3) 매우 두꺼움   0.0500    0.0595   0.0570
            로그정규   0.0345    0.0510   0.0350
         이항(0/1)   0.0605    0.0420   0.0455

import numpy as np
from scipy import stats

rng = np.random.default_rng(4646)
B = 20_000
print("지수분포(평균 1, 치우침 2), 명목 95%")
print(f"{'n':>5s} {'평균의 t 구간':>12s} {'예측구간':>10s}")
for n in [10, 30, 100, 500]:
    a = b = 0
    for _ in range(B):
        x = rng.exponential(1, n)
        m, s = x.mean(), x.std(ddof=1)
        t = stats.t.ppf(0.975, n - 1)
        a += abs(m - 1.0) < t * s / np.sqrt(n)          # 평균의 구간
        b += abs(rng.exponential(1) - m) < t * s * np.sqrt(1 + 1 / n)  # 예측구간
    print(f"{n:5d} {a / B:12.4f} {b / B:10.4f}")

xs = rng.exponential(1, (200_000, 10))
print("\n점추정의 불편성 (지수, n=10, 20만 회)")
print(f"  표본평균의 평균 = {xs.mean(1).mean():.6f}  (참값 1)")
print(f"  표본분산의 평균 = {xs.var(1, ddof=1).mean():.6f}  (참값 1)")

지수분포(평균 1, 치우침 2), 명목 95%
    n     평균의 t 구간       예측구간
   10       0.9041     0.9312
   30       0.9260     0.9415
  100       0.9434     0.9450
  500       0.9497     0.9446

점추정의 불편성 (지수, n=10, 20만 회)
  표본평균의 평균 = 1.000069  (참값 1)
  표본분산의 평균 = 1.004299  (참값 1)

잔차의 치우침 = 0.31, 초과첨도 = 0.48
Q-Q 그림에서 뚜렷한 체계적 이탈 없음 (양 꼬리에 각 1점)
샤피로-윌크 W = 0.978, p = 0.14 (n = 60)
→ 정규성 가정에 큰 문제 없음. 꼬리의 두 점은 관측 오류가 아님을 확인.