집단 비교 보기¶

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# 그림에 한글을 쓰므로 한글 글꼴을 지정한다. 맥이면 'Apple SD Gothic Neo',
# 윈도우면 'Malgun Gothic', 리눅스면 'NanumGothic' 정도가 무난하다.
# 글꼴을 바꾸면 마이너스 기호가 깨지므로 unicode_minus 도 함께 꺼 준다.
plt.rcParams['font.family'] = 'Apple SD Gothic Neo'
plt.rcParams['axes.unicode_minus'] = False

np.random.seed(42)
airlines = ['American', 'Delta', 'Southwest', 'United']
n_obs = 100

# 감마분포는 0 이상이고 오른쪽으로 치우쳐 있어 지연율 모형에 적합하다.
# (shape, scale)로 모양과 척도를 따로 정한다. 평균은 shape x scale이다.
params = {
    'American':  (2,   3),    # 평균 6.0 — 지연이 잦은 편
    'Delta':     (1.5, 2.5),  # 평균 3.75 — 중간
    'Southwest': (1.2, 2),    # 평균 2.4 — 지연이 적음
    'United':    (1.8, 3.2),  # 평균 5.76 — American과 비슷하나 퍼짐이 크다
}

data_list = []
for airline in airlines:
    shape, scale = params[airline]
    delays = np.random.gamma(shape=shape, scale=scale, size=n_obs)
    for delay in delays:
        data_list.append({'airline': airline, 'pct_carrier_delay': delay})

airline_stats = pd.DataFrame(data_list)

# 그림을 그리기 전에 숫자로 확인한다
print(airline_stats.groupby('airline')['pct_carrier_delay']
      .agg(['count', 'mean', 'median', 'std', 'max']).round(2))

           count  mean  median   std    max
airline
American     100  5.76    5.09  3.68  21.89
Delta        100  3.77    3.08  2.82  14.69
Southwest    100  2.31    1.67  2.08   9.07
United       100  6.18    5.09  4.92  24.20

fig, ax = plt.subplots(figsize=(8, 5))

# by=     이 열의 값으로 집단을 나눈다
# column= 분포를 볼 열
airline_stats.boxplot(by='airline', column='pct_carrier_delay', ax=ax)

ax.set_xlabel('항공사')
ax.set_ylabel('지연 항공편의 일간 비율 (%)')
ax.set_title('항공사별 지연 비교: 상자그림')
plt.suptitle('')      # pandas가 자동으로 붙이는 제목을 지운다
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(8, 5))

# inner='quartile' 로 바이올린 안에 사분위수 선을 그려 넣으면
# 상자그림의 정보를 잃지 않으면서 밀도까지 함께 볼 수 있다.
sns.violinplot(data=airline_stats, x='airline', y='pct_carrier_delay',
               ax=ax, inner='quartile', color='lightblue')

ax.set_xlabel('항공사')
ax.set_ylabel('지연 항공편의 일간 비율 (%)')
ax.set_title('항공사별 지연 비교: 바이올린 그림')
plt.tight_layout()
plt.show()

from scipy import stats

samp = {a: airline_stats.loc[airline_stats.airline == a,
                             'pct_carrier_delay'].to_numpy()
        for a in airlines}

print(f"{'항공사':<11}{'참 평균':>9}{'참 sd':>9}{'참 중앙값':>11}"
      f"{'표본 평균':>11}{'표본 sd':>10}{'표본 중앙값':>12}")
for a in airlines:
    sh, sc = params[a]
    x = samp[a]
    print(f"{a:<11}{sh*sc:>9.4f}{np.sqrt(sh)*sc:>9.4f}"
          f"{stats.gamma.ppf(0.5, sh, scale=sc):>11.4f}"
          f"{x.mean():>11.4f}{x.std(ddof=1):>10.4f}{np.median(x):>12.4f}")

