카이제곱 적합도 검정 (수동)¶

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

observed_counts = np.array([4, 13, 7])
# H0가 균등분포이므로 기대도수는 모두 n/k, 즉 관측도수의 평균과 같다.
expected_counts = np.ones(3) * observed_counts.mean()
degrees_of_freedom = observed_counts.shape[0] - 1   # 총합이 고정이라 하나를 잃는다

# 분모가 관측도수가 아니라 **기대도수**인 것에 주의하라.
# H0 아래에서 각 칸 도수의 분산이 근사적으로 E_i이기 때문이다(연습문제 4).
chi_square_statistic = np.sum(
    (observed_counts - expected_counts) ** 2 / expected_counts
)
# 카이제곱 적합도 검정은 언제나 우측검정이므로 sf(위쪽 꼬리)를 쓴다.
p_value = stats.chi2(degrees_of_freedom).sf(chi_square_statistic)

print(f"Chi-square Statistic = {chi_square_statistic:.4f}")
print(f"p-value = {p_value:.4f}")

Chi-square Statistic = 5.2500
p-value = 0.0724

fig, ax = plt.subplots(figsize=(12, 4))

# 왼쪽 구간 — 기각하지 않는 쪽
x_left = np.linspace(0, chi_square_statistic, 100)
y_left = stats.chi2(degrees_of_freedom).pdf(x_left)
ax.plot(x_left, y_left, linewidth=3)
x_fill_left = np.concatenate([[0], x_left, [chi_square_statistic], [0]])
y_fill_left = np.concatenate([[0], y_left, [0], [0]])
ax.fill(x_fill_left, y_fill_left, alpha=0.1)

# 오른쪽 꼬리 — 기각역
x_right = np.linspace(chi_square_statistic, 20, 100)
y_right = stats.chi2(degrees_of_freedom).pdf(x_right)
ax.plot(x_right, y_right, linewidth=3)
x_fill_right = np.concatenate(
    [[chi_square_statistic], x_right, [20], [chi_square_statistic]]
)
y_fill_right = np.concatenate([[0], y_right, [0], [0]])
ax.fill(x_fill_right, y_fill_right, alpha=0.1)

# p-값을 표시한다
ax.annotate(
    f"p-value = {p_value:.02%}",
    xy=((12.5 + 15.0) / 2, 0.01),
    xytext=(16.5, 0.10),
    fontsize=15,
    arrowprops=dict(width=0.2, headwidth=8),
)

ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)
ax.spines["bottom"].set_position("zero")
ax.spines["left"].set_position("zero")
plt.tight_layout()
plt.show()

import numpy as np
from scipy import stats

obs = np.array([4, 13, 7])
n = obs.sum()
exp = np.full(3, n / 3)
term = (obs - exp)**2 / exp

print("범주별 기여")
for i, (o, e, t) in enumerate(zip(obs, exp, term), 1):
    print(f"  범주{i}: O={o:3d}  E={e:6.3f}  (O-E)={o - e:+7.3f}  "
          f"기여 {t:7.4f}  ({t / term.sum() * 100:5.1f}%)")
print(f"  합계 χ² = {term.sum():.4f}\n")

resid = (obs - exp) / np.sqrt(exp)                     # 표준화 잔차
adj = (obs - exp) / np.sqrt(exp * (1 - 1 / 3))         # 조정 표준화 잔차
print(f"표준화 잔차      {np.round(resid, 4).tolist()}")
print(f"조정 표준화 잔차  {np.round(adj, 4).tolist()}")
print(f"본페로니 임계값(범주 3개): {stats.norm.ppf(1 - 0.025 / 3):.4f}")

범주별 기여
  범주1: O=  4  E= 8.000  (O-E)= -4.000  기여  2.0000  ( 38.1%)
  범주2: O= 13  E= 8.000  (O-E)= +5.000  기여  3.1250  ( 59.5%)
  범주3: O=  7  E= 8.000  (O-E)= -1.000  기여  0.1250  (  2.4%)
  합계 χ² = 5.2500

