독립성 검정 (수동 계산과 그림)¶

import numpy as np
from scipy import stats

def compute_expected(observed_counts: np.ndarray) -> np.ndarray:
    """독립이라는 가정 아래의 기대도수. E_ij = (행합 x 열합) / 전체."""
    # 이 공식이 곧 독립의 정의다. P(A and B) = P(A)P(B)의 양변에 n을 곱하면
    # n * (행합/n) * (열합/n) = 행합 * 열합 / n 이 된다.
    row_totals = observed_counts.sum(axis=1, keepdims=True)
    col_totals = observed_counts.sum(axis=0, keepdims=True)
    total = observed_counts.sum()
    return (row_totals @ col_totals) / total


# 확인: 기대도수의 행합과 열합은 관측도수의 것과 정확히 같아야 한다.
demo = np.array([[934., 1070.], [113., 92.], [20., 8.]])
E = compute_expected(demo)
print(np.round(E, 2))
print("행합 일치:", np.allclose(E.sum(axis=1), demo.sum(axis=1)))
print("열합 일치:", np.allclose(E.sum(axis=0), demo.sum(axis=0)))
print("rank(E) =", np.linalg.matrix_rank(E))

[[ 955.86 1048.14]
 [  97.78  107.22]
 [  13.36   14.64]]
행합 일치: True
열합 일치: True
rank(E) = 1

# 행이 한 변수의 수준, 열이 다른 변수의 수준이다.
# 실수로 나눗셈을 하게 되므로 dtype=float 로 만들어 둔다.
observed_counts = np.array([[934, 1070],
                            [113,   92],
                            [ 20,    8]], dtype=float)

expected_counts = compute_expected(observed_counts)
df = (observed_counts.shape[0] - 1) * (observed_counts.shape[1] - 1)

chi2 = np.sum((observed_counts - expected_counts)**2 / expected_counts)
p_value = stats.chi2(df).sf(chi2)      # 언제나 우측검정

print(f"chi_squared_statistic = {chi2:.2f}")
print(f"p_value = {p_value:.2%}")

# 칸별 기여와, 자유도 2 의 닫힌 꼴
print(f"칸별 기여\n{np.round((observed_counts - expected_counts)**2 / expected_counts, 4)}")
print(f"exp(-chi2/2) = {np.exp(-chi2 / 2):.7f}   sf = {p_value:.7f}")
print(f"크래머 V = {np.sqrt(chi2 / observed_counts.sum()):.4f}")

# 같은 표를 scipy 에 넣는다
c2, p2, df2, _ = stats.chi2_contingency(observed_counts)
print(f"chi2_contingency: chi2 = {c2:.4f}, df = {df2}, p = {p2:.7f}")

chi_squared_statistic = 11.81
p_value = 0.27%
칸별 기여
[[0.5001 0.4561]
 [2.3689 2.1604]
 [3.3058 3.0148]]
exp(-chi2/2) = 0.0027311   sf = 0.0027311
크래머 V = 0.0726
chi2_contingency: chi2 = 11.8061, df = 2, p = 0.0027311

import matplotlib.pyplot as plt

# 카이제곱 분포를 그리고 통계량을 경계로 좌우를 나눠 칠한다.
# 오른쪽 넓이가 p-값이다. 이 검정은 방향이 없어 언제나 우측만 본다.
fig, ax = plt.subplots(figsize=(12, 4))

x_left = np.linspace(0, chi2, 200)
y_left = stats.chi2(df).pdf(x_left)
ax.plot(x_left, y_left, linewidth=3)
x_fill = np.concatenate([[0], x_left, [chi2], [0]])
y_fill = np.concatenate([[0], y_left, [0], [0]])
ax.fill(x_fill, y_fill, alpha=0.1)

x_right = np.linspace(chi2, max(20, chi2 + 5), 200)
y_right = stats.chi2(df).pdf(x_right)
ax.plot(x_right, y_right, linewidth=3)
x_fill_r = np.concatenate([[chi2], x_right, [max(20, chi2 + 5)], [chi2]])
y_fill_r = np.concatenate([[0], y_right, [0], [0]])
ax.fill(x_fill_r, y_fill_r, alpha=0.1)