sd_a, sd_u = np.sqrt(2)*3, np.sqrt(1.8)*3.2
print(f"\n참 sd:  American {sd_a:.4f},  United {sd_u:.4f}  ->  비 {sd_u/sd_a:.4f}")
obs = samp['United'].std(ddof=1)/samp['American'].std(ddof=1)
print(f"표본 sd 비 = {obs:.4f}")

rng = np.random.default_rng(1)
ratios = np.array([rng.gamma(1.8, 3.2, 100).std(ddof=1)
                   / rng.gamma(2, 3, 100).std(ddof=1) for _ in range(20000)])
print(f"\n같은 모형에서 n=100 으로 20000 번 다시 뽑았을 때 표본 sd 비의 분포")
print(f"  중앙값 {np.median(ratios):.4f},  95% 구간 "
      f"[{np.quantile(ratios,0.025):.4f}, {np.quantile(ratios,0.975):.4f}]")
print(f"  관측된 {obs:.4f} 는 이 분포의 {100*(ratios<obs).mean():.1f} 분위수")

fig, (a1, a2) = plt.subplots(1, 2, figsize=(12, 3.6))
grid = np.linspace(0, 25, 500)
cols = {"American": "#1565C0", "United": "#E65100"}
for a in ["American", "United"]:
    sh, sc = params[a]
    a1.plot(grid, stats.gamma.pdf(grid, sh, scale=sc), color=cols[a], lw=2,
            label=f"{a} 참 분포 (sd = {np.sqrt(sh)*sc:.3f})")
    a1.hist(samp[a], bins=25, density=True, color=cols[a], alpha=0.25)
a1.set_xlabel("지연 비율 (%)")
a1.set_ylabel("밀도")
a1.set_title("참 분포는 퍼짐이 거의 같다", fontsize=11)
a1.legend(fontsize=9)

a2.hist(ratios, bins=60, color="#DCEBFB", edgecolor="#90A4AE", lw=0.3)
a2.axvline(sd_u/sd_a, color="#1565C0", lw=2, label=f"참 비 {sd_u/sd_a:.4f}")
a2.axvline(obs, color="#D32F2F", lw=2, label=f"이 표본의 비 {obs:.4f}")
for q in [0.025, 0.975]:
    a2.axvline(np.quantile(ratios, q), color="#37474F", ls=":", lw=1.2)
a2.set_xlabel("표본 sd 비 (United / American)")
a2.set_ylabel("빈도")
a2.set_title("n = 100 이면 비가 0.74 ~ 1.37 사이를 오간다", fontsize=11)
a2.legend(fontsize=9)
for ax in (a1, a2):
    ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
plt.show()

항공사             참 평균     참 sd      참 중앙값      표본 평균     표본 sd      표본 중앙값
American      6.0000   4.2426     5.0350     5.7575    3.6792      5.0924
Delta         3.7500   3.0619     2.9575     3.7739    2.8208      3.0805
Southwest     2.4000   2.1909     1.7759     2.3074    2.0814      1.6691
United        5.7600   4.2933     4.7355     6.1765    4.9154      5.0893

참 sd:  American 4.2426,  United 4.2933  ->  비 1.0119
표본 sd 비 = 1.3360

같은 모형에서 n=100 으로 20000 번 다시 뽑았을 때 표본 sd 비의 분포
  중앙값 1.0118,  95% 구간 [0.7402, 1.3745]
  관측된 1.3360 는 이 분포의 96.2 분위수

np.random.seed(123)
zip_codes = [98188, 98105, 98108, 98126]
n_homes = 150

data_list = []
for zip_code in zip_codes:
    # 두 우편번호는 30만 달러대, 나머지 둘은 45만 달러대로 설정한다
    base_price = 300_000 if zip_code in [98105, 98108] else 450_000
    prices = np.random.normal(base_price, 100_000, n_homes)
    # 음수나 비현실적인 값이 나오지 않도록 아래위를 잘라 낸다.
    # 정규분포는 양쪽으로 무한히 뻗으므로 가격 자료에는 이 처리가 필요하다.
    prices = np.clip(prices, 50_000, 2_000_000)
    for price in prices:
        data_list.append({'ZipCode': str(zip_code), 'TaxAssessedValue': price})

