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

from scipy import stats

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

# H0(균등분포) 아래의 기대도수.
# f_exp의 합은 f_obs의 합과 같아야 한다. 다르면 scipy가 오류를 낸다.
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

from math import factorial
from scipy import stats

n, k = 24, 3
obs_stat = 5.25

# H0 가 참인 세상에서 가능한 (승, 패, 무) 조합을 전부 열거한다.
dist = {}
for a in range(n + 1):
    for b in range(n + 1 - a):
        c = n - a - b
        prob = (factorial(n) / (factorial(a) * factorial(b) * factorial(c))
                / k ** n)
        stat = sum((x - n / k) ** 2 / (n / k) for x in (a, b, c))
        dist[round(stat, 9)] = dist.get(round(stat, 9), 0.0) + prob

crit = stats.chi2.ppf(0.95, k - 1)
print(f"확률의 총합              {sum(dist.values()):.12f}")
print(f"통계량이 가질 수 있는 값   {len(dist)} 개")
print(f"E[chi2]   정확 {sum(s * v for s, v in dist.items()):.6f}   이론 k-1 = {k - 1}")
print(f"p-값      정확 {sum(v for s, v in dist.items() if s >= obs_stat - 1e-9):.6f}"
      f"   근사 {stats.chi2.sf(obs_stat, k - 1):.6f}")
print(f"명목 5% 임계값 {crit:.4f} 에서의 실제 제1종 오류율 "
      f"{sum(v for s, v in dist.items() if s >= crit - 1e-9):.6f}")

확률의 총합              1.000000000000
통계량이 가질 수 있는 값   44 개
E[chi2]   정확 2.000000   이론 k-1 = 2
p-값      정확 0.080040   근사 0.072440
명목 5% 임계값 5.9915 에서의 실제 제1종 오류율 0.049902

statistic, p = stats.chisquare(f_obs=[4, 13, 7])
print(f"{statistic = }, {p = }")

statistic = 5.25, p = 0.07243975703425146

from scipy import stats

for obs in [[4, 13, 7], [1, 2, 3, 4], [7, 7, 7, 7, 7], [0, 5, 11]]:
    k = len(obs)
    uniform = [sum(obs) / k] * k          # scipy 가 속으로 넣는 값
    a = stats.chisquare(obs)              # f_exp 생략
    b = stats.chisquare(obs, f_exp=uniform)
    print(f"{str(obs):>18s}  n/k = {sum(obs) / k!r:<18s} "
          f"chi2 {a.statistic:8.4f}   같은가 "
          f"{a.statistic == b.statistic and a.pvalue == b.pvalue}")

        [4, 13, 7]  n/k = 8.0                chi2   5.2500   같은가 True
      [1, 2, 3, 4]  n/k = 2.5                chi2   2.0000   같은가 True
   [7, 7, 7, 7, 7]  n/k = 7.0                chi2   0.0000   같은가 True
        [0, 5, 11]  n/k = 5.333333333333333  chi2  11.3750   같은가 True

n = 300
proportions = [0.4, 0.3, 0.3]
# 비율이 아니라 **도수**를 넘겨야 한다. 비율 [0.4, 0.3, 0.3]을 그대로 넣으면
# 합이 f_obs의 합과 달라 오류가 나거나, 운이 나쁘면 엉뚱한 값이 나온다.
expected = [n * p for p in proportions]
stat, pval = stats.chisquare(f_obs=[130, 85, 85], f_exp=expected)
print(f"{stat = :.4f}, {pval = :.4f}")

stat = 1.3889, pval = 0.4994

from scipy import stats

n = 300
proportions = [0.4, 0.3, 0.3]
obs = [130, 85, 85]

# ① 도수로 바꾸어 넘긴다 — 올바른 사용
stat, pval = stats.chisquare(f_obs=obs, f_exp=[n * p for p in proportions])
print(f"도수로 넘김   chi2 = {stat:.6f},  p = {pval:.6f}")

# ② 비율을 그대로 넘긴다 — scipy 가 막아 준다
try:
    stats.chisquare(f_obs=obs, f_exp=proportions)
except ValueError as e:
    print(f"비율로 넘김   ValueError: {str(e).splitlines()[0][:60]}...")

