모자이크 그림과 도수분포표¶

import pandas as pd

# 네 조합의 개수를 그대로 펼쳐 원자료 형태로 만든다.
#   SUV·사고 28명 / 비SUV·사고 35명 / SUV·무사고 97명 / 비SUV·무사고 104명
data = {'SUV': 28*['예'] + 35*['아니오'] + 97*['예'] + 104*['아니오'],
        '사고': 28*['예'] + 35*['예']    + 97*['아니오'] + 104*['아니오']}
df = pd.DataFrame(data)

# crosstab이 두 범주형 변수를 교차하여 도수를 센다
dg = pd.crosstab(df.SUV, df.사고, rownames=['SUV'], colnames=['사고'])

dg.loc['합계', :] = dg.sum()         # 열 방향 합계를 맨 아래 줄에 추가
dg.loc[:, '합계'] = dg.sum(axis=1)   # 행 방향 합계를 맨 오른쪽 열에 추가
dg = dg.astype(int)                  # 합계를 더하며 실수가 되었으므로 정수로 되돌린다
print(dg)

# (1) 주변합을 고정하고 n11 하나만 바꾸면 나머지 세 칸이 따라온다.
O = pd.crosstab(df.SUV, df.사고).values.astype(int)
r1, r2 = O.sum(1)          # 행합 (아니오, 예)
c1, c2 = O.sum(0)          # 열합 (아니오, 예)
n = O.sum()
print(f"\n행합 {r1}, {r2}   열합 {c1}, {c2}   총합 {n}")

lo = max(0, r1 + c1 - n)
hi = min(r1, c1)
print(f"n11 이 놓일 수 있는 범위 [{lo}, {hi}]  →  주변합이 같은 표가 {hi - lo + 1} 개")

print(f"\n{'n11':>5}{'n12':>5}{'n21':>5}{'n22':>5}{'합':>6}")
for n11 in [lo, 100, int(O[0, 0]), 110, hi]:
    n12, n21 = r1 - n11, c1 - n11
    n22 = n - r1 - c1 + n11
    mark = "   <- 실제 자료" if n11 == O[0, 0] else ""
    print(f"{n11:>5}{n12:>5}{n21:>5}{n22:>5}{n11+n12+n21+n22:>6}{mark}")

사고   아니오   예   합계
SUV
아니오  104  35  139
예     97  28  125
합계   201  63  264

행합 139, 125   열합 201, 63   총합 264
n11 이 놓일 수 있는 범위 [76, 139]  →  주변합이 같은 표가 64 개

  n11  n12  n21  n22     합
   76   63  125    0   264
  100   39  101   24   264
  104   35   97   28   264   <- 실제 자료
  110   29   91   34   264
  139    0   62   63   264

# 모든 칸을 전체 합으로 나누면 상대도수분포표가 된다.
# 모자이크 그림은 바로 이 비율을 넓이로 옮겨 그린 것이다.
dh = dg / dg.loc['합계', '합계']
print(dh)

ct = pd.crosstab(df.SUV, df.사고)          # 합계 줄이 없는 2x2 표
n = ct.values.sum()

print("\n행으로 나눈 표 — P(사고 | SUV)")
print((ct.div(ct.sum(axis=1), axis=0)).round(6))
print("\n열로 나눈 표 — P(SUV | 사고)")
print((ct.div(ct.sum(axis=0), axis=1)).round(6))

# 베이즈 정리로 두 조건부확률을 서로 옮긴다.
p_acc_given_suv = ct.loc['예', '예'] / ct.loc['예'].sum()
p_suv = ct.loc['예'].sum() / n
p_acc = ct['예'].sum() / n
p_suv_given_acc = p_acc_given_suv * p_suv / p_acc

print(f"\nP(사고 | SUV) = {ct.loc['예','예']}/{ct.loc['예'].sum()} = {p_acc_given_suv:.6f}")
print(f"P(SUV)       = {ct.loc['예'].sum()}/{n} = {p_suv:.6f}")
print(f"P(사고)       = {ct['예'].sum()}/{n} = {p_acc:.6f}")
print(f"베이즈:  P(SUV | 사고) = {p_acc_given_suv:.6f} * {p_suv:.6f} / {p_acc:.6f}"
      f" = {p_suv_given_acc:.6f}")
