자유도와 점근 이론¶

import numpy as np
from scipy import stats

rng = np.random.default_rng(1010)
lam, n, M, K = 3.0, 200, 20_000, 8      # 범주: 0,1,...,6 그리고 7 이상

def chi2_stat(x):
    """λ를 자료에서 추정한 포아송 적합도 통계량."""
    obs = np.bincount(np.minimum(x, K - 1), minlength=K).astype(float)
    lam_hat = x.mean()                   # ← 모수를 추정했다
    p = stats.poisson.pmf(np.arange(K - 1), lam_hat)
    p = np.r_[p, 1 - p.sum()]
    exp = n * p
    return ((obs - exp)**2 / exp).sum(), exp.min()

stat = np.empty(M)
emin = np.empty(M)
for i in range(M):
    stat[i], emin[i] = chi2_stat(rng.poisson(lam, n))

print(f"기대도수 최솟값의 평균  {emin.mean():.3f}")
print(f"통계량의 평균          {stat.mean():.4f}")
print(f"  (자유도가 df 이면 평균이 df 여야 한다: k-1={K - 1}, k-2={K - 2})\n")
for df, label in [(K - 1, "df = k-1  (조정 안 함)"), (K - 2, "df = k-2  (올바름)")]:
    print(f"{label:24s} 기각률 {np.mean(stat > stats.chi2.ppf(0.95, df)):.4f}")

기대도수 최솟값의 평균  6.665
통계량의 평균          5.9963
  (자유도가 df 이면 평균이 df 여야 한다: k-1=7, k-2=6)

df = k-1  (조정 안 함)       기각률 0.0273
df = k-2  (올바름)          기각률 0.0471

import numpy as np
from scipy import stats

rng = np.random.default_rng(4040)
M, n, K = 20_000, 200, 8
qs = np.linspace(0, 1, K + 1)[1:-1]      # 등확률 구간의 경계 분위수

stat = np.empty(M)
for i in range(M):
    x = rng.standard_normal(n)           # H0 가 참 (정규자료)
    mu, sd = x.mean(), x.std(ddof=1)     # ← 원자료에서 추정
    edges = np.r_[-np.inf, mu + sd * stats.norm.ppf(qs), np.inf]
    obs = np.histogram(x, bins=edges)[0].astype(float)
    exp = np.full(K, n / K)
    stat[i] = ((obs - exp)**2 / exp).sum()

print(f"통계량의 평균 {stat.mean():.4f}   "
      f"(k-1 = {K - 1},  k-1-2 = {K - 3})\n")
for df in [K - 1, K - 3]:
    print(f"  df={df} 를 쓰면 실제 수준 "
          f"{np.mean(stat > stats.chi2.ppf(0.95, df)):.4f}")

통계량의 평균 5.4538   (k-1 = 7,  k-1-2 = 5)

  df=7 를 쓰면 실제 수준 0.0195
  df=5 를 쓰면 실제 수준 0.0605

import numpy as np
from scipy import stats

rng = np.random.default_rng(3030)
M, n = 20_000, 100
p0 = np.array([.70, .10, .08, .07, .05])     # H0
p1 = np.array([.60, .10, .08, .12, .10])     # 대립: 작은 범주들이 커짐

def run(p_true, merge):
    obs = rng.multinomial(n, p_true, size=M).astype(float)
    if merge is None:
        exp = n * p0
    else:
        obs = np.column_stack([obs[:, s].sum(1) for s in merge])
        exp = n * np.array([p0[s].sum() for s in merge])
    k = len(exp)
    stat = ((obs - exp)**2 / exp).sum(1)
    return np.mean(stat > stats.chi2.ppf(0.95, k - 1)), exp.min(), k

plans = [
    ("병합 없음", None),
    ("뒤 2개 병합", [np.array([0]), np.array([1]), np.array([2]), np.array([3, 4])]),
    ("뒤 3개 병합", [np.array([0]), np.array([1]), np.array([2, 3, 4])]),
]
print(f"{'설정':>12s} {'k':>3s} {'E_min':>7s} {'수준':>8s} {'검정력':>8s}")
for label, merge in plans:
    lvl, emin, k = run(p0, merge)
    pw, _, _ = run(p1, merge)
    print(f"{label:>12s} {k:3d} {emin:7.2f} {lvl:8.4f} {pw:8.4f}")

          설정   k   E_min       수준      검정력
       병합 없음   5    5.00   0.0469   0.7053
     뒤 2개 병합   4    8.00   0.0505   0.7029
     뒤 3개 병합   3   10.00   0.0491   0.6034

