기대 칸 도수와 타당성 조건¶

import numpy as np
from scipy import stats

rng = np.random.default_rng(606)

def mc_pvalue(obs, B, rng):
    """주변합을 고정한 몬테카를로 p 값."""
    chi2_obs = stats.chi2_contingency(obs, correction=False)[0]
    r = obs.sum(1).astype(int)
    c = obs.sum(0).astype(int)
    exp = np.outer(r, c) / obs.sum()
    rows = np.repeat(np.arange(len(r)), r)
    cols = np.repeat(np.arange(len(c)), c)
    cnt = 1
    for _ in range(B):
        rng.shuffle(cols)                     # 짝짓기를 무작위로
        t = np.zeros_like(obs, dtype=float)
        np.add.at(t, (rows, cols), 1)
        cnt += ((t - exp)**2 / exp).sum() >= chi2_obs - 1e-9
    return cnt / (B + 1)

M = 2_000
print(f"{'상황':>18s} {'E_min':>6s} {'χ² 근사':>9s} {'몬테카를로':>11s}")
setups = [(np.array([[.25, .25], [.25, .25]]), 20, "2×2 균등 n=20"),
          (np.outer([.6, .4], [.85, .15]), 40, "2×2 치우침 n=40"),
          (np.outer([.5, .3, .2], [.7, .2, .1]), 40, "3×3 치우침 n=40")]
for P, n, label in setups:
    a = b = used = 0
    for _ in range(M):
        obs = rng.multinomial(n, P.ravel()).reshape(P.shape).astype(float)
        if (obs.sum(0) == 0).any() or (obs.sum(1) == 0).any():
            continue
        used += 1
        a += stats.chi2_contingency(obs, correction=False)[1] < 0.05
        b += mc_pvalue(obs, 499, rng) < 0.05
    print(f"{label:>18s} {n * P.min():6.2f} {a / used:9.4f} {b / used:11.4f}")

                상황  E_min     χ² 근사       몬테카를로
       2×2 균등 n=20   5.00    0.0605      0.0240
      2×2 치우침 n=40   2.40    0.0410      0.0190
      3×3 치우침 n=40   0.80    0.0527      0.0542

import numpy as np
from scipy import stats

rng = np.random.default_rng(707)
M = 5_000

print(f"{'군집 수':>7s} {'크기':>5s} {'ICC':>5s} {'DEFF':>6s} "
      f"{'보정 없음':>10s} {'DEFF 로 나눔':>13s}")