ax.annotate(f"p_value = {p_value:.2%}",
            xy=(chi2 * 0.8, y_left.max() * 0.15),
            xytext=(chi2 * 0.9 + 5, y_left.max() * 0.6),
            fontsize=12, 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([[934., 1070.], [113., 92.], [20., 8.]])
exp = np.outer(obs.sum(1), obs.sum(0)) / obs.sum()
dev = obs - exp

print("O − E 의 행 합", np.round(dev.sum(1), 10).tolist())
print("O − E 의 열 합", np.round(dev.sum(0), 10).tolist())

resid = dev / np.sqrt(exp)
print(f"\n표준화 잔차\n{np.round(resid, 4)}")
print(f"R 의 행 합 {np.round(resid.sum(1), 4).tolist()}  ← 0 이 아니다")

O − E 의 행 합 [0.0, 0.0, 0.0]
O − E 의 열 합 [0.0, 0.0]

표준화 잔차
[[-0.7072  0.6753]
 [ 1.5391 -1.4698]
 [ 1.8182 -1.7363]]
R 의 행 합 [-0.0318, 0.0693, 0.0819]  ← 0 이 아니다

free = np.array([[5.0, 0.0], [0.0, 0.0], [0.0, 0.0]])   # (0,0) 칸만 +5
print("(0,0) 칸을 5 늘리고 주변합을 유지하려면")
adj = free.copy()
adj[0, 1] = -5      # 같은 행에서 빼고
adj[2, 0] = -5      # 같은 열에서 빼고
adj[2, 1] = +5      # 그 교차점에서 다시 더한다
print(adj.astype(int))
print(f"  행 합 {adj.sum(1).astype(int).tolist()},  "
      f"열 합 {adj.sum(0).astype(int).tolist()}   모두 0")

(0,0) 칸을 5 늘리고 주변합을 유지하려면
[[ 5 -5]
 [ 0  0]
 [-5  5]]
  행 합 [0, 0, 0],  열 합 [0, 0]   모두 0

import numpy as np
from scipy import stats

rng = np.random.default_rng(4321)
M = 20_000

print(f"{'표':>7s} {'n':>6s} {'df':>4s} {'통계량 평균':>12s} "
      f"{'통계량 분산':>12s} {'2·df':>6s}")
for (r, c), n in [((3, 2), 2237), ((4, 5), 1000), ((2, 2), 500)]:
    # 주변분포를 임의로 잡되 독립이 참이 되도록 만든다
    P = np.outer(rng.dirichlet(np.ones(r) * 5),
                 rng.dirichlet(np.ones(c) * 5))
    stat = []
    for _ in range(M):
        T = rng.multinomial(n, P.ravel()).reshape(r, c).astype(float)
        if (T.sum(0) > 0).all() and (T.sum(1) > 0).all():
            stat.append(stats.chi2_contingency(T, correction=False)[0])
    stat = np.array(stat)
    df = (r - 1) * (c - 1)
    print(f"{f'{r}×{c}':>7s} {n:6d} {df:4d} {stat.mean():12.4f} "
          f"{stat.var(ddof=1):12.4f} {2 * df:6d}")

      표      n   df       통계량 평균       통계량 분산   2·df
    3×2   2237    2       1.9995       3.9104      4
    4×5   1000   12      12.0108      23.7417     24
    2×2    500    1       0.9893       1.9736      2

P = np.outer([0.4, 0.35, 0.25], [0.55, 0.45])
stat = []
for _ in range(M):
    T = rng.multinomial(2237, P.ravel()).reshape(3, 2).astype(float)
    if (T.sum(0) > 0).all() and (T.sum(1) > 0).all():
        stat.append(stats.chi2_contingency(T, correction=False)[0])
stat = np.array(stat)
print(f"통계량 평균 {stat.mean():.4f}")
for df in [1, 2, 3, 5]:
    print(f"  df={df} 로 읽으면 기각률 "
          f"{np.mean(stat > stats.chi2.ppf(0.95, df)):.4f}")

통계량 평균 1.9783
  df=1 로 읽으면 기각률 0.1442
  df=2 로 읽으면 기각률 0.0500
  df=3 로 읽으면 기각률 0.0197
  df=5 로 읽으면 기각률 0.0029

import numpy as np
from scipy import stats

print("① 도수의 규모가 커지면")
base = np.array([[934., 1070.], [113., 92.], [20., 8.]])
for scale in [1, 1e3, 1e6]:
    c, p, df, _ = stats.chi2_contingency(base * scale, correction=False)
    print(f"  ×{scale:>8.0e}:  χ² = {c:14.4f},  p = {p:.6g}")

print("\n② 0 인 칸이 있으면")
cases = [np.array([[10., 0.], [5., 15.]]),
         np.array([[10., 0.], [0., 15.]]),
         np.array([[10., 5., 0.], [5., 10., 0.]])]
for T in cases:
    try:
        c, p, df, _ = stats.chi2_contingency(T, correction=False)
        print(f"  {T.astype(int).tolist()}:  χ² = {c:.4f}, df = {df}, "
              f"p = {p:.4f}")
    except Exception as e:
        print(f"  {T.astype(int).tolist()}:  {type(e).__name__}")
        print(f"    {e}")

print("\n③ 대수적으로 같은 두 표현")
T = np.array([[1e6 + 3, 1e6], [1e6, 1e6]])
exp = np.outer(T.sum(1), T.sum(0)) / T.sum()
form_a = np.sum((T - exp)**2 / exp)
form_b = np.sum(T**2 / exp) - T.sum()
print(f"  (O−E)²/E 의 합  = {form_a:.12f}")
print(f"  O²/E − n        = {form_b:.12f}")
print(f"  차이 = {abs(form_a - form_b):.3e}")

① 도수의 규모가 커지면
  ×   1e+00:  χ² =        11.8061,  p = 0.00273105
  ×   1e+03:  χ² =     11806.1347,  p = 0
  ×   1e+06:  χ² =  11806134.6670,  p = 0

② 0 인 칸이 있으면
  [[10, 0], [5, 15]]:  χ² = 15.0000, df = 1, p = 0.0001
  [[10, 0], [0, 15]]:  χ² = 25.0000, df = 1, p = 0.0000
  [[10, 5, 0], [5, 10, 0]]:  ValueError
    The internally computed table of expected frequencies has a zero element at (0, 2).

③ 대수적으로 같은 두 표현
  (O−E)²/E 의 합  = 0.000002249995
  O²/E − n        = 0.000002250075
  차이 = 8.040e-11

import numpy as np

obs = np.array([[934., 1070.], [113., 92.], [20., 8.]])
n = obs.sum()

row_keep = obs.sum(axis=1, keepdims=True)
col_keep = obs.sum(axis=0, keepdims=True)
row_flat = obs.sum(axis=1)
col_flat = obs.sum(axis=0)
print(f"keepdims=True 일 때  행합 모양 {row_keep.shape}, "
      f"열합 모양 {col_keep.shape}")
print(f"keepdims 없이       행합 모양 {row_flat.shape}, "
      f"열합 모양 {col_flat.shape}\n")

good = (row_keep @ col_keep) / n
print(f"올바른 기대도수 (모양 {good.shape})\n{np.round(good, 2)}")

try:
    bad = (row_flat @ col_flat) / n
    print(f"\nkeepdims 없이: {bad}   ← 스칼라!  모양 {np.shape(bad)}")
except Exception as e:
    print(f"\nkeepdims 없이: {type(e).__name__}: {e}")

keepdims=True 일 때  행합 모양 (3, 1), 열합 모양 (1, 2)
keepdims 없이       행합 모양 (3,), 열합 모양 (2,)

올바른 기대도수 (모양 (3, 2))
[[ 955.86 1048.14]
 [  97.78  107.22]
 [  13.36   14.64]]

keepdims 없이: ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 3)

