분산분석 진단¶

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

rng = np.random.default_rng(42)
n = 20
clean = np.concatenate([
    rng.normal(10.0, 1.0, n), rng.normal(10.8, 1.3, n), rng.normal(12.0, 1.6, n)])
dirty = clean.copy()
dirty[-1] = 20.0                                  # 이 쪽이 쓰는 자료
N, k = 3 * n, 3

# 제곱합 분해
grand = dirty.mean()
mi = [dirty[i * n:(i + 1) * n].mean() for i in range(k)]
SST = ((dirty - grand) ** 2).sum()
SSB = n * sum((a - grand) ** 2 for a in mi)
SSE = sum(((dirty[i * n:(i + 1) * n] - mi[i]) ** 2).sum() for i in range(k))
print(f"SST = {SST:.6f}")
print(f"SSB = {SSB:.6f}   (자유도 {k - 1})")
print(f"SSE = {SSE:.6f}   (자유도 {N - k})")
print(f"SSB + SSE = {SSB + SSE:.6f}   차 = {SST - SSB - SSE:.2e}")
print(f"eta^2 = SSB/SST = {SSB / SST:.6f}")
print(f"F = (SSB/{k - 1})/(SSE/{N - k}) = {(SSB / (k - 1)) / (SSE / (N - k)):.6f}")

# 원자료에 검정하면 안 되는 까닭
print(f"\n{'자료':>12}{'원자료 W':>11}{'원자료 p':>12}{'잔차 W':>10}{'잔차 p':>12}")
for lab, arr in [("이상점 있음", dirty), ("이상점 없음", clean)]:
    dd = pd.DataFrame({"group": np.repeat(["A", "B", "C"], n), "response": arr})
    res = ols("response ~ C(group)", data=dd).fit().resid
    wr, pr = shapiro(arr)
    we, pe = shapiro(res)
    print(f"{lab:>12}{wr:>11.4f}{pr:>12.4g}{we:>10.4f}{pe:>12.4g}")

SST = 182.674007
SSB = 66.024886   (자유도 2)
SSE = 116.649122   (자유도 57)
SSB + SSE = 182.674007   차 = -1.42e-14
eta^2 = SSB/SST = 0.361436
F = (SSB/2)/(SSE/57) = 16.131362

          자료      원자료 W       원자료 p      잔차 W        잔차 p
      이상점 있음     0.8449   2.165e-06    0.8116   2.658e-07
      이상점 없음     0.9834      0.5886    0.9870      0.7711

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

# 이 페이지의 진단은 모두 아래 모형 하나를 놓고 수행한다.
# 집단마다 표준편차를 1.0, 1.3, 1.6으로 다르게 주었고,
# 집단 C에 이상점을 하나 심어 두었다.
rng = np.random.default_rng(42)
n = 20
response = np.concatenate([
    rng.normal(10.0, 1.0, n),
    rng.normal(10.8, 1.3, n),
    rng.normal(12.0, 1.6, n),
])
response[-1] = 20.0                     # 마지막 관측값을 이상점으로 만든다
data = pd.DataFrame({
    "group": np.repeat(["A", "B", "C"], n),
    "response": response,
})
group1 = data.loc[data["group"] == "A", "response"]
group2 = data.loc[data["group"] == "B", "response"]
group3 = data.loc[data["group"] == "C", "response"]

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}")

       count    mean    std
group                      
A         20   9.967  0.870
B         20  10.942  1.034
C         20  12.513  2.077

F = 16.1314, p = 0.0000

from scipy.stats import shapiro
import statsmodels.api as sm

# 정규성은 자료가 아니라 잔차에 요구되는 가정이다. 집단마다 평균이 다르므로
# 자료 전체를 한 번에 검정하면 안 된다.
resid = model.resid
stat, p_value = shapiro(resid)
print(f"Shapiro-Wilk: W = {stat:.4f}, p = {p_value:.4f}")
sm.qqplot(resid, line='s')

Shapiro-Wilk: W = 0.8116, p = 0.0000

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

rng = np.random.default_rng(42)
n = 20
clean = np.concatenate([
    rng.normal(10.0, 1.0, n), rng.normal(10.8, 1.3, n), rng.normal(12.0, 1.6, n)])