for n_cluster, m, icc in [(40, 5, 0.0), (40, 5, 0.1), (40, 5, 0.3),
                          (20, 10, 0.3), (10, 20, 0.3)]:
    deff = 1 + (m - 1) * icc
    # 베타·이항으로 군집 내 상관을 만든다: icc = 1/(1+a+b)
    a = b = 0.5 * (1 / icc - 1) if icc > 0 else None
    r_raw = r_adj = 0
    for _ in range(M):
        t = np.zeros((2, 2))
        for g in range(2):                       # 두 군
            for _ in range(n_cluster // 2):
                p_c = rng.beta(a, b) if icc > 0 else 0.5
                k = rng.binomial(m, p_c)
                t[g, 0] += k
                t[g, 1] += m - k
        if t.sum(0).min() > 0:
            chi2 = stats.chi2_contingency(t, correction=False)[0]
            r_raw += stats.chi2.sf(chi2, 1) < 0.05
            r_adj += stats.chi2.sf(chi2 / deff, 1) < 0.05
    print(f"{n_cluster:7d} {m:5d} {icc:5.1f} {deff:6.2f} "
          f"{r_raw / M:10.4f} {r_adj / M:13.4f}")

   군집 수    크기   ICC   DEFF      보정 없음     DEFF 로 나눔
     40     5   0.0   1.00     0.0606        0.0606
     40     5   0.1   1.40     0.1132        0.0530
     40     5   0.3   2.20     0.1984        0.0556
     20    10   0.3   3.70     0.3218        0.0534
     10    20   0.3   6.70     0.4744        0.0566

import numpy as np
from scipy import stats

rng = np.random.default_rng(808)
M = 3_000

print("① 수준 (독립이 참, 주변이 균등)")
print(f"{'표 크기':>8s} {'n':>6s} {'칸당 기대':>10s} {'0인 칸 평균':>12s} "
      f"{'기각률':>8s}")
for (r, c), n in [((5, 5), 50), ((5, 5), 100), ((5, 5), 250),
                  ((10, 10), 100), ((10, 10), 500), ((20, 20), 400)]:
    P = np.ones(r * c) / (r * c)
    rej = used = 0
    zeros = []
    for _ in range(M):
        obs = rng.multinomial(n, P).reshape(r, c).astype(float)
        zeros.append((obs == 0).sum())
        if (obs.sum(0) == 0).any() or (obs.sum(1) == 0).any():
            continue
        used += 1
        rej += stats.chi2_contingency(obs, correction=False)[1] < 0.05
    print(f"{f'{r}×{c}':>8s} {n:6d} {n / (r * c):10.2f} "
          f"{np.mean(zeros):12.1f} {rej / used:8.4f}")

① 수준 (독립이 참, 주변이 균등)
    표 크기      n      칸당 기대      0인 칸 평균      기각률
     5×5     50       2.00          3.3   0.0437
     5×5    100       4.00          0.4   0.0507
     5×5    250      10.00          0.0   0.0483
   10×10    100       1.00         36.6   0.0353
   10×10    500       5.00          0.7   0.0510
   20×20    400       1.00        146.9   0.0460

print("\n② 검정력 (같은 연관을 잘게 쪼갠 표들, n=200)")
base = np.array([[.30, .20], [.20, .30]])       # 2×2, V=0.20
for k in [1, 2, 3, 5]:
    P = np.kron(base, np.ones((k, k)) / (k * k))  # 행·열을 k 개로 쪼갬
    r, c = P.shape
    p_row, p_col = P.sum(1), P.sum(0)
    v = np.sqrt((((P - np.outer(p_row, p_col))**2)
                 / np.outer(p_row, p_col)).sum() / (min(r, c) - 1))
    rej = used = 0
    for _ in range(M):
        obs = rng.multinomial(200, P.ravel()).reshape(r, c).astype(float)
        if (obs.sum(0) == 0).any() or (obs.sum(1) == 0).any():
            continue
        used += 1
        rej += stats.chi2_contingency(obs, correction=False)[1] < 0.05
    print(f"  {r}×{c} 표 (V={v:.4f}): 칸당 기대 {200 / (r * c):5.2f}, "
          f"검정력 {rej / used:.4f}")

② 검정력 (같은 연관을 잘게 쪼갠 표들, n=200)
  2×2 표 (V=0.2000): 칸당 기대 50.00, 검정력 0.8097
  4×4 표 (V=0.1155): 칸당 기대 12.50, 검정력 0.4413
  6×6 표 (V=0.0894): 칸당 기대  5.56, 검정력 0.2793
  10×10 표 (V=0.0667): 칸당 기대  2.00, 검정력 0.1393

① 독립성 (가장 중요, 검정으로 확인 불가)
    자료가 어떻게 수집되었는가?
    ├─ 같은 대상의 반복 측정 → McNemar / 코크런 Q
    ├─ 군집 표집 → DEFF 보정 / 혼합모형 (연습문제 8)
    └─ 단순 확률표본 → 통과
          ↓
② 범주의 정의
    ├─ 상호배타적이고 빠짐없는가
    ├─ 자료를 보고 범주를 정하지 않았는가
    └─ 구조적 0 이 있는가 (연습문제 9)
          ↓
③ 기대도수
    ├─ 모두 5 이상 → 카이제곱
    ├─ 일부 작지만 주변이 고르다 → 카이제곱도 대체로 무방
    └─ 작고 주변이 심하게 치우침 → 정확검정 / 몬테카를로
          ↓
④ 표본크기
    희소한 큰 표에서 유의하지 않으면 검정력을 계산해 본다