housing = pd.DataFrame(data_list)

# 우편번호별 요약. 금액만 천 달러 단위로 바꾸고 count는 개수 그대로 둔다.
# 표 전체를 1000으로 나누면 count(150)까지 0.15가 되어 버리므로 열을 나누어 처리한다.
summary = (housing.groupby('ZipCode')['TaxAssessedValue']
           .agg(['count', 'mean', 'median', 'std']))
summary[['mean', 'median', 'std']] = (summary[['mean', 'median', 'std']] / 1000).round(1)
print(summary)

fig, ax = plt.subplots(figsize=(8, 5))
housing.boxplot(by='ZipCode', column='TaxAssessedValue', ax=ax)
ax.set_xlabel('우편번호')
ax.set_ylabel('과세평가액 (달러)')
ax.set_title('동네별 주택 가치')
plt.suptitle('')
plt.tight_layout()
plt.show()

         count   mean  median    std
ZipCode
98105      150  288.9   290.9   93.9
98108      150  290.8   287.2   93.2
98126      150  459.3   464.0  101.9
98188      150  455.8   454.0  109.4

from scipy import stats

z = stats.norm.ppf(0.975)
sd, n = 100_000.0, 150
SE = sd / np.sqrt(n)
print(f"z_0.975 = {z:.6f},  sd = {sd:,.0f},  n = {n}  ->  SE = {SE:,.1f}")
print(f"구간이 겹치지 않을 조건  |d| > z(SE+SE) = {2*z:.4f} SE = {2*z*SE:,.0f}")
print(f"차의 검정이 기각할 조건  |d| > z sqrt(SE^2+SE^2) = {z*np.sqrt(2):.4f} SE "
      f"= {z*np.sqrt(2)*SE:,.0f}")
print(f"두 문턱의 비 = sqrt(2) = {2*z/(z*np.sqrt(2)):.6f}")
print(f"'겹치면 차이 없음' 규칙의 실효 유의수준 = "
      f"{2*(1 - stats.norm.cdf(2*z/np.sqrt(2))):.6f}")

rng = np.random.default_rng(7)
R = 20000
print(f"\n참 차이를 바꾸어 가며 {R} 번 모의")
print(f"{'참 차이':>9}{'/SE':>7}{'구간 겹침':>11}{'t 기각':>9}{'겹치는데 기각':>14}")
for delta in [0, 10_000, 23_000, 25_000, 32_000, 45_000]:
    A = rng.normal(300_000, sd, (R, n))
    B = rng.normal(300_000 + delta, sd, (R, n))
    ma, mb = A.mean(1), B.mean(1)
    sa, sb = A.std(1, ddof=1), B.std(1, ddof=1)
    ha, hb = z*sa/np.sqrt(n), z*sb/np.sqrt(n)
    overlap = np.abs(ma - mb) < ha + hb
    t = (mb - ma)/np.sqrt((sa**2 + sb**2)/n)
    rej = np.abs(t) > stats.t.ppf(0.975, 2*n - 2)
    print(f"{delta:>9,}{delta/SE:>7.2f}{overlap.mean():>11.4f}{rej.mean():>9.4f}"
          f"{(overlap & rej).mean():>14.4f}")

# 한 번의 보기
rng2 = np.random.default_rng(6)
a = rng2.normal(300_000, sd, n)
b = rng2.normal(325_000, sd, n)
ha = z*a.std(ddof=1)/np.sqrt(n)
hb = z*b.std(ddof=1)/np.sqrt(n)
dm = b.mean() - a.mean()
se_d = np.sqrt(a.var(ddof=1)/n + b.var(ddof=1)/n)
tt = stats.ttest_ind(a, b, equal_var=False)
print(f"\n한 표본의 예 (참 차이 25,000)")
print(f"  A 평균 {a.mean():,.0f}  95% CI [{a.mean()-ha:,.0f}, {a.mean()+ha:,.0f}]")
print(f"  B 평균 {b.mean():,.0f}  95% CI [{b.mean()-hb:,.0f}, {b.mean()+hb:,.0f}]")
print(f"  두 구간이 겹치는가: {abs(dm) < ha + hb}")
print(f"  차이 {dm:,.0f},  차의 95% CI "
      f"[{dm - z*se_d:,.0f}, {dm + z*se_d:,.0f}],  p = {tt.pvalue:.4f}")