response = clean.copy()
response[-1] = 20.0
data = pd.DataFrame({"group": np.repeat(["A", "B", "C"], n), "response": response})
model = ols("response ~ C(group)", data=data).fit()
e = model.resid.values
N = len(e)

pp = sm.ProbPlot(e)
tq, sq = pp.theoretical_quantiles, pp.sample_quantiles
s = e.std(ddof=1)
print(f"가장 큰 점의 이론 분위수 Phi^-1(60/61) = {tq[-1]:.4f}")
print(f"line='s' 가 그 자리에 기대하는 높이 = s * q = {s:.4f} * {tq[-1]:.4f} = {s * tq[-1]:.4f}")
print(f"실제 높이 = {sq[-1]:.4f}   ->  {sq[-1] / (s * tq[-1]):.2f} 배")
print(f"둘째로 큰 점: 기대 {s * tq[-2]:.4f}, 실제 {sq[-2]:.4f}  ({sq[-2] / (s * tq[-2]):.2f} 배)")

print(f"\nShapiro W = {stats.shapiro(e).statistic:.6f}")
print(f"정렬 잔차와 정규점수의 상관계수 제곱 = {np.corrcoef(sq, tq)[0, 1] ** 2:.6f}")

print(f"\n{'자료':>16}{'W':>10}{'p':>12}")
for lab, arr in [("전체 60 개", e), ("이상점 뺀 59 개", np.sort(e)[:-1])]:
    w, p = stats.shapiro(arr)
    print(f"{lab:>16}{w:>10.4f}{p:>12.4g}")
print(f"\n왜도 = {stats.skew(e):.4f},  이상점을 빼면 {stats.skew(np.sort(e)[:-1]):.4f}")
print(f"첨도(초과) = {stats.kurtosis(e):.4f},  이상점을 빼면 {stats.kurtosis(np.sort(e)[:-1]):.4f}")

가장 큰 점의 이론 분위수 Phi^-1(60/61) = 2.1347
line='s' 가 그 자리에 기대하는 높이 = s * q = 1.4061 * 2.1347 = 3.0016
실제 높이 = 7.4865   ->  2.49 배
둘째로 큰 점: 기대 2.5891, 실제 2.6418  (1.02 배)

Shapiro W = 0.811595
정렬 잔차와 정규점수의 상관계수 제곱 = 0.778751

              자료         W           p
         전체 60 개    0.8116   2.658e-07
      이상점 뺀 59 개    0.9906       0.931

왜도 = 2.4350,  이상점을 빼면 -0.0123
첨도(초과) = 11.7232,  이상점을 빼면 0.1740

from scipy.stats import levene, bartlett

# 두 검정을 함께 돌려 결론이 갈리는지 본다. 갈린다면 정규성이 의심스럽다는
# 뜻이므로, 정규성을 덜 타는 Levene 쪽을 믿는다.
groups = [data[data['group'] == g]['response'].values for g in data['group'].unique()]
stat_lev, p_lev = levene(*groups)
stat_bart, p_bart = bartlett(*groups)
print(f"Levene:   W = {stat_lev:.4f}, p = {p_lev:.4f}")
print(f"Bartlett: chi2 = {stat_bart:.4f}, p = {p_bart:.4f}")

Levene:   W = 1.1666, p = 0.3188
Bartlett: chi2 = 16.6837, p = 0.0002

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
n, k = 20, 3
clean = np.concatenate([
    rng.normal(10.0, 1.0, n), rng.normal(10.8, 1.3, n), rng.normal(12.0, 1.6, n)])
dirty = clean.copy()
dirty[-1] = 20.0
N = n * k
gs = [dirty[i * n:(i + 1) * n] for i in range(k)]

# Bartlett 를 공식으로
s2 = np.array([g.var(ddof=1) for g in gs])
sp2 = ((n - 1) * s2).sum() / (N - k)
C_bart = 1 + (1 / (3 * (k - 1))) * (k / (n - 1) - 1 / (N - k))
chi2 = ((N - k) * np.log(sp2) - ((n - 1) * np.log(s2)).sum()) / C_bart
print(f"s_i^2 = {np.round(s2, 4)},  합동 s_p^2 = {sp2:.6f}")
print(f"보정계수 C = {C_bart:.6f}")
print(f"공식  chi2 = {chi2:.6f},  p = {stats.chi2.sf(chi2, k - 1):.6f}")
print(f"scipy chi2 = {stats.bartlett(*gs).statistic:.6f},  p = {stats.bartlett(*gs).pvalue:.6f}")