print(f"표에서 직접: {ct.loc['예','예']}/{ct['예'].sum()} = {ct.loc['예','예']/ct['예'].sum():.6f}")
print(f"\n두 조건부확률의 비 = P(SUV)/P(사고) = {p_suv/p_acc:.6f}"
      f"   (실제 비 {p_suv_given_acc/p_acc_given_suv:.6f})")

사고        아니오         예        합계
SUV
아니오  0.393939  0.132576  0.526515
예    0.367424  0.106061  0.473485
합계   0.761364  0.238636  1.000000

행으로 나눈 표 — P(사고 | SUV)
사고        아니오         예
SUV
아니오  0.748201  0.251799
예    0.776000  0.224000

열로 나눈 표 — P(SUV | 사고)
사고        아니오         예
SUV
아니오  0.517413  0.555556
예    0.482587  0.444444

P(사고 | SUV) = 28/125 = 0.224000
P(SUV)       = 125/264 = 0.473485
P(사고)       = 63/264 = 0.238636
베이즈:  P(SUV | 사고) = 0.224000 * 0.473485 / 0.238636 = 0.444444
표에서 직접: 28/63 = 0.444444

두 조건부확률의 비 = P(SUV)/P(사고) = 1.984127   (실제 비 1.984127)

import matplotlib.pyplot as plt
import pandas as pd
from statsmodels.graphics.mosaicplot import mosaic

data = {'SUV': 28*['예'] + 35*['아니오'] + 97*['예'] + 104*['아니오'],
        '사고': 28*['예'] + 35*['예']    + 97*['아니오'] + 104*['아니오']}
df = pd.DataFrame(data)

fig, ax = plt.subplots(figsize=(7, 4.5))

# 변수를 준 순서대로 나눈다.
#   첫 변수(SUV)   -> 세로로 분할. 열의 너비가 그 범주의 비율
#   두 번째(사고)  -> 각 열 안에서 가로로 분할. 조건부확률이 높이가 된다
# gap은 사각형 사이의 간격
_, rects = mosaic(df, ['SUV', '사고'], ax=ax, gap=0.02,
                  title='모자이크 그림: SUV와 사고')
plt.tight_layout()
plt.show()

# 넓이를 재려면 gap 을 0 으로 두어야 한다. gap 은 너비와 높이를 깎는다.
fig0, ax0 = plt.subplots()
_, tiles = mosaic(df, ['SUV', '사고'], ax=ax0, gap=0.0)
plt.close(fig0)

prop = pd.crosstab(df.SUV, df.사고) / len(df)
print(f"{'타일':>18}{'너비':>11}{'높이':>11}{'넓이':>11}{'상대도수':>12}")
total = 0.0
for suv in ['예', '아니오']:
    for acc in ['예', '아니오']:
        x, y, w, h = tiles[(suv, acc)]
        total += w * h
        print(f"{'SUV=' + suv + ', 사고=' + acc:>18}{w:>11.6f}{h:>11.6f}"
              f"{w * h:>11.6f}{prop.loc[suv, acc]:>12.6f}")
print(f"{'합':>18}{'':>11}{'':>11}{total:>11.6f}{prop.values.sum():>12.6f}")
print(f"넓이의 합에서 1 을 뺀 값 = {total - 1:.3e}")

# gap 을 넣으면 넓이가 깎인다
w0, h0 = tiles[('예', '예')][2], tiles[('예', '예')][3]
wg, hg = rects[('예', '예')][2], rects[('예', '예')][3]
print(f"\n(SUV=예, 사고=예) 타일")
print(f"  gap=0.00: 너비 {w0:.6f} x 높이 {h0:.6f} = {w0*h0:.6f}")
print(f"  gap=0.02: 너비 {wg:.6f} x 높이 {hg:.6f} = {wg*hg:.6f}"
      f"  ({wg*hg/(w0*h0) - 1:+.2%})")

                타일         너비         높이         넓이        상대도수
       SUV=예, 사고=예   0.473485   0.224000   0.106061    0.106061
     SUV=예, 사고=아니오   0.473485   0.776000   0.367424    0.367424
     SUV=아니오, 사고=예   0.526515   0.251799   0.132576    0.132576
   SUV=아니오, 사고=아니오   0.526515   0.748201   0.393939    0.393939
                 합                         1.000000    1.000000