fig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 3.8),
                             gridspec_kw={"width_ratios": [1.3, 1]})
a1.boxplot([a/1000, b/1000], widths=0.5)
a1.set_xticklabels(["A 동네", "B 동네"])
a1.set_ylabel("과세평가액 (천 달러)")
a1.set_title("상자가 거의 완전히 겹친다", fontsize=11)

a2.errorbar([0, 1], [a.mean()/1000, b.mean()/1000],
            yerr=[ha/1000, hb/1000], fmt="o", color="#1565C0", capsize=6, ms=7)
a2.axhline(a.mean()/1000, color="#90A4AE", ls=":", lw=1)
a2.axhline(b.mean()/1000, color="#90A4AE", ls=":", lw=1)
a2.set_xlim(-0.6, 1.9)
a2.set_xticks([0, 1])
a2.set_xticklabels(["A 평균", "B 평균"])
a2.axhspan((b.mean()-hb)/1000, (a.mean()+ha)/1000, color="#FFE0B2", zorder=0)
a2.annotate(f"두 구간이 겹치는 띠\n({(a.mean()+ha-b.mean()+hb)/1000:.1f} 천 달러)",
            (1.15, (a.mean()+ha+b.mean()-hb)/2000), fontsize=9, color="#E65100",
            va="center")
a2.set_title(f"그래도 $p$ = {tt.pvalue:.4f}", fontsize=11)
for ax in (a1, a2):
    ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
plt.show()

z_0.975 = 1.959964,  sd = 100,000,  n = 150  ->  SE = 8,165.0
구간이 겹치지 않을 조건  |d| > z(SE+SE) = 3.9199 SE = 32,006
차의 검정이 기각할 조건  |d| > z sqrt(SE^2+SE^2) = 2.7718 SE = 22,632
두 문턱의 비 = sqrt(2) = 1.414214
'겹치면 차이 없음' 규칙의 실효 유의수준 = 0.005575

참 차이를 바꾸어 가며 20000 번 모의
     참 차이    /SE      구간 겹침     t 기각       겹치는데 기각
        0   0.00     0.9931   0.0527        0.0458
   10,000   1.22     0.9714   0.1364        0.1077
   23,000   2.82     0.7802   0.5060        0.2862
   25,000   3.06     0.7351   0.5726        0.3078
   32,000   3.92     0.4951   0.7867        0.2819
   45,000   5.51     0.1281   0.9747        0.1028

한 표본의 예 (참 차이 25,000)
  A 평균 300,444  95% CI [284,446, 316,442]
  B 평균 330,285  95% CI [314,190, 346,381]
  두 구간이 겹치는가: True
  차이 29,841,  차의 95% CI [7,147, 52,535],  p = 0.0104

np.random.seed(456)
grades = ['A', 'B', 'C', 'D', 'E', 'F', 'G']
n_per_grade = 100

data_list = []
for grade in grades:
    grade_idx = ord(grade) - ord('A')          # A=0, B=1, ... G=6
    # 등급이 낮아질수록 평균 소득은 내려가고(-8천/등급)
    # 퍼짐은 커진다(+5천/등급). 두 변화를 동시에 준 것이 이 보기의 요점이다.
    base_income = 80_000 - grade_idx * 8_000
    income_std  = 15_000 + grade_idx * 5_000
    incomes = np.random.normal(base_income, income_std, n_per_grade)
    incomes = np.clip(incomes, 10_000, 200_000)
    for income in incomes:
        data_list.append({'grade': grade, 'income': income})

loans = pd.DataFrame(data_list)

