적합도 검정¶

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

# 관측도수와, 고르게 나온다고 볼 때의 기대도수
observed_counts = np.array([4, 13, 7])
expected_counts = np.ones(3) * observed_counts.mean()
degrees_of_freedom = observed_counts.shape[0] - 1

# 검정통계량과 p-값 계산
chi_square_statistic = np.sum((observed_counts - expected_counts) ** 2 / expected_counts)
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}\n")

# 그림 준비
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, color='b', 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, color='b', 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, color='r', linewidth=3)

# 오른쪽을 칠한다. 이 넓이가 곧 p-값이다.
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, color='r', alpha=0.1)

# 화살표로 p-값을 가리킨다
annotation_xy = ((12.5 + 15.0) / 2, 0.01)
annotation_xytext = (16.5, 0.10)
arrow_properties = dict(color='k', width=0.2, headwidth=8)
ax.annotate(f'p-value = {p_value:.02%}', annotation_xy, xytext=annotation_xytext,
            fontsize=15, arrowprops=arrow_properties)

# 축과 테두리를 다듬는다
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.show()

Chi-square Statistic = 5.2500
p-value = 0.0724

from scipy import stats

# 결과별 관측도수: 승, 패, 무
observed_frequencies = [4, 13, 7]

# 세 결과가 고르게 나온다고 볼 때의 기대도수
total_games = sum(observed_frequencies)
expected_frequencies = [total_games / 3] * 3

# 카이제곱 적합도 검정
chi_square_statistic, p_value = stats.chisquare(f_obs=observed_frequencies, f_exp=expected_frequencies)

# 결과 출력
print(f"{chi_square_statistic = }")
print(f"{p_value = }")

chi_square_statistic = 5.25
p_value = 0.07243975703425146

import numpy as np
import scipy.stats as stats

def main():
    p_0 = 1/6
    p_hat = 17 / 60
    n = 60

    # 눈 3 하나만 놓고 보는 일표본 비율 z-검정
    statistic = (p_hat - p_0) / np.sqrt(p_0 * (1 - p_0) / n)
    p_value = stats.norm().sf(abs(statistic)) * 2

    print(f"{statistic = :.02f}")
    print(f"{p_value   = :.02%}")

if __name__ == "__main__":
    main()

statistic = 2.42
p_value   = 1.53%

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
M = 200_000
n, k, p = 60, 6, 1 / 6
z_crit = stats.norm.ppf(0.975)

# 공정한 주사위를 60번 굴리는 일을 20만 번 되풀이한다.
X = rng.multinomial(n, [p] * k, size=M)
Z = np.abs(X - n * p) / np.sqrt(n * p * (1 - p))
any_sig = (Z > z_crit).any(axis=1).mean()

# 눈 하나만 보았을 때 유의해질 확률 q — 주변적으로는 이항분포다.
b = stats.binom(n, p)
q = b.cdf(4) + b.sf(15)

print(f"임계값  |O_i - 10| > {z_crit * np.sqrt(n * p * (1 - p)):.4f}"
      f"   곧 O_i <= 4 또는 O_i >= 16")
print(f"눈 하나가 유의할 확률  q          {q:.4f}")
print(f"여섯이 독립이라면  1 - (1-q)^6    {1 - (1 - q) ** k:.4f}")
print(f"본페로니 상한      6q            {k * q:.4f}")
print(f"여섯 중 적어도 하나 (모의)        {any_sig:.4f}   "
      f"오차 {np.sqrt(any_sig * (1 - any_sig) / M):.4f}")
print(f"corr(O_1, O_2)  모의 {np.corrcoef(X[:, 0], X[:, 1])[0, 1]:.4f}"
      f"   이론 -p/(1-p) = {-p / (1 - p):.4f}")

임계값  |O_i - 10| > 5.6579   곧 O_i <= 4 또는 O_i >= 16
눈 하나가 유의할 확률  q          0.0540
여섯이 독립이라면  1 - (1-q)^6    0.2835
본페로니 상한      6q            0.3242
여섯 중 적어도 하나 (모의)        0.2726   오차 0.0010
corr(O_1, O_2)  모의 -0.1993   이론 -p/(1-p) = -0.2000