넓이의 합에서 1 을 뺀 값 = 0.000e+00

(SUV=예, 사고=예) 타일
  gap=0.00: 너비 0.473485 x 높이 0.224000 = 0.106061
  gap=0.02: 너비 0.464201 x 높이 0.221053 = 0.102613  (-3.25%)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from statsmodels.graphics.mosaicplot import mosaic

rng = np.random.default_rng(7)
n = 400

# --- 독립: Y가 X와 무관하게 생성된다 ---
A = rng.choice(['예', '아니오'], n, p=[.5, .5])
B = rng.choice(['예', '아니오'], n, p=[.3, .7])
df_indep = pd.DataFrame({'X': A, 'Y': B})

# --- 종속: X가 '예'면 Y가 '예'일 확률이 크게 높아진다 ---
A2 = rng.choice(['예', '아니오'], n, p=[.5, .5])
B2 = np.where(A2 == '예',
              rng.choice(['예', '아니오'], n, p=[.7, .3]),
              rng.choice(['예', '아니오'], n, p=[.1, .9]))
df_dep = pd.DataFrame({'X': A2, 'Y': B2})

for name, d in [('독립', df_indep), ('종속', df_dep)]:
    ct = pd.crosstab(d.X, d.Y)
    print(f"=== {name} ===")
    print(ct)
    print("P(Y=예 | X=아니오) =", round(ct.loc['아니오', '예'] / ct.loc['아니오'].sum(), 3))
    print("P(Y=예 | X=예)    =", round(ct.loc['예', '예'] / ct.loc['예'].sum(), 3))

    O = ct.values.astype(float)
    chi2, p, dof, E = stats.chi2_contingency(O, correction=False)
    r = (O - E) / np.sqrt(E)                       # 피어슨 잔차
    pr, pc = O.sum(1) / O.sum(), O.sum(0) / O.sum()
    adj = r / np.sqrt(np.outer(1 - pr, 1 - pc))    # 조정잔차
    print("분할선 높이 차 =", round(abs(ct.loc['예', '예'] / ct.loc['예'].sum()
                                 - ct.loc['아니오', '예'] / ct.loc['아니오'].sum()), 3))
    print("기대도수 E =", np.round(E, 2).tolist())
    print("피어슨 잔차 =", np.round(r, 3).tolist())
    print("조정잔차    =", np.round(adj, 3).tolist())
    print(f"잔차 제곱합 {np.sum(r**2):.4f}  = chi2 {chi2:.4f}   (dof {dof}, p {p:.3g})")
    print(f"sqrt(chi2) = {np.sqrt(chi2):.4f}  = |조정잔차| {np.abs(adj).max():.4f}")
    print()

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
mosaic(df_indep, ['X', 'Y'], ax=axes[0], gap=0.02, title='독립')
mosaic(df_dep,   ['X', 'Y'], ax=axes[1], gap=0.02, title='종속')
plt.tight_layout()
plt.show()

# 앞의 SUV 자료도 같은 눈금으로 재어 본다.
O_suv = np.array([[104, 35], [97, 28]], float)        # 보기 1 의 표
chi2, p, dof, E = stats.chi2_contingency(O_suv, correction=False)
pr, pc = O_suv.sum(1) / O_suv.sum(), O_suv.sum(0) / O_suv.sum()
adj = ((O_suv - E) / np.sqrt(E)) / np.sqrt(np.outer(1 - pr, 1 - pc))
print(f"SUV 자료: chi2 {chi2:.4f}, p {p:.3f}, |조정잔차| {np.abs(adj).max():.4f}")

=== 독립 ===
Y    아니오   예
X
아니오  149  61
예    132  58
P(Y=예 | X=아니오) = 0.29
P(Y=예 | X=예)    = 0.305
분할선 높이 차 = 0.015
기대도수 E = [[147.52, 62.48], [133.48, 56.52]]
피어슨 잔차 = [[0.121, -0.187], [-0.128, 0.196]]
조정잔차    = [[0.323, -0.323], [-0.323, 0.323]]
잔차 제곱합 0.1044  = chi2 0.1044   (dof 1, p 0.747)
sqrt(chi2) = 0.3230  = |조정잔차| 0.3230