# 등급별 요약. 여기서는 count를 뽑지 않으므로 표 전체를 나눠도 안전하다.
print((loans.groupby('grade')['income']
       .agg(['mean', 'median', 'std']) / 1000).round(1))

fig, ax = plt.subplots(figsize=(10, 5))
sns.violinplot(data=loans, x='grade', y='income', ax=ax, color='lightgreen')
ax.set_xlabel('대출 등급 (A=최상, G=최하)')
ax.set_ylabel('연 소득 (달러)')
ax.set_title('신용등급별 소득 분포')
plt.tight_layout()
plt.show()

       mean  median   std
grade
A      82.2    81.2  13.9
B      69.8    69.5  21.9
C      67.5    68.8  24.7
D      52.5    52.1  29.2
E      53.5    49.1  32.6
F      43.1    38.2  32.4
G      36.2    23.9  32.2

from scipy.integrate import quad

mA, sA = 80_000., 15_000.
mB, sB = 72_000., 20_000.
delta = mA - mB
print(f"A등급 N({mA:,.0f}, {sA:,.0f}^2),  B등급 N({mB:,.0f}, {sB:,.0f}^2)")
print(f"참 차이 delta = {delta:,.0f}")

ovl = quad(lambda x: min(stats.norm.pdf(x, mA, sA), stats.norm.pdf(x, mB, sB)),
           0, 200_000, limit=200)[0]
print(f"두 분포의 겹침 계수 OVL = int min(f_A, f_B) dx = {ovl:.4f}")
print(f"  -> 두 분포가 {100*ovl:.1f}% 겹친다. 아무리 n 을 키워도 이 값은 그대로다.")
print(f"P(B 쪽에서 뽑은 사람이 A 쪽 사람보다 소득이 높을 확률) = "
      f"{stats.norm.cdf(-delta/np.sqrt(sA**2+sB**2)):.4f}")

print(f"\n{'n(집단당)':>10}{'차의 SE':>11}{'비중심 모수':>13}{'검정력':>9}{'겹침':>9}")
for n in [10, 25, 50, 100, 200, 400]:
    se = np.sqrt(sA**2/n + sB**2/n)
    ncp = delta/se
    dfw = (sA**2/n + sB**2/n)**2 / ((sA**2/n)**2/(n-1) + (sB**2/n)**2/(n-1))
    crit = stats.t.ppf(0.975, dfw)
    power = stats.nct.sf(crit, dfw, ncp) + stats.nct.cdf(-crit, dfw, ncp)
    print(f"{n:>10}{se:>11,.0f}{ncp:>13.4f}{power:>9.4f}{ovl:>9.4f}")

rng = np.random.default_rng(11)
print(f"\n모의로 확인 (20000 회)")
for n in [25, 100, 400]:
    R = 20000
    A = rng.normal(mA, sA, (R, n))
    B = rng.normal(mB, sB, (R, n))
    t, p = stats.ttest_ind(A, B, axis=1, equal_var=False)
    print(f"  n={n:>4}: 기각 비율 {np.mean(p < 0.05):.4f}")

fig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 3.6))
g = np.linspace(10_000, 150_000, 600)
a1.plot(g/1000, stats.norm.pdf(g, mA, sA), color="#1565C0", lw=2, label="A등급")
a1.plot(g/1000, stats.norm.pdf(g, mB, sB), color="#E65100", lw=2, label="B등급")
a1.fill_between(g/1000, np.minimum(stats.norm.pdf(g, mA, sA),
                                   stats.norm.pdf(g, mB, sB)),
                color="#CFD8DC", label=f"겹침 {100*ovl:.1f}%")
a1.set_xlabel("연 소득 (천 달러)")
a1.set_ylabel("밀도")
a1.set_title("두 분포는 70% 넘게 겹친다", fontsize=11)
a1.legend(fontsize=9)

