독립성 검정¶

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

def compute_expected(observed_counts):
    """관측도수로부터 독립 아래의 기대도수를 계산한다.

    독립의 정의 P(A and B) = P(A)P(B)를 그대로 옮긴 것이다.
    주변 확률의 곱으로 결합 확률을 만들고 전체 개수를 곱한다.
    """
    row_totals = observed_counts.sum(axis=1)
    col_totals = observed_counts.sum(axis=0)

    # reshape로 (r,1)과 (1,c)를 만들면 브로드캐스팅으로 (r,c) 곱이 나온다.
    row_pmf = row_totals.reshape((-1, 1)) / row_totals.sum()
    col_pmf = col_totals.reshape((1, -1)) / col_totals.sum()

    joint_pmf = row_pmf * col_pmf
    expected_counts = joint_pmf * row_totals.sum()

    return expected_counts

# 분할표 형태의 관측도수
observed_counts = np.array([[934, 1070], [113, 92], [20, 8]])
expected_counts = compute_expected(observed_counts)

# 자유도 = (행-1)(열-1)
degrees_of_freedom = (observed_counts.shape[0] - 1) * (observed_counts.shape[1] - 1)

# 검정통계량 계산
chi_squared_statistic = np.sum((observed_counts - expected_counts) ** 2 / expected_counts)
p_value = stats.chi2(degrees_of_freedom).sf(chi_squared_statistic)

# 통계량과 p-값 출력
print(f"chi_squared_statistic = {chi_squared_statistic:.02f}")
print(f"p_value = {p_value:.02%}")

# 카이제곱 분포를 그리고 관측된 통계량 자리를 표시한다
fig, ax = plt.subplots(figsize=(12, 4))

x_values = np.linspace(0, chi_squared_statistic, 100)
y_values = stats.chi2(degrees_of_freedom).pdf(x_values)
ax.plot(x_values, y_values, color='b', linewidth=3)

x_fill_left = np.concatenate([[0], x_values, [chi_squared_statistic], [0]])
y_fill_left = np.concatenate([[0], y_values, [0], [0]])
ax.fill(x_fill_left, y_fill_left, color='b', alpha=0.1)

x_values_right = np.linspace(chi_squared_statistic, 20, 100)
y_values_right = stats.chi2(degrees_of_freedom).pdf(x_values_right)
ax.plot(x_values_right, y_values_right, color='r', linewidth=3)

x_fill_right = np.concatenate([[chi_squared_statistic], x_values_right, [20], [chi_squared_statistic]])
y_fill_right = np.concatenate([[0], y_values_right, [0], [0]])
ax.fill(x_fill_right, y_fill_right, color='r', alpha=0.1)