=== 종속 ===
Y    아니오    예
X
아니오  178   16
예     61  145
P(Y=예 | X=아니오) = 0.082
P(Y=예 | X=예)    = 0.704
분할선 높이 차 = 0.621
기대도수 E = [[115.92, 78.08], [123.08, 82.92]]
피어슨 잔차 = [[5.767, -7.026], [-5.596, 6.818]]
조정잔차    = [[12.666, -12.666], [-12.666, 12.666]]
잔차 제곱합 160.4208  = chi2 160.4208   (dof 1, p 9.16e-37)
sqrt(chi2) = 12.6657  = |조정잔차| 12.6657

SUV 자료: chi2 0.2799, p 0.597, |조정잔차| 0.5291

import numpy as np
from scipy import stats

obs = np.array([[40, 60], [70, 30]], float)
chi2, p, dof, exp = stats.chi2_contingency(obs, correction=False)
resid = (obs - exp) / np.sqrt(exp)

print("관측 도수\n", obs.astype(int))
print("독립일 때의 기대 도수\n", np.round(exp, 2))
print("피어슨 잔차\n", np.round(resid, 3))
print(f"\nchi2 = {chi2:.4f}   잔차 제곱합 = {np.sum(resid ** 2):.4f}   p = {p:.3e}")

관측 도수
 [[40 60]
 [70 30]]
독립일 때의 기대 도수
 [[55. 45.]
 [55. 45.]]
피어슨 잔차
 [[-2.023  2.236]
 [ 2.023 -2.236]]

chi2 = 18.1818   잔차 제곱합 = 18.1818   p = 2.008e-05

import numpy as np
import pandas as pd

tab = pd.DataFrame({"성공": [45, 180], "실패": [15, 120]},
                   index=["소규모 병원", "대규모 병원"])
n = tab.values.sum()
row_n = tab.sum(axis=1)

print(tab.to_string())
print(f"\n{'집단':>12}{'표본 수':>9}{'비중':>9}{'성공률':>9}")
for name in tab.index:
    print(f"{name:>12}{row_n[name]:>9}{row_n[name]/n:>9.3f}"
          f"{tab.loc[name, '성공']/row_n[name]:>9.3f}")

print("\n100% 누적 막대: 두 막대의 폭이 같다 → 표본 수 정보 없음")
print(f"모자이크:       폭이 {row_n['소규모 병원']/n:.3f} : {row_n['대규모 병원']/n:.3f}"
      f" 로 표본 수를 반영한다")

         성공   실패
소규모 병원   45   15
대규모 병원  180  120

          집단     표본 수       비중      성공률
      소규모 병원       60    0.167    0.750
      대규모 병원      300    0.833    0.600

100% 누적 막대: 두 막대의 폭이 같다 → 표본 수 정보 없음
모자이크:       폭이 0.167 : 0.833 로 표본 수를 반영한다

import numpy as np
from scipy import stats

# [결석 크기][치료] -> (성공, 실패)
table = np.array([[[81, 6], [234, 36]],       # 작은 결석: A, B
                  [[192, 71], [55, 25]]], float)  # 큰 결석: A, B

print(f"{'층':>10}{'치료':>6}{'성공률':>10}{'환자 수':>9}")
for i, size in enumerate(["작은 결석", "큰 결석"]):
    for j, tr in enumerate(["A", "B"]):
        s, f = table[i, j]
        print(f"{size:>10}{tr:>6}{s/(s+f):>10.4f}{int(s+f):>9}")

agg = table.sum(axis=0)
print()
for j, tr in enumerate(["A", "B"]):
    s, f = agg[j]
    print(f"{'전체':>10}{tr:>6}{s/(s+f):>10.4f}{int(s+f):>9}")