import numpy as np

def main():
    observed = np.array([5, 7, 17, 14, 8, 9])
    expected = np.array([10] * 6)
    # 여섯 눈의 어긋남을 **한 숫자로 모은다**.
    # 눈 하나만 보던 앞의 검정과 여기서 갈린다.
    statistic = np.sum((observed - expected)**2 / expected)
    print(f'{statistic = }')

if __name__ == "__main__":
    main()

statistic = 10.4

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

def main():
    """적합도 검정통계량을 정의대로 구하고 p-값을 그림으로 보인다."""
    # 주사위를 60번 굴린 결과라고 하자. 고른 주사위라면 눈마다 10번씩 나온다.
    observed = np.array([5, 7, 17, 14, 8, 9])
    expected = np.array([10] * 6)

    # 자유도는 범주 수에서 1을 뺀다. 도수의 합이 60으로 묶여 있어
    # 다섯 칸을 알면 나머지 한 칸이 자동으로 정해지기 때문이다.
    df = observed.shape[0] - 1

    statistic = np.sum((observed - expected)**2 / expected)
    p_value = stats.chi2(df).sf(statistic)
    print(f"{statistic = :.02f}")
    print(f"{p_value    = :.02%}")

    # 파랑은 통계량보다 작은 쪽, 빨강은 그보다 큰 쪽이다.
    # 적합도 검정은 언제나 우측검정이므로 빨간 넓이가 곧 p-값이다.
    _, ax = plt.subplots(figsize=(12, 4))

    x = np.linspace(0, statistic)
    y = stats.chi2(df).pdf(x)
    ax.plot(x, y, color='b', linewidth=3)

    x = np.concatenate([[0], x, [statistic], [0]])
    y = np.concatenate([[0], y, [0], [0]])
    ax.fill(x, y, color='b', alpha=0.1)

    x = np.linspace(statistic, 20, 100)
    y = stats.chi2(df).pdf(x)
    ax.plot(x, y, color='r', linewidth=3)

    x = np.concatenate([[statistic], x, [20], [statistic]])
    y = np.concatenate([[0], y, [0], [0]])
    ax.fill(x, y, color='r', alpha=0.1)

    xy = ((12.5 + 15.0) / 2, 0.01)
    xytext = (16.5, 0.10)
    arrowprops = dict(color='k', width=0.2, headwidth=8)
    ax.annotate(f'{p_value = :.02%}', xy, xytext=xytext, fontsize=15, arrowprops=arrowprops)

    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.show()

if __name__ == "__main__":
    main()

statistic = 10.40
p_value    = 6.47%

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

def main():
    """앞의 손계산을 scipy 한 줄로 대신한다. 값이 같아야 한다."""
    observed = np.array([5, 7, 17, 14, 8, 9])
    expected = np.array([10] * 6)
    df = observed.shape[0] - 1

    # chisquare 는 통계량과 p-값을 한꺼번에 돌려준다. 자유도는 알아서 정한다.
    statistic, p_value = stats.chisquare(observed, f_exp=expected)
    print(f"{statistic = :.02f}")
    print(f"{p_value    = :.02%}")

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

    x = np.linspace(0, statistic)
    y = stats.chi2(df).pdf(x)
    ax.plot(x, y, color='b', linewidth=3)

    x = np.concatenate([[0], x, [statistic], [0]])
    y = np.concatenate([[0], y, [0], [0]])
    ax.fill(x, y, color='b', alpha=0.1)

    x = np.linspace(statistic, 20, 100)
    y = stats.chi2(df).pdf(x)
    ax.plot(x, y, color='r', linewidth=3)

    x = np.concatenate([[statistic], x, [20], [statistic]])
    y = np.concatenate([[0], y, [0], [0]])
    ax.fill(x, y, color='r', alpha=0.1)

    xy = ((12.5 + 15.0) / 2, 0.01)
    xytext = (16.5, 0.10)
    arrowprops = dict(color='k', width=0.2, headwidth=8)
    ax.annotate(f'{p_value = :.02%}', xy, xytext=xytext, fontsize=15, arrowprops=arrowprops)

    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.show()