ax.annotate(f'p_value = {p_value:.02%}', xy=(12.5, 0.01), xytext=(16.5, 0.10),
            fontsize=15, arrowprops=dict(color='k', 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.show()

chi_squared_statistic = 11.81
p_value = 0.27%

from fractions import Fraction as F

O = [[934, 1070], [113, 92], [20, 8]]
R = [sum(row) for row in O]
C = [sum(O[i][j] for i in range(3)) for j in range(2)]
n = sum(R)

# 정의대로 — 기대도수를 유리수로 두면 chi^2 도 유리수가 된다.
chi2_def = sum((O[i][j] - F(R[i] * C[j], n)) ** 2 / F(R[i] * C[j], n)
               for i in range(3) for j in range(2))

# (2) 에서 유도한 r x 2 닫힌 꼴
d = [O[i][0] - F(R[i] * C[0], n) for i in range(3)]
chi2_closed = F(n * n, C[0] * C[1]) * sum(di ** 2 / R[i] for i, di in enumerate(d))

print(f"정의대로 = {float(chi2_def):.10f}")
print(f"닫힌 꼴  = {float(chi2_closed):.10f}")
print(f"유리수로 똑같은가: {chi2_def == chi2_closed}")
print(f"d 의 합 = {sum(d)}  (열 합 보존)")

정의대로 = 11.8061346670
닫힌 꼴  = 11.8061346670
유리수로 똑같은가: True
d 의 합 = 0  (열 합 보존)

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

# 분할표 형태의 관측도수
observed_counts = np.array([[934, 1070], [113, 92], [20, 8]])

# chi2_contingency는 2x2 표에 한해 Yates 연속성 보정을 **기본으로 적용한다**.
# 여기는 3x2라 보정이 없으므로 위의 수동 계산과 정확히 같은 값이 나온다.
chi_squared_statistic, p_value, degrees_of_freedom, expected_counts = stats.chi2_contingency(observed_counts)

# 통계량과 p-값 출력
print(f"chi_squared_statistic = {chi_squared_statistic:.02f}")
print(f"p_value = {p_value:.02%}", end="\n\n")

print("expected_counts")
print(expected_counts, end="\n\n")

# 그림으로 확인한다
fig, ax = plt.subplots(figsize=(12, 4))

x_values = np.linspace(0, chi_squared_statistic, 100)
y_values = stats.chi2(degrees_of_freedom).pdf(x_values)
ax.plot(x_values, y_values, color='b', linewidth=3)

x_fill_left = np.concatenate([[0], x_values, [chi_squared_statistic], [0]])
y_fill_left = np.concatenate([[0], y_values, [0], [0]])
ax.fill(x_fill_left, y_fill_left, color='b', alpha=0.1)

x_values_right = np.linspace(chi_squared_statistic, 20, 100)
y_values_right = stats.chi2(degrees_of_freedom).pdf(x_values_right)
ax.plot(x_values_right, y_values_right, color='r', linewidth=3)

x_fill_right = np.concatenate([[chi_squared_statistic], x_values_right, [20], [chi_squared_statistic]])
y_fill_right = np.concatenate([[0], y_values_right, [0], [0]])
ax.fill(x_fill_right, y_fill_right, color='r', alpha=0.1)

ax.annotate(f'p_value = {p_value:.02%}', xy=(12.5, 0.01), xytext=(16.5, 0.10),
            fontsize=15, arrowprops=dict(color='k', 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.show()

chi_squared_statistic = 11.81
p_value = 0.27%

expected_counts
[[ 955.86410371 1048.13589629]
 [  97.78050961  107.21949039]
 [  13.35538668   14.64461332]]

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

def main():
    observed = np.array([[11, 3, 8], [2, 9, 14], [12, 13, 28]])

    statistic, p_value, df, expected = stats.chi2_contingency(observed)
    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 = (15.0, 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 = 11.94
p_value   = 1.78%

"""
Expected Frequency Calculator for Contingency Table
Computes expected frequencies under independence assumption
"""

import numpy as np
import pandas as pd
from scipy import stats

# 관측도수
observed = np.array([
    [10, 20, 30, 40, 20],
    [5, 15, 40, 50, 10],
    [10, 10, 20, 30, 19]
])

alpha = 0.05

# 행합과 열합, 곧 주변도수
row_totals = observed.sum(axis=1)
col_totals = observed.sum(axis=0)
grand_total = observed.sum()

print("=" * 70)
print("OBSERVED FREQUENCIES")
print("=" * 70)
obs_df = pd.DataFrame(observed,
                      index=['Row 1 (20)', 'Row 2 (30)', 'Row 3 (40)'],
                      columns=['Col 1', 'Col 2', 'Col 3', 'Col 4', 'Col 5'])
obs_df['Row Total'] = row_totals
print(obs_df)
print(f"\nColumn Totals: {col_totals}")
print(f"Grand Total: {grand_total}")

print("\n" + "=" * 70)
print("EXPECTED FREQUENCIES")
print("=" * 70)
print("Formula: E_ij = (Row_i_total × Column_j_total) / Grand_total\n")

# 주변도수의 곱으로 기대도수를 만든다
expected = np.zeros_like(observed, dtype=float)
for i in range(observed.shape[0]):
    for j in range(observed.shape[1]):
        expected[i, j] = (row_totals[i] * col_totals[j]) / grand_total

# 기대도수 출력
exp_df = pd.DataFrame(expected,
                      index=['Row 1 (20)', 'Row 2 (30)', 'Row 3 (40)'],
                      columns=['Col 1', 'Col 2', 'Col 3', 'Col 4', 'Col 5'])
exp_df['Row Total'] = exp_df.sum(axis=1)
print(exp_df)
print(f"\nColumn Totals: {expected.sum(axis=0)}")
print(f"Grand Total: {expected.sum()}")

print("\n" + "=" * 70)
print("DETAILED EXPECTED FREQUENCY CALCULATIONS")
print("=" * 70)
for i in range(observed.shape[0]):
    print(f"\nRow {i+1} (Row Total = {row_totals[i]}):")
    for j in range(observed.shape[1]):
        calculation = f"E[{i+1},{j+1}] = ({row_totals[i]} × {col_totals[j]}) / {grand_total}"
        result = f"= {row_totals[i] * col_totals[j]} / {grand_total} = {expected[i,j]:.4f}"
        print(f"  {calculation} {result}")

print("\n" + "=" * 70)
print("CHI-SQUARE CONTRIBUTIONS")
print("=" * 70)
print("Formula: (Observed - Expected)² / Expected\n")

chi_sq_contrib = (observed - expected)**2 / expected
chi_sq_df = pd.DataFrame(chi_sq_contrib,
                         index=['Row 1 (20)', 'Row 2 (30)', 'Row 3 (40)'],
                         columns=['Col 1', 'Col 2', 'Col 3', 'Col 4', 'Col 5'])
print(chi_sq_df)
print(f"\nChi-square statistic: {chi_sq_contrib.sum():.4f}")
print(f"Degrees of freedom: {(observed.shape[0]-1) * (observed.shape[1]-1)}")

p_value = stats.chi2(df=(observed.shape[0]-1) * (observed.shape[1]-1)).sf(chi_sq_contrib.sum())
print(f"{p_value = :.4f}")
if p_value < alpha:
    print("We have enough evidence to reject the null hypothesis that X and Y are independent.")
else:
    print("We do not have enough evidence to reject the null hypothesis that X and Y are independent.")

======================================================================
OBSERVED FREQUENCIES
======================================================================
            Col 1  Col 2  Col 3  Col 4  Col 5  Row Total
Row 1 (20)     10     20     30     40     20        120
Row 2 (30)      5     15     40     50     10        120
Row 3 (40)     10     10     20     30     19         89

Column Totals: [ 25  45  90 120  49]
Grand Total: 329

======================================================================
EXPECTED FREQUENCIES
======================================================================
Formula: E_ij = (Row_i_total × Column_j_total) / Grand_total

               Col 1      Col 2      Col 3      Col 4      Col 5  Row Total
Row 1 (20)  9.118541  16.413374  32.826748  43.768997  17.872340      120.0
Row 2 (30)  9.118541  16.413374  32.826748  43.768997  17.872340      120.0
Row 3 (40)  6.762918  12.173252  24.346505  32.462006  13.255319       89.0

Column Totals: [ 25.  45.  90. 120.  49.]
Grand Total: 328.99999999999994

======================================================================
DETAILED EXPECTED FREQUENCY CALCULATIONS
======================================================================

Row 1 (Row Total = 120):
  E[1,1] = (120 × 25) / 329 = 3000 / 329 = 9.1185
  E[1,2] = (120 × 45) / 329 = 5400 / 329 = 16.4134
  E[1,3] = (120 × 90) / 329 = 10800 / 329 = 32.8267
  E[1,4] = (120 × 120) / 329 = 14400 / 329 = 43.7690
  E[1,5] = (120 × 49) / 329 = 5880 / 329 = 17.8723

Row 2 (Row Total = 120):
  E[2,1] = (120 × 25) / 329 = 3000 / 329 = 9.1185
  E[2,2] = (120 × 45) / 329 = 5400 / 329 = 16.4134
  E[2,3] = (120 × 90) / 329 = 10800 / 329 = 32.8267
  E[2,4] = (120 × 120) / 329 = 14400 / 329 = 43.7690
  E[2,5] = (120 × 49) / 329 = 5880 / 329 = 17.8723

Row 3 (Row Total = 89):
  E[3,1] = (89 × 25) / 329 = 2225 / 329 = 6.7629
  E[3,2] = (89 × 45) / 329 = 4005 / 329 = 12.1733
  E[3,3] = (89 × 90) / 329 = 8010 / 329 = 24.3465
  E[3,4] = (89 × 120) / 329 = 10680 / 329 = 32.4620
  E[3,5] = (89 × 49) / 329 = 4361 / 329 = 13.2553

======================================================================
CHI-SQUARE CONTRIBUTIONS
======================================================================
Formula: (Observed - Expected)² / Expected

               Col 1     Col 2     Col 3     Col 4     Col 5
Row 1 (20)  0.085208  0.783744  0.243414  0.324553  0.253293
Row 2 (30)  1.860208  0.121707  1.567488  0.887053  3.467579
Row 3 (40)  1.549435  0.387984  0.775968  0.186725  2.489669

Chi-square statistic: 14.9840
Degrees of freedom: 8
p_value = 0.0595
We do not have enough evidence to reject the null hypothesis that X and Y are independent.

import numpy as np
import pandas as pd
import random
from scipy import stats

# 제목 세 가지의 클릭 수
headlines = pd.DataFrame({
    'Click': [14, 8, 12],
    'No-click': [986, 992, 988],
    'Headline': ['Headline A', 'Headline B', 'Headline C']
})

# 결과를 행, 헤드라인을 열로 놓은 분할표를 만든다.
click_rate = headlines.copy()
clicks = click_rate.set_index('Headline')[['Click', 'No-click']].T

print("Observed Contingency Table:")
print(clicks)
print(f"\nTotal: {clicks.values.sum()}")

Observed Contingency Table:
Headline  Headline A  Headline B  Headline C
Click             14           8          12
No-click         986         992         988

Total: 3000

def chi2_stat(observed, expected):
    """카이제곱 통계량을 구한다. 칸마다 (관측-기대)^2/기대 를 더한 값이다."""
    pearson_residuals = []
    for row, expect in zip(observed, expected):
        pearson_residuals.append([(observe - expect) ** 2 / expect
                                  for observe in row])
    return np.sum(pearson_residuals)

# 관측된 카이제곱 통계량
row_average = clicks.mean(axis=1)
expected = np.array([[row_average['Click'], row_average['Click'], row_average['Click']],
                     [row_average['No-click'], row_average['No-click'], row_average['No-click']]])

chi2_obs = chi2_stat(clicks.values, row_average.values)
print(f"Observed chi-square: {chi2_obs:.4f}")

# 재표본추출 방법 — 분포를 가정하지 않는다
def perm_fun_chisq(box):
    """
    Generate permuted contingency table by random allocation.

    Parameters:
    -----------
    box : list
        Binary response (1 = click, 0 = no-click) for all users

    Returns:
    --------
    float : Chi-square statistic for permuted table
    """
    random.shuffle(box)
    # 앞 1000개를 A, 다음 1000개를 B, 마지막 1000개를 C에 배정한다
    sample_clicks = [sum(box[0:1000]),
                     sum(box[1000:2000]),
                     sum(box[2000:3000])]
    sample_noclicks = [1000 - n for n in sample_clicks]
    return chi2_stat([sample_clicks, sample_noclicks], row_average.values)

# 상자를 만든다. 클릭은 1, 비클릭은 0 이다.
box = [1] * 34 + [0] * 2966

# 순열검정 실행
random.seed(42)
perm_chi2 = [perm_fun_chisq(box) for _ in range(2000)]

p_value_resamp = sum(np.array(perm_chi2) >= chi2_obs) / len(perm_chi2)
print(f"Resampling p-value: {p_value_resamp:.4f}")

Observed chi-square: 1.6659
Resampling p-value: 0.4750

import numpy as np

n, c, m = 3000, 3, 1000
print(f"순열 귀무분포  평균 = {np.mean(perm_chi2):.4f},  표준편차 = {np.std(perm_chi2):.4f}")
print(f"이론  n(c-1)/(n-1) = {n * (c - 1) / (n - 1):.4f},  chi2_2 의 평균 = 2")
print(f"몬테카를로 오차 = {np.std(perm_chi2) / np.sqrt(len(perm_chi2)):.4f}")

순열 귀무분포  평균 = 1.9789,  표준편차 = 1.9640
이론  n(c-1)/(n-1) = 2.0007,  chi2_2 의 평균 = 2
몬테카를로 오차 = 0.0439

def sample_with_replacement(box):
    """
    Generate contingency table by sampling with replacement.
    """
    sample_clicks = [sum(random.choices(box, k=1000)),
                     sum(random.choices(box, k=1000)),
                     sum(random.choices(box, k=1000))]
    sample_noclicks = [1000 - n for n in sample_clicks]
    return chi2_stat([sample_clicks, sample_noclicks], row_average.values)

# 복원추출 방식으로 실행
random.seed(42)
perm_chi2_wr = [sample_with_replacement(box) for _ in range(2000)]

p_value_wr = sum(np.array(perm_chi2_wr) >= chi2_obs) / len(perm_chi2_wr)
print(f"Resampling (with replacement) p-value: {p_value_wr:.4f}")

Resampling (with replacement) p-value: 0.6745

print(f"복원 귀무분포  평균 = {np.mean(perm_chi2_wr):.4f},  표준편차 = {np.std(perm_chi2_wr):.4f}")
print("이론  E[chi2] = c = 3")
print(f"몬테카를로 오차 = {np.std(perm_chi2_wr) / np.sqrt(len(perm_chi2_wr)):.4f}")

복원 귀무분포  평균 = 3.0976,  표준편차 = 2.4859
이론  E[chi2] = c = 3
몬테카를로 오차 = 0.0556

# 모수적 카이제곱 검정
chi2_param, p_param, df, expected_param = stats.chi2_contingency(clicks.values)

print(f"\nComparison:")
print(f"Parametric chi-square: {chi2_param:.4f}, p-value: {p_param:.4f}")
print(f"Resampling (without repl): p-value: {p_value_resamp:.4f}")
print(f"Resampling (with repl): p-value: {p_value_wr:.4f}")

Comparison:
Parametric chi-square: 1.6659, p-value: 0.4348
Resampling (without repl): p-value: 0.4750
Resampling (with repl): p-value: 0.6745

from fractions import Fraction as F
from math import comb

n, m, c, k = 3000, 1000, 3, 34
const = F(c * n, k * (n - k))      # 보기 6 (1) 의 상수 2250/25211
obs = F(42000, 25211)              # 관측 통계량, 유리수로
denom = comb(n, k)

total = F(0)   # 확률의 총합 — 1 이어야 한다
tail = F(0)    # p-값
mean = F(0)    # 분포의 평균
for k1 in range(k + 1):
    for k2 in range(k + 1 - k1):
        k3 = k - k1 - k2
        prob = F(comb(m, k1) * comb(m, k2) * comb(m, k3), denom)
        stat = const * sum((F(a) - F(k, c)) ** 2 for a in (k1, k2, k3))
        total += prob
        mean += prob * stat
        if stat >= obs:
            tail += prob

se = (p_value_resamp * (1 - p_value_resamp) / 2000) ** 0.5
print(f"확률의 총합      = {total}   (유리수로 정확히)")
print(f"정확한 순열 p-값 = {float(tail):.7f}")
print(f"모의 순열 p-값   = {p_value_resamp:.7f}   (2,000 회)")
print(f"몬테카를로 오차  = {se:.4f}")
print(f"차이 / 오차      = {abs(float(tail) - p_value_resamp) / se:.2f}")
print(f"정확한 순열 평균 = {float(mean):.7f}   이론 n(c-1)/(n-1) = {n * (c - 1) / (n - 1):.7f}")

확률의 총합      = 1   (유리수로 정확히)
정확한 순열 p-값 = 0.4824141
모의 순열 p-값   = 0.4750000   (2,000 회)
몬테카를로 오차  = 0.0112
차이 / 오차      = 0.66
정확한 순열 평균 = 2.0006669   이론 n(c-1)/(n-1) = 2.0006669

import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 비복원 재표본 분포 — 순열검정에 해당한다
ax1.hist(perm_chi2, bins=40, alpha=0.7, color='steelblue', edgecolor='black')
ax1.axvline(chi2_obs, color='red', linewidth=2, label=f'Observed = {chi2_obs:.2f}')
ax1.set_xlabel('Chi-square Statistic')
ax1.set_ylabel('Frequency')
ax1.set_title(f'Resampling Distribution (without replacement)\np-value = {p_value_resamp:.4f}')
ax1.legend()
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)

# 복원 재표본 분포 — 붓스트랩에 해당한다
ax2.hist(perm_chi2_wr, bins=40, alpha=0.7, color='forestgreen', edgecolor='black')
ax2.axvline(chi2_obs, color='red', linewidth=2, label=f'Observed = {chi2_obs:.2f}')
ax2.set_xlabel('Chi-square Statistic')
ax2.set_ylabel('Frequency')
ax2.set_title(f'Resampling Distribution (with replacement)\np-value = {p_value_wr:.4f}')
ax2.legend()
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)

plt.tight_layout()
plt.show()

import numpy as np

for name, d, theory in [("비복원", perm_chi2, 3000 * 2 / 2999),
                        ("복원", perm_chi2_wr, 3.0)]:
    a = np.array(d)
    print(f"{name:>4}  평균 {a.mean():.4f} (이론 {theory:.4f})  "
          f"표준편차 {a.std():.4f}  중앙값 {np.median(a):.4f}  "
          f"p-값 {(a >= chi2_obs).mean():.4f}")

 비복원  평균 1.9789 (이론 2.0007)  표준편차 1.9640  중앙값 1.4874  p-값 0.4750
  복원  평균 3.0976 (이론 3.0000)  표준편차 2.4859  중앙값 2.4692  p-값 0.6745

import numpy as np
from scipy import stats

def report(T, title):
    print(f"[{title}]")
    marginal = T.sum(2)                      # Z 를 무시한 X×Y 표
    c, p, df, _ = stats.chi2_contingency(marginal, correction=False)
    print(f"  주변표 (Z 무시)\n{marginal.astype(int)}")
    print(f"    χ² = {c:.4f},  df = {df},  p = {p:.4g}")
    for z in range(T.shape[2]):
        sub = T[:, :, z]
        if sub.sum(0).min() > 0 and sub.sum(1).min() > 0:
            c2, p2, _, _ = stats.chi2_contingency(sub, correction=False)
            print(f"  Z={z} 층\n{sub.astype(int)}    "
                  f"χ² = {c2:.4f},  p = {p2:.4g}")
        else:
            print(f"  Z={z} 층\n{sub.astype(int)}    (한 칸이 0: 완전 종속)")
    print()

n = 800

# ① X 와 Y 는 주변적으로 독립이지만 Z 를 고정하면 완전히 종속
P1 = np.zeros((2, 2, 2))
P1[0, 0, 0] = P1[0, 1, 1] = P1[1, 0, 1] = P1[1, 1, 0] = 0.25
report(np.round(P1 * n), "① 주변 독립, 조건부 종속")

# ② Z 가 X 와 Y 의 공통 원인 — 주변적으로 종속, 조건부로 독립
P2 = np.zeros((2, 2, 2))
for z, (px, py) in enumerate([(0.2, 0.2), (0.8, 0.8)]):
    for x in range(2):
        for y in range(2):
            P2[x, y, z] = 0.5 * (px if x else 1 - px) * (py if y else 1 - py)
report(np.round(P2 * n), "② 주변 종속, 조건부 독립")

[① 주변 독립, 조건부 종속]
  주변표 (Z 무시)
[[200 200]
 [200 200]]
    χ² = 0.0000,  df = 1,  p = 1
  Z=0 층
[[200   0]
 [  0 200]]    χ² = 400.0000,  p = 5.507e-89
  Z=1 층
[[  0 200]
 [200   0]]    χ² = 400.0000,  p = 5.507e-89

[② 주변 종속, 조건부 독립]
  주변표 (Z 무시)
[[272 128]
 [128 272]]
    χ² = 103.6800,  df = 1,  p = 2.378e-24
  Z=0 층
[[256  64]
 [ 64  16]]    χ² = 0.0000,  p = 1
  Z=1 층
[[ 16  64]
 [ 64 256]]    χ² = 0.0000,  p = 1

① 충돌자 구조             ② 공통 원인 구조
   X → Z ← Y                  X ← Z → Y

Z 를 무시: X ⊥ Y            Z 를 무시: X 와 Y 가 연관
Z 로 층화: X 와 Y 가 연관    Z 로 층화: X ⊥ Y
            (선택 편향!)                (교란 제거)

import numpy as np
from scipy import stats

def linear_by_linear(T, u=None, v=None):
    """순서형 × 순서형 표의 선형-선형 연관 검정."""
    T = np.asarray(T, float)
    r, c = T.shape
    n = T.sum()
    u = np.arange(r, dtype=float) if u is None else np.asarray(u, float)
    v = np.arange(c, dtype=float) if v is None else np.asarray(v, float)
    ru, rv, w = np.repeat(u, c), np.tile(v, r), T.ravel()
    mu, mv = (w * ru).sum() / n, (w * rv).sum() / n
    su = np.sqrt((w * (ru - mu)**2).sum() / n)
    sv = np.sqrt((w * (rv - mv)**2).sum() / n)
    rho = (w * (ru - mu) * (rv - mv)).sum() / n / (su * sv)
    m2 = (n - 1) * rho**2
    return rho, m2, stats.chi2.sf(m2, 1)

tables = {
    "단조 연관": np.array([[30, 20, 10,  5],
                          [20, 30, 20, 10],
                          [10, 20, 30, 20],
                          [ 5, 10, 20, 30]], float),
    "대칭 U자": np.array([[40, 10, 10, 40],
                         [10, 40, 40, 10],
                         [10, 40, 40, 10],
                         [40, 10, 10, 40]], float),
}
for name, T in tables.items():
    c, p, df, _ = stats.chi2_contingency(T, correction=False)
    rho, m2, pm = linear_by_linear(T)
    print(f"[{name}]  n = {T.sum():.0f}")
    print(f"  일반 독립성   χ² = {c:8.4f},  df = {df},  p = {p:.4g}")
    print(f"  선형-선형     r = {rho:+.6f},  M² = {m2:8.4f},  df = 1,  "
          f"p = {pm:.4g}")

[단조 연관]  n = 290
  일반 독립성   χ² =  66.3332,  df = 9,  p = 7.935e-11
  선형-선형     r = +0.443609,  M² =  56.8720,  df = 1,  p = 4.651e-14
[대칭 U자]  n = 400
  일반 독립성   χ² = 144.0000,  df = 9,  p = 1.539e-26
  선형-선형     r = +0.000000,  M² =   0.0000,  df = 1,  p = 1

[[40 10 10 40]
 [10 40 40 10]
 [10 40 40 10]
 [40 10 10 40]]

import numpy as np
from scipy import stats

def g2_and_mi(T):
    T = np.asarray(T, float)
    n = T.sum()
    P = T / n
    exp = np.outer(T.sum(1), T.sum(0)) / n
    mask = T > 0
    g2 = 2 * np.sum(T[mask] * np.log(T[mask] / exp[mask]))
    indep = np.outer(P.sum(1), P.sum(0))
    mi = np.sum(P[mask] * np.log(P[mask] / indep[mask]))
    return g2, mi

tables = {
    "결제수단 × 요일": np.array([[30, 50, 20], [40, 60, 30]], float),
    "강한 연관 2×2": np.array([[40, 10], [10, 40]], float),
    "성별 × 손잡이": np.array([[934, 1070], [113, 92], [20, 8]], float),
}
for name, T in tables.items():
    g2, mi = g2_and_mi(T)
    chi2, p, df, _ = stats.chi2_contingency(T, correction=False)
    n = T.sum()
    print(f"[{name}]  n = {n:.0f},  df = {df}")
    print(f"  χ² = {chi2:8.4f}   G² = {g2:8.4f}   2n·MI = {2 * n * mi:8.4f}")
    print(f"  MI = {mi:.6f} nat = {mi / np.log(2):.6f} bit"
          f"    G² 의 p = {stats.chi2.sf(g2, df):.4f}")

[결제수단 × 요일]  n = 230,  df = 2
  χ² =   0.4320   G² =   0.4332   2n·MI =   0.4332
  MI = 0.000942 nat = 0.001358 bit    G² 의 p = 0.8053
[강한 연관 2×2]  n = 100,  df = 1
  χ² =  36.0000   G² =  38.5490   2n·MI =  38.5490
  MI = 0.192745 nat = 0.278072 bit    G² 의 p = 0.0000
[성별 × 손잡이]  n = 2237,  df = 2
  χ² =  11.8061   G² =  11.9606   2n·MI =  11.9606
  MI = 0.002673 nat = 0.003857 bit    G² 의 p = 0.0025

두 범주형 변수의 연관을 본다
    │
    ├─ 관측이 독립인가 (군집·반복측정 아님)
    │      └─ 아니오 → McNemar / 코크런 Q / 혼합모형
    │
    ├─ 변수의 척도는
    │      ├─ 둘 다 명목형 ──────→ 일반 χ² 독립성 검정
    │      ├─ 한쪽만 순서형 ─────→ 코크런·아미티지 추세검정
    │      └─ 둘 다 순서형 ──────→ 선형-선형 연관 검정 (연습문제 8)
    │
    ├─ 기대도수를 확인
    │      └─ 작으면 → 피셔 / 바너드 / 몬테카를로
    │
    ├─ 셋째 변수가 있는가
    │      ├─ 교란변수 → 층화 (맨텔·헨젤) / 로그선형모형
    │      └─ 충돌자  → 층화하면 안 됨 (연습문제 7)
    │
    └─ 언제나: 효과크기 + 잔차 + 원 도수표