# Levene 는 |y - 중앙값| 에 대한 일원배치 F
Z = [np.abs(g - np.median(g)) for g in gs]
print(f"\nZbar = {[round(z.mean(), 4) for z in Z]}")
print(f"공식  W = {stats.f_oneway(*Z).statistic:.6f},  p = {stats.f_oneway(*Z).pvalue:.6f}")
print(f"scipy W = {stats.levene(*gs).statistic:.6f},  p = {stats.levene(*gs).pvalue:.6f}")

# 이상점을 빼면
print(f"\n{'자료':>12}{'s_A':>8}{'s_B':>8}{'s_C':>8}{'Bartlett p':>12}{'Levene p':>11}")
for lab, arr in [("이상점 있음", dirty), ("이상점 없음", clean)]:
    g = [arr[i * n:(i + 1) * n] for i in range(k)]
    sd = [x.std(ddof=1) for x in g]
    print(f"{lab:>12}{sd[0]:>8.4f}{sd[1]:>8.4f}{sd[2]:>8.4f}"
          f"{stats.bartlett(*g).pvalue:>12.4f}{stats.levene(*g).pvalue:>11.4f}")

s_i^2 = [0.7572 1.0682 4.314 ],  합동 s_p^2 = 2.046476
보정계수 C = 1.023392
공식  chi2 = 16.683675,  p = 0.000238
scipy chi2 = 16.683675,  p = 0.000238

Zbar = [0.6952, 0.818, 1.1998]
공식  W = 1.166587,  p = 0.318752
scipy W = 1.166587,  p = 0.318752

          자료     s_A     s_B     s_C  Bartlett p   Levene p
      이상점 있음  0.8702  1.0335  2.0770      0.0002     0.3188
      이상점 없음  0.8702  1.0335  1.1450      0.4996     0.6662

from statsmodels.stats.stattools import durbin_watson

# 통계량은 0 에서 4 사이이고 2 가 무상관에 해당한다. 2 보다 뚜렷이 작으면
# 양의 자기상관, 크면 음의 자기상관이다.
dw = durbin_watson(model.resid)
print(f"Durbin-Watson: {dw:.4f}")

Durbin-Watson: 1.7586

import numpy as np
import pandas as pd
from statsmodels.formula.api import ols
from statsmodels.stats.stattools import durbin_watson

rng = np.random.default_rng(42)
n = 20
response = np.concatenate([
    rng.normal(10.0, 1.0, n), rng.normal(10.8, 1.3, n), rng.normal(12.0, 1.6, n)])
response[-1] = 20.0
data = pd.DataFrame({"group": np.repeat(["A", "B", "C"], n), "response": response})
e = ols("response ~ C(group)", data=data).fit().resid.values
S = (e ** 2).sum()

rho = (e[1:] * e[:-1]).sum() / S
edge = (e[0] ** 2 + e[-1] ** 2) / S
print(f"rho_1 = {rho:.6f}")
print(f"2(1 - rho_1)                = {2 * (1 - rho):.6f}")
print(f"끝점 보정 (e_1^2+e_N^2)/sum = {edge:.6f}")
print(f"2(1-rho) - 끝점 보정        = {2 * (1 - rho) - edge:.6f}")
print(f"durbin_watson               = {durbin_watson(e):.6f}")
print(f"\n끝점 보정 중 마지막 관측 혼자의 몫: e_N^2/sum = {e[-1] ** 2 / S:.6f}")

# 행 순서를 바꾸면 통계량이 어떻게 되는가
rg = np.random.default_rng(7)
dws = np.array([durbin_watson(rg.permutation(e)) for _ in range(20_000)])
print(f"\n행을 무작위로 섞어 2만 번 다시 잰 DW")
print(f"  평균 {dws.mean():.4f},  표준편차 {dws.std(ddof=1):.4f}")
print(f"  2.5% ~ 97.5% 분위: [{np.quantile(dws, 0.025):.4f}, {np.quantile(dws, 0.975):.4f}]")
print(f"  관측값 {durbin_watson(e):.4f} 의 백분위 = {np.mean(dws < durbin_watson(e)):.3f}")
print(f"  경험칙 (1.5, 2.5) 안에 들어가는 비율 = {np.mean((dws > 1.5) & (dws < 2.5)):.3f}")