# 막지 않았다면 어떤 값이 나왔을까. 보기 1 (1) 의 식으로 예측한다.
c = sum(proportions)
predicted = (n / c) * stat + (n - c) ** 2 / c
actual = sum((o - c * p) ** 2 / (c * p) for o, p in zip(obs, proportions))
print(f"식의 예측 (n/c)chi2 + (n-c)^2/c = {predicted:.4f}")
print(f"직접 계산                        = {actual:.4f}")

도수로 넘김   chi2 = 1.388889,  p = 0.499352
비율로 넘김   ValueError: For each axis slice, the sum of the observed frequencies mus...
식의 예측 (n/c)chi2 + (n-c)^2/c = 89817.6667
직접 계산                        = 89817.6667

from scipy import stats
stat, p = stats.chisquare(f_obs=[90, 70, 40], f_exp=[100, 60, 40])
print(f"{stat = :.4f}, {p = :.4f}")

stat = 2.6667, p = 0.2636

from scipy import stats

obs = [90, 70, 40]
p = [0.5, 0.3, 0.2]
n = sum(obs)
exp = [n * q for q in p]

print(f"{'lambda_':>8s} {'이름':>14s} {'통계량':>10s} {'p 값':>9s}")
for lam, name in [(1, "피어슨 X²"), (0, "우도비 G²"), (-0.5, "프리먼·튜키"),
                  (2 / 3, "크레시·리드"), (-1, "수정 우도비"),
                  (-2, "네이만 수정")]:
    r = stats.power_divergence(f_obs=obs, f_exp=exp, lambda_=lam)
    print(f"{lam:8.3f} {name:>14s} {r.statistic:10.4f} {r.pvalue:9.4f}")

 lambda_             이름        통계량       p 값
   1.000         피어슨 X²     2.6667    0.2636
   0.000         우도비 G²     2.6162    0.2703
  -0.500         프리먼·튜키     2.5941    0.2733
   0.667         크레시·리드     2.6489    0.2659
  -1.000         수정 우도비     2.5740    0.2761
  -2.000         네이만 수정     2.5397    0.2809

for name in ["pearson", "log-likelihood", "freeman-tukey",
             "cressie-read", "neyman", "mod-log-likelihood"]:
    r = stats.power_divergence(f_obs=obs, f_exp=exp, lambda_=name)
    print(f"{name:>20s}  통계량 {r.statistic:.4f}")

             pearson  통계량 2.6667
      log-likelihood  통계량 2.6162
       freeman-tukey  통계량 2.5941
        cressie-read  통계량 2.6489
              neyman  통계량 2.5397
  mod-log-likelihood  통계량 2.5740

import numpy as np
from scipy import stats

table = np.array([[30, 20], [15, 35]])          # 2×2 분할표

r1 = stats.chi2_contingency(table, correction=False)
print(f"chi2_contingency (올바름)    χ²={r1[0]:.4f}  df={r1[2]}  "
      f"p={r1[1]:.6f}")
print(f"  기대도수\n{np.round(r1[3], 2)}")

# 잘못된 사용: 분할표를 1차원으로 펴서 적합도 검정으로
r2 = stats.chisquare(f_obs=table.ravel())
print(f"\nchisquare (잘못)             χ²={r2.statistic:.4f}  df=3  "
      f"p={r2.pvalue:.6f}")
print(f"  → 이것은 '네 칸이 모두 같은가'를 검정한 것이다")

odds, p_f = stats.fisher_exact(table)
print(f"\nfisher_exact                 오즈비={odds:.4f}  p={p_f:.6f}")

chi2_contingency (올바름)    χ²=9.0909  df=1  p=0.002569
  기대도수
[[22.5 27.5]
 [22.5 27.5]]

chisquare (잘못)             χ²=10.0000  df=3  p=0.018566
  → 이것은 '네 칸이 모두 같은가'를 검정한 것이다

fisher_exact                 오즈비=3.5000  p=0.004635

for corr in [True, False]:
    r = stats.chi2_contingency(table, correction=corr)
    print(f"correction={str(corr):5s}:  χ²={r[0]:.4f}  p={r[1]:.6f}")

correction=True :  χ²=7.9192  p=0.004891
correction=False:  χ²=9.0909  p=0.002569

import numpy as np
from scipy import stats

obs = [90, 70, 40]
p_exp = [0.5, 0.3, 0.2]           # 비율만 주어진 경우

# ① 틀린 방법: 비율을 그대로 f_exp 에 넣는다
try:
    stats.chisquare(f_obs=obs, f_exp=p_exp)
except ValueError as e:
    print(f"① 비율을 그대로:  ValueError\n   {e}\n")