표준화 잔차      [-1.4142, 1.7678, -0.3536]
조정 표준화 잔차  [-1.7321, 2.1651, -0.433]
본페로니 임계값(범주 3개): 2.3940

import numpy as np
from math import lgamma
from scipy import stats

n, k = 24, 3
p = np.ones(k) / k
exp = n * p

# 가능한 모든 (a, b, c) 를 열거해 χ² 값별 확률을 모은다
table = {}
for a in range(n + 1):
    for b in range(n - a + 1):
        c = n - a - b
        t = round((((np.array([a, b, c]) - exp)**2) / exp).sum(), 9)
        log_p = (lgamma(n + 1) - lgamma(a + 1) - lgamma(b + 1)
                 - lgamma(c + 1) + n * np.log(1 / 3))
        table[t] = table.get(t, 0.0) + np.exp(log_p)

ts = np.array(sorted(table))
pr = np.array([table[t] for t in ts])
sf = np.cumsum(pr[::-1])[::-1]              # 오른쪽 꼬리 확률

print(f"가능한 도수 벡터 {(n + 1) * (n + 2) // 2}개, "
      f"서로 다른 χ² 값 {len(ts)}개")
print(f"통계량의 평균 {(ts * pr).sum():.4f}   (df = {k - 1})\n")

obs_stat = (((np.array([4, 13, 7]) - exp)**2) / exp).sum()
i = int(np.argmin(np.abs(ts - obs_stat)))
print(f"보기 1 [4, 13, 7]:  χ² = {obs_stat:.4f}")
print(f"   정확 p = {sf[i]:.4f}   근사 p = {stats.chi2.sf(obs_stat, k - 1):.4f}\n")

print("0.05 근처에서 도달 가능한 p 값")
mask = (sf > 0.01) & (sf < 0.15)
for t, s in zip(ts[mask], sf[mask]):
    print(f"   χ² = {t:7.4f}   정확 p = {s:.4f}   "
          f"근사 p = {stats.chi2.sf(t, k - 1):.4f}")

가능한 도수 벡터 325개, 서로 다른 χ² 값 44개
통계량의 평균 2.0000   (df = 2)

보기 1 [4, 13, 7]:  χ² = 5.2500
   정확 p = 0.0800   근사 p = 0.0724

0.05 근처에서 도달 가능한 p 값
   χ² =  4.0000   정확 p = 0.1481   근사 p = 0.1353
   χ² =  4.7500   정확 p = 0.1197   근사 p = 0.0930
   χ² =  5.2500   정확 p = 0.0800   근사 p = 0.0724
   χ² =  6.2500   정확 p = 0.0499   근사 p = 0.0439
   χ² =  6.7500   정확 p = 0.0412   근사 p = 0.0342
   χ² =  7.0000   정확 p = 0.0338   근사 p = 0.0302
   χ² =  7.7500   정확 p = 0.0213   근사 p = 0.0208
   χ² =  9.0000   정확 p = 0.0134   근사 p = 0.0111
   χ² =  9.2500   정확 p = 0.0115   근사 p = 0.0098

import numpy as np
from scipy import stats

print(f"{'k/n':>10s} {'χ²':>12s} {'z²':>12s} {'차이':>10s} "
      f"{'χ² 의 p':>9s} {'z 의 양측 p':>12s}")
for k, n in [(13, 24), (60, 100), (7, 20), (505, 1000)]:
    p0 = 0.5
    obs = np.array([k, n - k])
    exp = np.array([n * p0, n * (1 - p0)])
    chi2 = ((obs - exp)**2 / exp).sum()
    z = (k / n - p0) / np.sqrt(p0 * (1 - p0) / n)
    print(f"{f'{k}/{n}':>10s} {chi2:12.6f} {z**2:12.6f} "
          f"{abs(chi2 - z**2):10.2e} {stats.chi2.sf(chi2, 1):9.4f} "
          f"{2 * stats.norm.sf(abs(z)):12.4f}")

       k/n           χ²           z²         차이    χ² 의 p     z 의 양측 p
     13/24     0.166667     0.166667   3.05e-16    0.6831       0.6831
    60/100     4.000000     4.000000   1.78e-15    0.0455       0.0455
      7/20     1.800000     1.800000   6.66e-16    0.1797       0.1797
  505/1000     0.100000     0.100000   1.80e-16    0.7518       0.7518