if __name__ == "__main__":
    main()

statistic = 10.40
p_value    = 6.47%

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
M = 200_000
n, k = 60, 6
obs = 10.4

# H0 가 참인 세상에서 주사위를 60번 굴리는 일을 20만 번 되풀이한다.
X = rng.multinomial(n, np.ones(k) / k, size=M)
chi2_sim = ((X - n / k) ** 2 / (n / k)).sum(axis=1)

print(f"모의 {M}회")
print(f"  E[chi2]  이론 k-1 = {k - 1}        모의 {chi2_sim.mean():.4f}   "
      f"오차 {chi2_sim.std() / np.sqrt(M):.4f}")
print(f"  p-값     카이제곱 근사 {stats.chi2.sf(obs, k - 1):.4f}   "
      f"모의 {(chi2_sim >= obs - 1e-9).mean():.4f}   "
      f"오차 {np.sqrt(0.066 * 0.934 / M):.4f}")

모의 200000회
  E[chi2]  이론 k-1 = 5        모의 5.0039   오차 0.0070
  p-값     카이제곱 근사 0.0647   모의 0.0659   오차 0.0006

import numpy as np
from scipy import stats

obs = np.array([109, 65, 22, 3, 1])
k = np.arange(5)
n = obs.sum()
lam = (obs * k).sum() / n                 # ① 모수를 자료에서 추정
print(f"n = {n},  λ̂ = {lam:.4f}\n")

p = stats.poisson.pmf(np.arange(4), lam)  # ② 기대도수 (마지막은 4 이상)
p = np.r_[p, 1 - p.sum()]
exp = n * p
print("범주        0       1       2       3      4+")
print(f"관측  {obs[0]:7d} {obs[1]:7d} {obs[2]:7d} {obs[3]:7d} {obs[4]:7d}")
print("기대  " + " ".join(f"{v:7.3f}" for v in exp))

chi2 = ((obs - exp)**2 / exp).sum()
print(f"\n병합 전:  χ² = {chi2:.4f},  E_min = {exp.min():.3f}  ← 5 미만")
for df in [len(obs) - 1, len(obs) - 2]:
    print(f"   df={df} 이면 p = {stats.chi2.sf(chi2, df):.4f}")

# ③ 기대도수가 작은 꼬리를 병합한다
obs2 = np.r_[obs[:3], obs[3:].sum()]
p2 = stats.poisson.pmf(np.arange(3), lam)
p2 = np.r_[p2, 1 - p2.sum()]
exp2 = n * p2
chi2b = ((obs2 - exp2)**2 / exp2).sum()
df2 = len(obs2) - 1 - 1                   # 범주 4개, λ 하나 추정
print(f"\n병합 후 (0, 1, 2, 3+)")
print(f"  관측 {obs2.tolist()},  기대 {np.round(exp2, 3).tolist()}")
print(f"  χ² = {chi2b:.4f},  E_min = {exp2.min():.3f},  df = {df2},  "
      f"p = {stats.chi2.sf(chi2b, df2):.4f}")
print(f"  표준화 잔차 {np.round((obs2 - exp2) / np.sqrt(exp2), 3).tolist()}")

n = 200,  λ̂ = 0.6100

범주        0       1       2       3      4+
관측      109      65      22       3       1
기대  108.670  66.289  20.218   4.111   0.712

병합 전:  χ² = 0.5999,  E_min = 0.712  ← 5 미만
   df=4 이면 p = 0.9631
   df=3 이면 p = 0.8964

병합 후 (0, 1, 2, 3+)
  관측 [109, 65, 22, 4],  기대 [108.67, 66.289, 20.218, 4.823]
  χ² = 0.3235,  E_min = 4.823,  df = 2,  p = 0.8506
  표준화 잔차 [0.032, -0.158, 0.396, -0.375]

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