ns = np.array([10, 25, 50, 100, 200, 400])
pw = []
for n in ns:
    se = np.sqrt(sA**2/n + sB**2/n)
    ncp = delta/se
    dfw = (sA**2/n + sB**2/n)**2/((sA**2/n)**2/(n-1) + (sB**2/n)**2/(n-1))
    crit = stats.t.ppf(0.975, dfw)
    pw.append(stats.nct.sf(crit, dfw, ncp) + stats.nct.cdf(-crit, dfw, ncp))
a2.plot(ns, pw, "o-", color="#1565C0", lw=2, label="t 검정의 검정력")
a2.axhline(ovl, color="#E65100", ls="--", lw=2, label=f"겹침 {ovl:.3f} (n 과 무관)")
a2.axhline(0.8, color="#90A4AE", ls=":", lw=1.2)
a2.set_xlabel("집단당 표본 크기 $n$")
a2.set_ylabel("확률")
a2.set_ylim(0, 1.05)
a2.set_title("겹침은 그대로인데 검정력만 오른다", fontsize=11)
a2.legend(fontsize=9)
for ax in (a1, a2):
    ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
plt.show()

A등급 N(80,000, 15,000^2),  B등급 N(72,000, 20,000^2)
참 차이 delta = 8,000
두 분포의 겹침 계수 OVL = int min(f_A, f_B) dx = 0.7848
  -> 두 분포가 78.5% 겹친다. 아무리 n 을 키워도 이 값은 그대로다.
P(B 쪽에서 뽑은 사람이 A 쪽 사람보다 소득이 높을 확률) = 0.3745

    n(집단당)      차의 SE       비중심 모수      검정력       겹침
        10      7,906       1.0119   0.1591   0.7848
        25      5,000       1.6000   0.3468   0.7848
        50      3,536       2.2627   0.6098   0.7848
       100      2,500       3.2000   0.8894   0.7848
       200      1,768       4.5255   0.9947   0.7848
       400      1,250       6.4000   1.0000   0.7848

모의로 확인 (20000 회)
  n=  25: 기각 비율 0.3475
  n= 100: 기각 비율 0.8856
  n= 400: 기각 비율 1.0000

import numpy as np

rng = np.random.default_rng(0)
B, n = 20_000, 30

print("참 평균이 모두 0 으로 같은 K개 집단에서, 관측 평균 1위 집단의 값")
print(f"{'K':>5}{'1위의 평균':>12}{'표준오차 단위':>15}")
for K in (3, 5, 10, 20, 50):
    means = rng.normal(0, 1, (B, K, n)).mean(axis=2)
    best = means.max(axis=1)
    print(f"{K:>5}{best.mean():>+12.4f}{best.mean() * np.sqrt(n):>+15.2f}")

참 평균이 모두 0 으로 같은 K개 집단에서, 관측 평균 1위 집단의 값
    K      1위의 평균        표준오차 단위
    3     +0.1556          +0.85
    5     +0.2122          +1.16
   10     +0.2801          +1.53
   20     +0.3413          +1.87
   50     +0.4114          +2.25

import numpy as np

rng = np.random.default_rng(0)
K, n = 12, 25
data = rng.normal(0, 1, (K, n))        # 12개 집단, 참 평균은 모두 0
order = np.argsort(-data.mean(1))

print("참 평균이 모두 같은 12개 집단을 관측 평균으로 정렬하면")
print(f"  {np.round(data.mean(1)[order], 3)}")
gap = data.mean(1)[order][0] - data.mean(1)[order][-1]
print(f"\n  1위와 꼴찌의 차이 {gap:.4f}")
print(f"  차이의 표준오차 {np.sqrt(2) / np.sqrt(n):.4f}")
print(f"  → 표준오차 {gap / (np.sqrt(2) / np.sqrt(n)):.2f} 배")

참 평균이 모두 같은 12개 집단을 관측 평균으로 정렬하면
  [ 0.331  0.087  0.071 -0.015 -0.018 -0.021 -0.042 -0.047 -0.073 -0.166
 -0.214 -0.324]

  1위와 꼴찌의 차이 0.6553
  차이의 표준오차 0.2828
  → 표준오차 2.32 배