square = np.array([[30., 20.], [10., 40.]])
m = square.sum()
r_flat, c_flat = square.sum(axis=1), square.sum(axis=0)
r_keep = square.sum(axis=1, keepdims=True)
c_keep = square.sum(axis=0, keepdims=True)

print(f"올바른 기대도수\n{np.round((r_keep @ c_keep) / m, 4)}")
print(f"\nkeepdims 없이 @ 를 쓰면: {(r_flat @ c_flat) / m:.4f}  ← 내적(스칼라)")
print(f"안전한 대안 np.outer:\n{np.round(np.outer(r_flat, c_flat) / m, 4)}")

올바른 기대도수
[[20. 30.]
 [20. 30.]]

keepdims 없이 @ 를 쓰면: 50.0000  ← 내적(스칼라)
안전한 대안 np.outer:
[[20. 30.]
 [20. 30.]]

r, c = square.sum(axis=1), square.sum(axis=0)
print("① np.outer        :", np.allclose(np.outer(r, c) / m,
                                         (r_keep @ c_keep) / m))
print("② 명시적 reshape  :", np.allclose(r.reshape(-1, 1) * c.reshape(1, -1) / m,
                                         (r_keep @ c_keep) / m))
print("③ np.einsum       :", np.allclose(np.einsum('i,j->ij', r, c) / m,
                                         (r_keep @ c_keep) / m))