# ② 올바른 방법: n 을 곱해 도수로 바꾼다
n = sum(obs)
exp = np.array(p_exp) * n
r = stats.chisquare(f_obs=obs, f_exp=exp)
print(f"② n 을 곱한 뒤:  χ²={r.statistic:.4f}  p={r.pvalue:.4f}")
print(f"   f_exp = {exp.tolist()},  합 = {exp.sum():.1f} = n\n")

# ③ 비율의 합이 1 이 아닐 때는 먼저 정규화한다
p_raw = [5, 3, 2]                  # 5:3:2 비율
p_norm = np.array(p_raw) / sum(p_raw)
exp3 = p_norm * n
r3 = stats.chisquare(f_obs=obs, f_exp=exp3)
print(f"③ 5:3:2 비율:  정규화 {np.round(p_norm, 4).tolist()}")
print(f"   f_exp = {exp3.tolist()},  χ²={r3.statistic:.4f}  p={r3.pvalue:.4f}")

① 비율을 그대로:  ValueError
   For each axis slice, the sum of the observed frequencies must agree with the sum of the expected frequencies to a relative tolerance of 1e-08, but the percent differences are:
199.0

② n 을 곱한 뒤:  χ²=2.6667  p=0.2636
   f_exp = [100.0, 60.0, 40.0],  합 = 200.0 = n

③ 5:3:2 비율:  정규화 [0.5, 0.3, 0.2]
   f_exp = [100.0, 60.0, 40.0],  χ²=2.6667  p=0.2636

def gof(obs, ratios=None, n_estimated=0):
    """기대비율(정규화 여부 무관)을 받아 적합도 검정을 수행한다."""
    obs = np.asarray(obs, float)
    n, k = obs.sum(), len(obs)
    if ratios is None:
        p = np.full(k, 1 / k)
    else:
        p = np.asarray(ratios, float)
        p = p / p.sum()                       # 알아서 정규화
    exp = n * p
    r = stats.chisquare(f_obs=obs, f_exp=exp, ddof=n_estimated)
    return {"chi2": r.statistic, "df": k - 1 - n_estimated,
            "p": r.pvalue, "E_min": exp.min()}

for ratios in [None, [0.5, 0.3, 0.2], [5, 3, 2], [50, 30, 20]]:
    out = gof(obs, ratios)
    print(f"ratios={str(ratios):>18s} → χ²={out['chi2']:7.4f}  "
          f"p={out['p']:.4f}  E_min={out['E_min']:.1f}")

ratios=              None → χ²=19.0000  p=0.0001  E_min=66.7
ratios=   [0.5, 0.3, 0.2] → χ²= 2.6667  p=0.2636  E_min=40.0
ratios=         [5, 3, 2] → χ²= 2.6667  p=0.2636  E_min=40.0
ratios=      [50, 30, 20] → χ²= 2.6667  p=0.2636  E_min=40.0

import numpy as np
from scipy import stats

rng = np.random.default_rng(4242)
n = 200

def compare(x, label):
    # ① 카이제곱: 등확률 구간 8개로 나눈다 (μ, σ 를 자료에서 추정)
    k = 8
    mu, sd = x.mean(), x.std(ddof=1)
    qs = np.linspace(0, 1, k + 1)[1:-1]
    edges = np.r_[-np.inf, mu + sd * stats.norm.ppf(qs), np.inf]
    obs = np.histogram(x, bins=edges)[0]
    chi2 = stats.chisquare(f_obs=obs, f_exp=np.full(k, len(x) / k),
                           ddof=2)
    print(f"{label}")
    print(f"  카이제곱(구간 8, ddof=2)  통계량 {chi2.statistic:7.4f}  "
          f"p={chi2.pvalue:.4f}")
    print(f"  샤피로·윌크              W={stats.shapiro(x).statistic:.4f}  "
          f"p={stats.shapiro(x).pvalue:.4f}")
    ad = stats.anderson(x, dist='norm')
    crit5 = ad.critical_values[list(ad.significance_level).index(5.0)]
    print(f"  앤더슨·달링              A²={ad.statistic:.4f}  "
          f"5% 임계값 {crit5:.4f}  "
          f"{'기각' if ad.statistic > crit5 else '기각 못 함'}")
    ks = stats.kstest((x - mu) / sd, 'norm')
    print(f"  콜모고로프·스미르노프      D={ks.statistic:.4f}  "
          f"p={ks.pvalue:.4f}  (모수 추정을 무시한 값이라 보수적)")