rho_1 = -0.120015
2(1 - rho_1)                = 2.240031
끝점 보정 (e_1^2+e_N^2)/sum = 0.481463
2(1-rho) - 끝점 보정        = 1.758567
durbin_watson               = 1.758567

끝점 보정 중 마지막 관측 혼자의 몫: e_N^2/sum = 0.480486

행을 무작위로 섞어 2만 번 다시 잰 DW
  평균 1.9959,  표준편차 0.2409
  2.5% ~ 97.5% 분위: [1.4866, 2.4449]
  관측값 1.7586 의 백분위 = 0.154
  경험칙 (1.5, 2.5) 안에 들어가는 비율 = 0.959

# Cook 의 거리는 그 관측값 하나를 뺐을 때 적합값 전체가 얼마나 움직이는지를
# 잰다. 크다는 것은 결론이 그 한 점에 기대고 있다는 뜻이다.
influence = model.get_influence()
cooks_d = influence.cooks_distance[0]

# 4/n 은 널리 쓰이는 어림 기준일 뿐 검정이 아니다. 넘는 점은 지울 대상이
# 아니라 들여다볼 대상이다.
threshold = 4 / len(cooks_d)
flagged = np.where(cooks_d > threshold)[0]
print(f"threshold = {threshold:.4f}")
print(f"flagged observations = {flagged}")
print(f"max Cook's D = {cooks_d.max():.4f} (obs {cooks_d.argmax()})")

threshold = 0.0667
flagged observations = [52 59]
max Cook's D = 0.5058 (obs 59)

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

rng = np.random.default_rng(42)
n, k = 20, 3
response = np.concatenate([
    rng.normal(10.0, 1.0, n), rng.normal(10.8, 1.3, n), rng.normal(12.0, 1.6, n)])
response[-1] = 20.0
data = pd.DataFrame({"group": np.repeat(["A", "B", "C"], n), "response": response})
model = ols("response ~ C(group)", data=data).fit()
N = 3 * n

inf = model.get_influence()
D = inf.cooks_distance[0]
r = inf.resid_studentized_internal
h = inf.hat_matrix_diag

coef = (1 / k) * h[0] / (1 - h[0])
print(f"h_ii 는 모두 {h[0]:.4f},  (1/p)*h/(1-h) = {coef:.8f} = 1/{1 / coef:.4f}")
print(f"D 와 r^2/(N-k) 의 최대 차이 = {np.abs(D - r ** 2 / (N - k)).max():.2e}")

print(f"\n{'관측':>6}{'r':>10}{'r^2/57':>10}{'Cook D':>10}")
for i in np.argsort(-D)[:4]:
    print(f"{i:>6}{r[i]:>10.4f}{r[i] ** 2 / (N - k):>10.6f}{D[i]:>10.6f}")

thr = 4 / N
print(f"\n문턱 4/N = {thr:.6f}")
print(f"이에 대응하는 |r| 문턱 = sqrt(4(N-k)/N) = {np.sqrt(4 * (N - k) / N):.6f}")
print(f"D > 4/N 인 관측 = {np.where(D > thr)[0]}")
print(f"|r| > 2  인 관측 = {np.where(np.abs(r) > 2)[0]}")
print(f"두 집합이 같은가: {set(np.where(D > thr)[0]) == set(np.where(np.abs(r) > 2)[0])}")

# D 순위와 |r| 순위가 같은가
print(f"D 의 순위와 |r| 의 순위가 완전히 같은가: "
      f"{np.array_equal(np.argsort(D), np.argsort(np.abs(r)))}")

h_ii 는 모두 0.0500,  (1/p)*h/(1-h) = 0.01754386 = 1/57.0000
D 와 r^2/(N-k) 의 최대 차이 = 2.22e-16

    관측         r    r^2/57    Cook D
    59    5.3693  0.505775  0.505775
    52   -2.0403  0.073035  0.073035
    30    1.8947  0.062981  0.062981
     4   -1.3756  0.033200  0.033200