① np.outer        : True
② 명시적 reshape  : True
③ np.einsum       : True

import numpy as np
from scipy import stats

def chi2_independence(obs, alpha=0.05, verbose=True):
    """분할표의 독립성 검정을 처음부터 계산한다.

    진단 정보(기대도수 최솟값, 효과크기, 조정 잔차)를 함께 돌려준다.
    """
    obs = np.asarray(obs, float)
    if obs.ndim != 2:
        raise ValueError("2차원 분할표가 필요하다")
    r, c = obs.shape
    n = obs.sum()
    row, col = obs.sum(1), obs.sum(0)
    if row.min() == 0 or col.min() == 0:
        raise ValueError("합이 0 인 행 또는 열이 있다. 먼저 제거하라")

    exp = np.outer(row, col) / n                      # np.outer 가 안전
    assert np.allclose(exp.sum(1), row)               # 검산
    assert np.allclose(exp.sum(0), col)

    chi2 = np.sum((obs - exp)**2 / exp)               # 차이를 먼저
    df = (r - 1) * (c - 1)
    p = stats.chi2(df).sf(chi2)                       # 오른쪽 꼬리
    v = np.sqrt(chi2 / (n * (min(r, c) - 1)))
    adj = (obs - exp) / np.sqrt(exp * np.outer(1 - row / n, 1 - col / n))

    if verbose:
        print(f"χ² = {chi2:.4f},  df = {df},  p = {p:.6f}")
        print(f"크라메르 V = {v:.4f}"
              f"   (독립일 때 기댓값 ≈ "
              f"{np.sqrt(df / (n * (min(r, c) - 1))):.4f})")
        print(f"E_min = {exp.min():.3f}"
              + ("   ⚠ 5 미만" if exp.min() < 5 else ""))
        if p < alpha:
            print(f"조정 잔차 (|·| 최대 {np.abs(adj).max():.3f})\n"
                  f"{np.round(adj, 3)}")
    return {"chi2": chi2, "df": df, "p": p, "V": v,
            "expected": exp, "adj_resid": adj}

_ = chi2_independence([[934, 1070], [113, 92], [20, 8]])
print()
try:
    chi2_independence([[10, 5, 0], [5, 10, 0]])
except ValueError as e:
    print(f"ValueError: {e}")

χ² = 11.8061,  df = 2,  p = 0.002731
크라메르 V = 0.0726   (독립일 때 기댓값 ≈ 0.0299)
E_min = 13.355
조정 잔차 (|·| 최대 3.030)
[[-3.03   3.03 ]
 [ 2.233 -2.233]
 [ 2.53  -2.53 ]]

ValueError: 합이 0 인 행 또는 열이 있다. 먼저 제거하라