# 이탈의 방향이 반대인 경우
p2 = np.array([.70, .10, .08, .12, .00])     # 4번 ↑, 5번 ↓
for label, merge in plans:
    pw, _, k = run(p2, merge)
    print(f"{label:>12s} (k={k})  검정력 {pw:.4f}")

       병합 없음 (k=5)  검정력 0.5988
     뒤 2개 병합 (k=4)  검정력 0.0471
     뒤 3개 병합 (k=3)  검정력 0.0470

import numpy as np
from scipy import stats

rng = np.random.default_rng(2020)
M = 40_000
cases = [
    (np.ones(5) / 5, 25), (np.ones(5) / 5, 10),
    (np.array([.90, .05, .04, .01]), 50),
    (np.array([.90, .05, .04, .01]), 100),
    (np.array([.90, .05, .04, .01]), 500),
    (np.array([.97, .01, .01, .01]), 100),
    (np.array([.97, .01, .01, .01]), 500),
    (np.array([.50, .30, .15, .04, .01]), 100),
]
print(f"{'p0':>32s} {'n':>5s} {'E_min':>7s} {'실제 수준':>10s}")
for p, n in cases:
    k, exp = len(p), n * p
    obs = rng.multinomial(n, p, size=M).astype(float)
    stat = ((obs - exp)**2 / exp).sum(1)
    lvl = np.mean(stat > stats.chi2.ppf(0.95, k - 1))
    print(f"{str(np.round(p, 3)):>32s} {n:5d} {exp.min():7.2f} {lvl:10.4f}")

                              p0     n   E_min      실제 수준
           [0.2 0.2 0.2 0.2 0.2]    25    5.00     0.0481
           [0.2 0.2 0.2 0.2 0.2]    10    2.00     0.0394
           [0.9  0.05 0.04 0.01]    50    0.50     0.0524
           [0.9  0.05 0.04 0.01]   100    1.00     0.0526
           [0.9  0.05 0.04 0.01]   500    5.00     0.0498
           [0.97 0.01 0.01 0.01]   100    1.00     0.0658
           [0.97 0.01 0.01 0.01]   500    5.00     0.0491
      [0.5  0.3  0.15 0.04 0.01]   100    1.00     0.0556

import numpy as np
from scipy import stats

rng = np.random.default_rng(5050)
M, n = 5_000, 200

print(f"{'범주 k':>7s} {'df':>4s} {'E_min':>7s} {'수준':>8s} "
      f"{'검정력(평균 이동)':>17s} {'검정력(두꺼운 꼬리)':>19s}")
for K in [2, 4, 8, 16, 32]:
    qs = np.linspace(0, 1, K + 1)[1:-1]
    edges = np.r_[-np.inf, stats.norm.ppf(qs), np.inf]   # 등확률 구간
    exp = np.full(K, n / K)
    crit = stats.chi2.ppf(0.95, K - 1)
    out = []
    for gen in [lambda: rng.standard_normal(n),                 # H0
                lambda: rng.standard_normal(n) + 0.3,           # 위치 이동
                lambda: rng.standard_t(4, n) / np.sqrt(2.0)]:   # 분산 1 로 맞춘 t(4)
        hit = 0
        for _ in range(M):
            obs = np.histogram(gen(), bins=edges)[0].astype(float)
            hit += ((obs - exp)**2 / exp).sum() > crit
        out.append(hit / M)
    print(f"{K:7d} {K - 1:4d} {n / K:7.2f} {out[0]:8.4f} "
          f"{out[1]:17.4f} {out[2]:19.4f}")

   범주 k   df   E_min       수준        검정력(평균 이동)         검정력(두꺼운 꼬리)
      2    1  100.00   0.0574            0.9190              0.0626
      4    3   50.00   0.0510            0.9174              0.7114
      8    7   25.00   0.0526            0.8690              0.6354
     16   15   12.50   0.0470            0.7708              0.4952
     32   31    6.25   0.0538            0.6234              0.3698

자유도 = (자유로운 칸의 수) - (자료로 추정한 모수의 수)
    │
    ├─ 적합도 검정 (범주 k 개)
    │      자유로운 칸 = k - 1        ← 합이 n 으로 고정
    │      모수 r 개 추정 → df = k - 1 - r
    │      ※ 원자료에서 추정했으면 이 값이 정확하지 않다 (연습문제 6)
    │
    ├─ 독립성·동질성 (r × c)
    │      df = (r-1)(c-1)
    │      ※ 행·열 합이 모두 고정되므로
    │
    ├─ 맥니마 검정 (2×2 대응)
    │      df = 1                     ← 불일치 쌍만 본다
    │
    └─ 코크런의 Q (k 개 조건, 대응)
           df = k - 1