문턱 4/N = 0.066667
이에 대응하는 |r| 문턱 = sqrt(4(N-k)/N) = 1.949359
D > 4/N 인 관측 = [52 59]
|r| > 2  인 관측 = [52 59]
두 집합이 같은가: True
D 의 순위와 |r| 의 순위가 완전히 같은가: True

import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.formula.api import ols

def run_full_diagnostics(data, response_col, group_col):
    """앞의 진단 넷을 한 번에 돌려 2x2 격자로 보여 준다.

    실제 분석에서는 이 네 그림을 늘 함께 본다. 하나만 보고 판단하면
    다른 쪽에서 드러날 문제를 놓치기 쉽다.
    """
    formula = f'{response_col} ~ {group_col}'
    model = ols(formula, data=data).fit()
    anova_table = sm.stats.anova_lm(model, typ=2)
    print(anova_table)

    fig, axes = plt.subplots(2, 2, figsize=(12, 8))
    # 왼쪽 위: Q-Q 그림 — 정규성
    sm.qqplot(model.resid, line='s', ax=axes[0, 0])
    # 오른쪽 위: 잔차 히스토그램 — 치우침과 봉우리
    axes[0, 1].hist(model.resid, bins=15, density=True, alpha=0.7, edgecolor='black')
    # 왼쪽 아래: 잔차 대 적합값 — 등분산성과 남은 구조
    axes[1, 0].scatter(model.fittedvalues, model.resid, alpha=0.6)
    axes[1, 0].axhline(y=0, color='r', linestyle='--')
    # 오른쪽 아래: Cook 의 거리 — 영향점
    cooks_d = model.get_influence().cooks_distance[0]
    axes[1, 1].stem(range(len(cooks_d)), cooks_d, markerfmt=",")
    axes[1, 1].axhline(y=4 / len(cooks_d), color='r', linestyle='--')
    plt.tight_layout()
    plt.show()

run_full_diagnostics(data, "response", "C(group)")

              sum_sq    df          F    PR(>F)
C(group)   66.024886   2.0  16.131362  0.000003
Residual  116.649122  57.0        NaN       NaN

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

rng = np.random.default_rng(42)
n, k = 20, 3
response = np.concatenate([
    rng.normal(10.0, 1.0, n), rng.normal(10.8, 1.3, n), rng.normal(12.0, 1.6, n)])
response[-1] = 20.0
data = pd.DataFrame({"group": np.repeat(["A", "B", "C"], n), "response": response})
model = ols("response ~ C(group)", data=data).fit()
e = model.resid.values
inf = model.get_influence()
r = inf.resid_studentized_internal
D = inf.cooks_distance[0]
N = 3 * n

# 왼쪽 위 Q-Q: 보기 2 에서 잰 값
pp = sm.ProbPlot(e)
s = e.std(ddof=1)
print("[왼쪽 위] Q-Q")
print(f"  가장 큰 점: 기대 {s * pp.theoretical_quantiles[-1]:.4f}, 실제 {pp.sample_quantiles[-1]:.4f}")
print(f"  나머지 59 개의 선으로부터의 최대 거리 = "
      f"{np.abs(pp.sample_quantiles[:-1] - s * pp.theoretical_quantiles[:-1]).max():.4f}")

# 오른쪽 위 히스토그램
print("\n[오른쪽 위] 히스토그램 bins=15")
cnt, edges = np.histogram(e, bins=15)
w = edges[1] - edges[0]
print(f"  구간 폭 = {w:.4f},  칸별 도수 = {cnt}")
print(f"  마지막 칸의 도수 = {cnt[-1]},  그 앞의 연속된 빈 칸 수 = {int((cnt[:-1][::-1] != 0).argmax())}")
print(f"  밀도 한 칸 높이 = 1/(N*w) = {1 / (N * w):.4f}")

# 왼쪽 아래 잔차 대 적합값
print("\n[왼쪽 아래] 잔차 대 적합값")
for i, g in enumerate("ABC"):
    ei = e[i * n:(i + 1) * n]
    print(f"  띠 {g}: 위치 {model.fittedvalues.values[i * n]:.4f}, "
          f"범위 [{ei.min():.4f}, {ei.max():.4f}], 폭 {ei.max() - ei.min():.4f}")