import numpy as np
from scipy import stats

def gof_by_hand(obs, p_exp=None, n_estimated=0, alpha=0.05):
    """카이제곱 적합도 검정을 직접 계산한다.

    p_exp        : 기대비율. None 이면 균등분포.
    n_estimated  : 자료에서 추정한 모수의 개수 (자유도에서 뺀다).
    """
    obs = np.asarray(obs, float)
    n, k = obs.sum(), len(obs)
    p_exp = np.full(k, 1 / k) if p_exp is None else np.asarray(p_exp, float)
    exp = n * p_exp
    chi2 = ((obs - exp)**2 / exp).sum()
    df = k - 1 - n_estimated
    return {"chi2": chi2, "df": df, "p": stats.chi2.sf(chi2, df),
            "exp": exp, "E_min": exp.min(),
            "reject": stats.chi2.sf(chi2, df) < alpha}

cases = [
    ("보기 1 (균등)", [4, 13, 7], None, 0),
    ("불균등 기대비율", [4, 13, 7], [0.5, 0.3, 0.2], 0),
    ("주사위 120회", [15, 22, 18, 25, 20, 20], None, 0),
    ("포아송 적합(λ 추정)", [109, 65, 22, 4],
     [0.5434, 0.3314, 0.1011, 0.0241], 1),
]
for label, obs, p_exp, r in cases:
    out = gof_by_hand(obs, p_exp, r)
    chk = stats.chisquare(f_obs=obs,
                          f_exp=None if p_exp is None
                          else np.array(p_exp) * sum(obs),
                          ddof=r)
    print(f"{label}")
    print(f"  손계산   χ²={out['chi2']:8.4f}  df={out['df']}  "
          f"p={out['p']:.4f}  E_min={out['E_min']:.3f}")
    print(f"  scipy    χ²={chk.statistic:8.4f}  p={chk.pvalue:.4f}")
    print(f"  일치? {np.isclose(out['chi2'], chk.statistic)} / "
          f"{np.isclose(out['p'], chk.pvalue)}")

보기 1 (균등)
  손계산   χ²=  5.2500  df=2  p=0.0724  E_min=8.000
  scipy    χ²=  5.2500  p=0.0724
  일치? True / True
불균등 기대비율
  손계산   χ²= 11.0139  df=2  p=0.0041  E_min=4.800
  scipy    χ²= 11.0139  p=0.0041
  일치? True / True
주사위 120회
  손계산   χ²=  2.9000  df=5  p=0.7154  E_min=20.000
  scipy    χ²=  2.9000  p=0.7154
  일치? True / True
포아송 적합(λ 추정)
  손계산   χ²=  0.3219  df=2  p=0.8514  E_min=4.820
  scipy    χ²=  0.3219  p=0.8514
  일치? True / True

① n 과 k 를 센다
     n = 관측도수의 합,  k = 범주의 수
          ↓
② 기대비율 p_i 를 정한다
     ├─ 균등 → p_i = 1/k
     ├─ 이론값 → 문제에서 주어짐 (예: 3:1 분리비)
     └─ 모형 적합 → 추정한 모수로 계산 (추정 개수 r 을 기록)
          ↓
③ E_i = n p_i 를 계산한다
     └─ 검산: Σ E_i = n 이어야 한다
          ↓
④ E_min 을 확인한다
     └─ 5 미만이면 병합 규칙을 적용
          ↓
⑤ 항별로 (O_i - E_i)² / E_i 를 계산해 더한다
     └─ 검산: Σ(O_i - E_i) = 0 이어야 한다
          ↓
⑥ df = k - 1 - r
          ↓
⑦ p = P(χ²_df > 통계량)   ← 오른쪽 꼬리
          ↓
⑧ 잔차를 보고 어디가 어긋났는지 확인한다