print(f"\n주변화한 2x2 표의 카이제곱 p = "
      f"{stats.chi2_contingency(agg, correction=False)[1]:.4f}")

         층    치료       성공률     환자 수
     작은 결석     A    0.9310       87
     작은 결석     B    0.8667      270
      큰 결석     A    0.7300      263
      큰 결석     B    0.6875       80

        전체     A    0.7800      350
        전체     B    0.8257      350

주변화한 2x2 표의 카이제곱 p = 0.1285

import numpy as np
import pandas as pd

# 한 범주 조합에만 자료가 거의 없는 표
obs = pd.DataFrame([[120, 90], [85, 3]],
                   index=["A", "B"], columns=["X", "Y"])
n = obs.values.sum()
print("관측도수"); print(obs.to_string())

# 모자이크 그림에서 각 칸이 차지하는 넓이 = 결합비율
print("\n칸이 차지하는 넓이 (결합비율)")
print((obs / n).round(4).to_string())

# 폭(행 비율)과 높이(행 안에서의 열 비율)
print("\n폭 = 행 비율")
print((obs.sum(axis=1) / n).round(4).to_string())
print("\n높이 = 행 안에서의 열 비율")
print(obs.div(obs.sum(axis=1), axis=0).round(4).to_string())

관측도수
     X   Y
A  120  90
B   85   3

칸이 차지하는 넓이 (결합비율)
        X       Y
A  0.4027  0.3020
B  0.2852  0.0101

폭 = 행 비율
A    0.7047
B    0.2953

높이 = 행 안에서의 열 비율
        X       Y
A  0.5714  0.4286
B  0.9659  0.0341

import numpy as np
import pandas as pd

# 만족도(순서형) × 연령대(순서형)
data = pd.DataFrame({
    "매우 불만": [30, 20, 10, 5],
    "불만":     [40, 35, 25, 15],
    "보통":     [50, 55, 60, 50],
    "만족":     [30, 45, 60, 70],
    "매우 만족": [10, 20, 35, 60],
}, index=["20대", "30대", "40대", "50대"])

print(data.to_string())
print("\n행별 조건부확률 (연령대별 만족도 분포)")
print((data.div(data.sum(axis=1), axis=0) * 100).round(1).to_string())

top2 = data[["만족", "매우 만족"]].sum(axis=1) / data.sum(axis=1)
print(f"\n'만족 이상' 비율: {dict(zip(data.index, top2.round(3)))}")

     매우 불만  불만  보통  만족  매우 만족
20대     30  40  50  30     10
30대     20  35  55  45     20
40대     10  25  60  60     35
50대      5  15  50  70     60

행별 조건부확률 (연령대별 만족도 분포)
     매우 불만    불만    보통    만족  매우 만족
20대   18.8  25.0  31.2  18.8    6.2
30대   11.4  20.0  31.4  25.7   11.4
40대    5.3  13.2  31.6  31.6   18.4
50대    2.5   7.5  25.0  35.0   30.0

'만족 이상' 비율: {'20대': 0.25, '30대': 0.371, '40대': 0.5, '50대': 0.65}

import numpy as np
import pandas as pd

rng = np.random.default_rng(6)

for r, c in [(2, 2), (3, 4), (5, 6), (8, 10)]:
    n = 500
    tab = rng.multinomial(n, np.ones(r * c) / (r * c)).reshape(r, c).astype(float)
    empty = np.sum(tab == 0)
    small = np.sum(tab < 5)
    widths = tab.sum(axis=1) / n
    print(f"{r}x{c} 표 ({r*c}칸, n={n}): 빈 칸 {empty}개, 5 미만 {small}개, "
          f"평균 칸 도수 {n/(r*c):.1f}, 최소 열 폭 {widths.min():.3f}")

2x2 표 (4칸, n=500): 빈 칸 0개, 5 미만 0개, 평균 칸 도수 125.0, 최소 열 폭 0.480
3x4 표 (12칸, n=500): 빈 칸 0개, 5 미만 0개, 평균 칸 도수 41.7, 최소 열 폭 0.288
5x6 표 (30칸, n=500): 빈 칸 0개, 5 미만 0개, 평균 칸 도수 16.7, 최소 열 폭 0.188
8x10 표 (80칸, n=500): 빈 칸 0개, 5 미만 24개, 평균 칸 도수 6.2, 최소 열 폭 0.094