# 오른쪽 아래 Cook
print("\n[오른쪽 아래] Cook")
top = np.sort(D)[::-1][:3]
print(f"  가장 큰 셋 = {np.round(top, 4)},  최대/둘째 = {top[0] / top[1]:.2f}")
print(f"  문턱 4/N = {4 / N:.4f},  넘는 관측 = {np.where(D > 4 / N)[0]}")

# 네 그림이 서로 얼마나 겹치는가
print("\n네 그림의 재료가 서로 얼마나 같은가")
print(f"  Cook D 와 |표준화 잔차| 의 순위상관 = "
      f"{np.corrcoef(np.argsort(np.argsort(D)), np.argsort(np.argsort(np.abs(r))))[0, 1]:.4f}")
print(f"  Q-Q 와 히스토그램은 같은 60 개 잔차의 주변분포를 두 방식으로 그린 것")
print(f"  잔차 그림의 세로축 x {1 / (np.sqrt(model.mse_resid) * np.sqrt(1 - 1 / n)):.4f} = 표준화 잔차")

[왼쪽 위] Q-Q
  가장 큰 점: 기대 3.0016, 실제 7.4865
  나머지 59 개의 선으로부터의 최대 거리 = 0.7161

[오른쪽 위] 히스토그램 bins=15
  구간 폭 = 0.6888,  칸별 도수 = [ 1  4 13 12 15 10  3  1  0  0  0  0  0  0  1]
  마지막 칸의 도수 = 1,  그 앞의 연속된 빈 칸 수 = 6
  밀도 한 칸 높이 = 1/(N*w) = 0.0242

[왼쪽 아래] 잔차 대 적합값
  띠 A: 위치 9.9671, 범위 [-1.9181, 1.1602], 폭 3.0783
  띠 B: 위치 10.9423, 범위 [-1.2345, 2.6418], 폭 3.8763
  띠 C: 위치 12.5135, 범위 [-2.8449, 7.4865], 폭 10.3314

[오른쪽 아래] Cook
  가장 큰 셋 = [0.5058 0.073  0.063 ],  최대/둘째 = 6.93
  문턱 4/N = 0.0667,  넘는 관측 = [52 59]

네 그림의 재료가 서로 얼마나 같은가
  Cook D 와 |표준화 잔차| 의 순위상관 = 1.0000
  Q-Q 와 히스토그램은 같은 60 개 잔차의 주변분포를 두 방식으로 그린 것
  잔차 그림의 세로축 x 0.7172 = 표준화 잔차

import warnings
warnings.filterwarnings("ignore")

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

rng = np.random.default_rng(8001)
B = 3_000
print("모든 관측이 같은 정규분포 — 이상점이 하나도 없는 자료")
print(f"{'k':>3s} {'집단당 n':>8s} {'N':>5s} {'4/N':>8s} "
      f"{'D>4/N 인 관측 비율':>16s} {'적어도 하나':>10s} {'D>1 비율':>9s}")
for k, n in [(3, 10), (3, 20), (4, 15), (5, 30), (3, 50)]:
    N = k * n
    thr = 4 / N
    frac, any_, big = [], 0, 0
    for _ in range(B):
        rows = [pd.DataFrame({"y": rng.normal(0, 1, n), "g": f"G{g}"})
                for g in range(k)]
        df = pd.concat(rows, ignore_index=True)
        D = ols("y ~ C(g)", data=df).fit().get_influence().cooks_distance[0]
        frac.append((D > thr).mean())
        any_ += (D > thr).any()
        big += (D > 1).any()
    print(f"{k:3d} {n:8d} {N:5d} {thr:8.4f} {np.mean(frac):16.4f} "
          f"{any_ / B:10.4f} {big / B:9.4f}")

모든 관측이 같은 정규분포 — 이상점이 하나도 없는 자료
  k    집단당 n     N      4/N    D>4/N 인 관측 비율     적어도 하나    D>1 비율
  3       10    30   0.1333           0.0560     0.9303    0.0000
  3       20    60   0.0667           0.0502     0.9933    0.0000
  4       15    60   0.0667           0.0526     0.9937    0.0000
  5       30   150   0.0267           0.0487     1.0000    0.0000
  3       50   150   0.0267           0.0473     1.0000    0.0000