compare(rng.standard_normal(n), "① 정규자료")
print()
compare(rng.standard_t(4, n), "② t(4) — 꼬리가 두껍다")
print()
compare(rng.lognormal(0, 0.6, n), "③ 로그정규 — 오른쪽으로 치우침")

① 정규자료
  카이제곱(구간 8, ddof=2)  통계량  7.6000  p=0.1797
  샤피로·윌크              W=0.9961  p=0.8959
  앤더슨·달링              A²=0.2088  5% 임계값 0.7720  기각 못 함
  콜모고로프·스미르노프      D=0.0313  p=0.9862  (모수 추정을 무시한 값이라 보수적)

② t(4) — 꼬리가 두껍다
  카이제곱(구간 8, ddof=2)  통계량 18.4000  p=0.0025
  샤피로·윌크              W=0.8983  p=0.0000
  앤더슨·달링              A²=1.6305  5% 임계값 0.7720  기각
  콜모고로프·스미르노프      D=0.0767  p=0.1806  (모수 추정을 무시한 값이라 보수적)

③ 로그정규 — 오른쪽으로 치우침
  카이제곱(구간 8, ddof=2)  통계량 40.8800  p=0.0000
  샤피로·윌크              W=0.9027  p=0.0000
  앤더슨·달링              A²=4.2416  5% 임계값 0.7720  기각
  콜모고로프·스미르노프      D=0.0996  p=0.0352  (모수 추정을 무시한 값이라 보수적)

x = rng.standard_t(4, 300)
mu, sd = x.mean(), x.std(ddof=1)
print(f"{'구간 수':>7s} {'χ²':>9s} {'df':>4s} {'p':>9s}")
for k in [4, 6, 8, 12, 20]:
    qs = np.linspace(0, 1, k + 1)[1:-1]
    edges = np.r_[-np.inf, mu + sd * stats.norm.ppf(qs), np.inf]
    obs = np.histogram(x, bins=edges)[0]
    r = stats.chisquare(f_obs=obs, f_exp=np.full(k, len(x) / k), ddof=2)
    print(f"{k:7d} {r.statistic:9.4f} {k - 3:4d} {r.pvalue:9.4f}")
print(f"참고: 샤피로·윌크 p = {stats.shapiro(x).pvalue:.6f}")

   구간 수        χ²   df         p
      4    8.6667    1    0.0032
      6   10.0800    3    0.0179
      8   16.3200    5    0.0060
     12   17.3600    9    0.0434
     20   27.2000   17    0.0552
참고: 샤피로·윌크 p = 0.000007

from scipy import stats

obs = [90, 70, 40]
print(f"f_exp 생략     → p = {stats.chisquare(f_obs=obs).pvalue:.6f}")
print(f"f_exp 명시     → p = "
      f"{stats.chisquare(f_obs=obs, f_exp=[100, 60, 40]).pvalue:.6f}")

f_exp 생략     → p = 0.000075
f_exp 명시     → p = 0.263597

def run_gof(obs, ratios, n_estimated=0, label=""):
    """의도를 코드에 드러내는 적합도 검정 래퍼."""
    import numpy as np
    obs = np.asarray(obs, float)
    p = np.asarray(ratios, float)
    p = p / p.sum()
    exp = obs.sum() * p
    r = stats.chisquare(f_obs=obs, f_exp=exp, ddof=n_estimated)
    df = len(obs) - 1 - n_estimated
    print(f"{label}: χ²={r.statistic:.4f}, df={df}, p={r.pvalue:.4f}, "
          f"E_min={exp.min():.2f}")
    if exp.min() < 5:
        print("   ⚠ 기대도수가 5 미만인 칸이 있다. 병합이나 정확검정을 고려하라.")
    return r

run_gof([90, 70, 40], [0.5, 0.3, 0.2], label="주머니 구슬")
run_gof([4, 13, 7], [1, 1, 1], label="가위바위보")
run_gof([30, 12, 5, 2], [0.6, 0.25, 0.1, 0.05], label="소표본 예")

주머니 구슬: χ²=2.6667, df=2, p=0.2636, E_min=40.00
가위바위보: χ²=5.2500, df=2, p=0.0724, E_min=8.00
소표본 예: χ²=0.1020, df=3, p=0.9916, E_min=2.45
   ⚠ 기대도수가 5 미만인 칸이 있다. 병합이나 정확검정을 고려하라.