def multinomial_exact(obs, p):
    """모든 도수 벡터를 열거하는 다항 정확검정. k 와 n 이 작을 때만 쓴다."""
    obs = np.asarray(obs, float)
    p = np.asarray(p, float)
    n, k = int(obs.sum()), len(p)
    exp = n * p
    t_obs = ((obs - exp)**2 / exp).sum()
    log_n_fact = lgamma(n + 1)
    total = 0.0

    def walk(i, left, counts, log_p):
        nonlocal total
        if i == k - 1:
            c = np.array(counts + [left], float)
            lp = log_p + left * np.log(p[i]) - lgamma(left + 1)
            if ((c - exp)**2 / exp).sum() >= t_obs - 1e-9:
                total += np.exp(log_n_fact + lp)
            return
        for x in range(left + 1):
            walk(i + 1, left - x, counts + [x],
                 log_p + x * np.log(p[i]) - lgamma(x + 1))

    walk(0, n, [], 0.0)
    return t_obs, total

cases = [(np.array([2, 3, 7]), [1 / 3] * 3),
         (np.array([0, 4, 8]), [1 / 3] * 3),
         (np.array([5, 5, 2, 0]), [0.25] * 4),
         (np.array([9, 6, 3, 2]), [0.4, 0.3, 0.2, 0.1])]
for obs, p in cases:
    t, p_exact = multinomial_exact(obs, p)
    p_approx = stats.chi2.sf(t, len(p) - 1)
    print(f"관측 {str(obs.tolist()):>14s}  χ² = {t:6.4f}   "
          f"정확 p = {p_exact:.4f}   근사 p = {p_approx:.4f}")

관측      [2, 3, 7]  χ² = 3.5000   정확 p = 0.2672   근사 p = 0.1738
관측      [0, 4, 8]  χ² = 8.0000   정확 p = 0.0172   근사 p = 0.0183
관측   [5, 5, 2, 0]  χ² = 6.0000   정확 p = 0.1300   근사 p = 0.1116
관측   [9, 6, 3, 2]  χ² = 0.3750   정확 p = 0.9634   근사 p = 0.9454

import numpy as np
from scipy import stats
from scipy.optimize import brentq

def w_upper(chi2, df, n, alpha=0.05):
    """효과크기 w 의 단측 95% 상한 (비중심 카이제곱 반전)."""
    if stats.ncx2.cdf(chi2, df, 0.0) <= alpha:
        return 0.0
    lam = brentq(lambda l: stats.ncx2.cdf(chi2, df, l) - alpha, 0.0, 1e7)
    return np.sqrt(lam / n)

shape = np.array([8, 12, 9, 11, 10, 10]) / 60     # 이탈의 모양을 고정
p0, df = np.ones(6) / 6, 5
print(f"{'n':>7s} {'χ²':>10s} {'p':>9s} {'ŵ':>8s} {'w 의 95% 상한':>15s}")
for n in [60, 600, 6_000, 60_000]:
    obs = n * shape
    exp = n * p0
    chi2 = ((obs - exp)**2 / exp).sum()
    print(f"{n:7d} {chi2:10.4f} {stats.chi2.sf(chi2, df):9.4f} "
          f"{np.sqrt(chi2 / n):8.4f} {w_upper(chi2, df, n):15.4f}")

      n         χ²         p        ŵ      w 의 95% 상한
     60     1.0000    0.9626   0.1291          0.0000
    600    10.0000    0.0752   0.1291          0.1738
   6000   100.0000    0.0000   0.1291          0.1479
  60000  1000.0000    0.0000   0.1291          0.1356

① 자료의 형태를 본다
    ├─ 범주형 (도수) ──────────→ 카이제곱 적합도
    └─ 연속형 ────────────────→ 구간화하지 말고 전용 검정
                                 (KS, 앤더슨·달링, 샤피로·윌크)
          ↓
② 귀무분포를 정한다
    ├─ 모수가 완전히 주어짐 ──→ df = k - 1
    └─ 모수를 추정 ──────────→ df = k - 1 - r
                                 ※ 원자료 추정이면 근사일 뿐
          ↓
③ 기대도수를 계산하고 확인한다
    └─ 작은 칸은 사전에 정한 규칙대로 병합
          ↓
④ 검정 + 잔차 + 효과크기
          ↓
⑤ 해석: "기각 못 함"은 "맞음"이 아니다