import warnings
warnings.filterwarnings("ignore")

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

rng = np.random.default_rng(8002)
rows = [pd.DataFrame({"y": rng.normal(mu, 1.2, n), "g": f"G{g}"})
        for g, (n, mu) in enumerate([(8, 10.0), (12, 12.0), (6, 11.0)])]
df = pd.concat(rows, ignore_index=True)
df.loc[len(df) - 1, "y"] = 18.0          # 작은 집단에 이상점 하나

fit = ols("y ~ C(g)", data=df).fit()
inf = fit.get_influence()
t = pd.DataFrame({
    "g": df.g, "y": df.y.round(3),
    "잔차": fit.resid.round(3),
    "h": inf.hat_matrix_diag.round(4),
    "표준화": inf.resid_studentized_internal.round(3),
    "스튜던트화(삭제)": inf.resid_studentized_external.round(3),
    "CookD": inf.cooks_distance[0].round(4),
    "DFFITS": inf.dffits[0].round(3)})
print(t.tail(8).to_string(index=False))

N, p = len(df), 3
print(f"\nN = {N}, p = {p}")
print(f"  h 의 값: {sorted(set(np.round(inf.hat_matrix_diag, 4)))}")
print(f"  집단 크기: {df.g.value_counts().sort_index().to_dict()}")
print(f"  Σh = {inf.hat_matrix_diag.sum():.4f}  (이론 p = {p})")
print(f"\n문턱: 4/N = {4 / N:.4f},  DFFITS 2√(p/N) = {2 * np.sqrt(p / N):.4f}, "
      f"h 의 2p/N = {2 * p / N:.4f}")

i = len(df) - 1
r = inf.resid_studentized_internal[i]
h = inf.hat_matrix_diag[i]
te = inf.resid_studentized_external[i]
print(f"\n관계 확인 (마지막 관측):")
print(f"  Cook D = r²/p · h/(1-h) = {r**2 / p * h / (1 - h):.6f}   "
      f"실제 {inf.cooks_distance[0][i]:.6f}")
print(f"  DFFITS = t·√(h/(1-h))  = {te * np.sqrt(h / (1 - h)):.6f}   "
      f"실제 {inf.dffits[0][i]:.6f}")

 g      y     잔차      h    표준화  스튜던트화(삭제)  CookD  DFFITS
G1 12.268 -0.238 0.0833 -0.169     -0.165 0.0009  -0.050
G1 12.343 -0.164 0.0833 -0.116     -0.113 0.0004  -0.034
G2 11.720 -1.141 0.1667 -0.847     -0.841 0.0478  -0.376
G2 10.570 -2.291 0.1667 -1.701     -1.779 0.1929  -0.796
G2 12.256 -0.605 0.1667 -0.449     -0.441 0.0134  -0.197
G2 13.587  0.726 0.1667  0.539      0.531 0.0194   0.237
G2 11.031 -1.830 0.1667 -1.359     -1.385 0.1230  -0.620
G2 18.000  5.140 0.1667  3.816      6.162 0.9709   2.756

N = 26, p = 3
  h 의 값: [0.0833, 0.125, 0.1667]
  집단 크기: {'G0': 8, 'G1': 12, 'G2': 6}
  Σh = 3.0000  (이론 p = 3)

문턱: 4/N = 0.1538,  DFFITS 2√(p/N) = 0.6794, h 의 2p/N = 0.2308

관계 확인 (마지막 관측):
  Cook D = r²/p · h/(1-h) = 0.970882   실제 0.970882
  DFFITS = t·√(h/(1-h))  = 2.755922   실제 2.755922

import warnings
warnings.filterwarnings("ignore")

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

rng = np.random.default_rng(8003)
B = 3_000
print("자료에 이상점이 전혀 없을 때, 적어도 하나가 표시될 확률")
print(f"{'N':>5s} {'|t*|>2':>9s} {'|t*|>3':>9s} {'본페로니 보정':>12s} {'보정 임계값':>11s}")
for k, n in [(3, 10), (3, 20), (4, 25), (5, 40)]:
    N = k * n
    tb = stats.t.ppf(1 - 0.05 / (2 * N), N - k - 1)
    a = b = c = 0
    for _ in range(B):
        rows = [pd.DataFrame({"y": rng.normal(0, 1, n), "g": f"G{g}"})
                for g in range(k)]
        df = pd.concat(rows, ignore_index=True)
        te = (ols("y ~ C(g)", data=df).fit()
              .get_influence().resid_studentized_external)
        a += (np.abs(te) > 2).any()
        b += (np.abs(te) > 3).any()
        c += (np.abs(te) > tb).any()
    print(f"{N:5d} {a / B:9.4f} {b / B:9.4f} {c / B:12.4f} {tb:11.4f}")

자료에 이상점이 전혀 없을 때, 적어도 하나가 표시될 확률
    N    |t*|>2    |t*|>3      본페로니 보정      보정 임계값
   30    0.9340    0.1733       0.0533      3.5069
   60    0.9923    0.2417       0.0513      3.5322
  100    1.0000    0.3057       0.0453      3.6047
  200    1.0000    0.4857       0.0533      3.7314

p_i = 2 · P(t_{N-p-1} > |t_i|)          각 관측의 원 p-값
p_i^adj = min(1, N · p_i)               본페로니 보정
→ 가장 작은 p_i^adj 만 보고

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

def welch_p(gs):
    n = np.array([len(g) for g in gs], float)
    m = np.array([g.mean() for g in gs])
    v = np.array([g.var(ddof=1) for g in gs])
    k = len(n)
    w = n / v
    W = w.sum()
    mt = (w * m).sum() / W
    lam = ((1 - w / W)**2 / (n - 1)).sum()
    F = ((w * (m - mt)**2).sum() / (k - 1)
         / (1 + 2 * (k - 2) / (k**2 - 1) * lam))
    return stats.f.sf(F, k - 1, (3 / (k**2 - 1) * lam)**-1)

rng = np.random.default_rng(8004)
B = 10_000
NS = [50, 12, 45]
print("n = (50, 12, 45), 모든 평균이 같음, 명목 0.05")
print(f"{'시나리오':>28s} {'표준 F':>8s} {'웰치':>8s}")
for lab, sds in [("역페어링: 작은 집단에 큰 분산 σ=(1,3,1)", [1, 3, 1]),
                 ("정페어링: 작은 집단에 작은 분산 σ=(3,1,3)", [3, 1, 3]),
                 ("큰 집단끼리 다름 σ=(1,1,3)", [1, 1, 3])]:
    a = b = 0
    for _ in range(B):
        gs = [rng.normal(0, s, n) for n, s in zip(NS, sds)]
        a += stats.f_oneway(*gs).pvalue < 0.05
        b += welch_p(gs) < 0.05
    print(f"{lab:>28s} {a / B:8.4f} {b / B:8.4f}")

n = (50, 12, 45), 모든 평균이 같음, 명목 0.05
                        시나리오     표준 F       웰치
 역페어링: 작은 집단에 큰 분산 σ=(1,3,1)   0.2656   0.0541
정페어링: 작은 집단에 작은 분산 σ=(3,1,3)   0.0259   0.0515
         큰 집단끼리 다름 σ=(1,1,3)   0.0363   0.0500

[설계]
  □ 관측이 독립인가 (반복측정·군집·시계열?)
  □ 집단별 n 이 균형에 가까운가

[잔차]
  □ 잔차 대 적합값 그림 — 깔때기 모양?
  □ Q-Q 그림 — 꼬리, 치우침, 계단?
  □ 관측 순서에 대한 잔차 — 추세·주기?

[이상점·영향점]
  □ 스튜던트화(삭제) 잔차의 최댓값 — 본페로니 보정 임계값과 비교
  □ 쿡 거리 — 크기순 정렬, 뚜렷하게 튀는 점?
  □ 표시된 점을 빼고 다시 적합 — 결론이 바뀌는가?

[검정]
  □ 브라운-포사이드 (참고용)
  □ 정규성 검정 (n 이 크면 참고만)

1. 기록 오류인가?  ──→ 고치거나 결측 처리
2. 다른 모집단에서 왔는가?  ──→ 제외하고 그 사실을 보고
3. 그냥 극단값인가?  ──→ 남긴다
     ↓
4. 있을 때와 없을 때의 결과를 모두 계산
     ↓
5. 결론이 바뀌면 둘 다 보고, 바뀌지 않으면 그 사실을 보고