평균 차이 신뢰구간 모의실험¶

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import t, norm

np.random.seed(42)          # 아래 출력과 그림을 재현하려면 고정한다

n_simulations = 100
# 표본크기도 분산도 서로 다르게 잡았다. 이런 설정에서 합동 t가 무너지고
# Welch가 버티는지 보려는 것이다.
n1, n2 = 12, 10
mu1, mu2 = 0.0, 0.5
sigma1, sigma2 = 1.0, 1.5
alpha = 0.05
method = "welch"  # 'welch' | 'pooled' | 'z_known' | 'z_plugin'

delta_true = mu1 - mu2
lowers = np.empty(n_simulations)
uppers = np.empty(n_simulations)
centers = np.empty(n_simulations)

for i in range(n_simulations):
    x = np.random.normal(mu1, sigma1, n1)
    y = np.random.normal(mu2, sigma2, n2)
    xbar, ybar = x.mean(), y.mean()
    s1, s2 = x.std(ddof=1), y.std(ddof=1)
    diff_hat = xbar - ybar
    centers[i] = diff_hat

    if method == "welch":
        se = np.sqrt(s1**2 / n1 + s2**2 / n2)
        num = (s1**2 / n1 + s2**2 / n2)**2
        den = (s1**2 / n1)**2 / (n1 - 1) + (s2**2 / n2)**2 / (n2 - 1)
        df = num / den
        crit = t.ppf(1 - alpha / 2, df=df)
    elif method == "pooled":
        df = n1 + n2 - 2
        sp2 = ((n1 - 1) * s1**2 + (n2 - 1) * s2**2) / df
        se = np.sqrt(sp2 * (1 / n1 + 1 / n2))
        crit = t.ppf(1 - alpha / 2, df=df)
    elif method == "z_known":
        se = np.sqrt(sigma1**2 / n1 + sigma2**2 / n2)
        crit = norm.ppf(1 - alpha / 2)
    else:  # z_plugin
        se = np.sqrt(s1**2 / n1 + s2**2 / n2)
        crit = norm.ppf(1 - alpha / 2)

    lowers[i] = diff_hat - crit * se
    uppers[i] = diff_hat + crit * se

covered = (lowers <= delta_true) & (delta_true <= uppers)
coverage_pct = 100.0 * covered.mean()
print(f"{method} coverage: {coverage_pct:.1f}%")

welch coverage: 97.0%

# 구간 하나를 가로선 하나로 그린다. 참값을 담은 구간은 검정, 놓친 구간은
# 빨강이다. 세로 점선이 참값이고, 빨간 선이 몇 개인지 세는 것이 곧 포함확률을
# 재는 일이다. 구간마다 길이가 다른 까닭은 표본마다 s 가 다르기 때문이다.
fig, ax = plt.subplots(figsize=(12, 12))
for i in range(n_simulations):
    color = "k" if covered[i] else "r"
    ax.plot([lowers[i], uppers[i]], [i, i], lw=2, color=color)
    ax.plot(centers[i], i, marker="o", ms=3, color=color)

ax.axvline(delta_true, linestyle="--", linewidth=1.5)
n_fail = int((~covered).sum())
ax.set_title(f"{n_simulations} {method} CIs | n1={n1}, n2={n2}, CL=95%")
ax.set_yticks([])
ax.set_xlabel("Difference of means")
plt.tight_layout()
plt.show()

from scipy.stats import t as t_dist, norm
np.random.seed(0)
n1 = n2 = 10
sigma1, sigma2 = 1.0, 3.0
mu1 = mu2 = 0.0
delta = 0.0
n_sim = 10_000
alpha = 0.05
welch_cov = pooled_cov = 0
for _ in range(n_sim):
    x = np.random.normal(mu1, sigma1, n1)
    y = np.random.normal(mu2, sigma2, n2)
    s1, s2 = x.std(ddof=1), y.std(ddof=1)
    diff = x.mean() - y.mean()
    # Welch 방법 — 두 분산이 다를 수 있다고 본다
    se_w = np.sqrt(s1**2/n1 + s2**2/n2)
    num = (s1**2/n1 + s2**2/n2)**2
    den = (s1**2/n1)**2/9 + (s2**2/n2)**2/9
    df_w = num/den
    crit_w = t_dist.ppf(0.975, df_w)
    if diff - crit_w*se_w <= delta <= diff + crit_w*se_w:
        welch_cov += 1
    # 합동 방법 — 두 분산이 같다고 본다
    sp2 = (9*s1**2 + 9*s2**2)/18
    se_p = np.sqrt(sp2*(1/10+1/10))
    crit_p = t_dist.ppf(0.975, 18)
    if diff - crit_p*se_p <= delta <= diff + crit_p*se_p:
        pooled_cov += 1
print(f"Welch: {100*welch_cov/n_sim:.1f}%")
print(f"Pooled: {100*pooled_cov/n_sim:.1f}%")

Welch: 95.0%
Pooled: 94.3%

import numpy as np
from scipy import stats

rng = np.random.default_rng(17)
M = 20_000
print(f"{'n1':>4s} {'n2':>4s} {'σ1':>4s} {'σ2':>4s} "
      f"{'합동':>8s} {'웰치':>8s}")
cases = [(10, 10, 1, 3), (10, 30, 1, 3), (30, 10, 1, 3),
         (10, 30, 3, 1), (20, 20, 1, 1), (5, 50, 1, 4)]
for n1, n2, s1, s2 in cases:
    x = rng.normal(0, s1, (M, n1))
    y = rng.normal(0, s2, (M, n2))
    v1, v2 = x.var(1, ddof=1), y.var(1, ddof=1)
    d = x.mean(1) - y.mean(1)

    sp2 = ((n1 - 1) * v1 + (n2 - 1) * v2) / (n1 + n2 - 2)
    hp = stats.t.ppf(0.975, n1 + n2 - 2) * np.sqrt(sp2 * (1 / n1 + 1 / n2))

    se = np.sqrt(v1 / n1 + v2 / n2)
    df = (v1 / n1 + v2 / n2)**2 / ((v1 / n1)**2 / (n1 - 1)
                                   + (v2 / n2)**2 / (n2 - 1))
    hw = stats.t.ppf(0.975, df) * se

    print(f"{n1:4d} {n2:4d} {s1:4d} {s2:4d} "
          f"{np.mean(np.abs(d) <= hp):8.4f} {np.mean(np.abs(d) <= hw):8.4f}")

  n1   n2   σ1   σ2       합동       웰치
  10   10    1    3   0.9417   0.9498
  10   30    1    3   0.9961   0.9528
  30   10    1    3   0.7898   0.9494
  10   30    3    1   0.7865   0.9465
  20   20    1    1   0.9528   0.9535
   5   50    1    4   1.0000   0.9496

import numpy as np
from scipy import stats

print(f"{'n':>5s} {'합동 df':>8s} {'웰치 평균 df':>12s} {'폭 비':>8s}")
rng = np.random.default_rng(1)
for n in [5, 10, 20, 50]:
    x = rng.normal(0, 1, (20_000, n))
    y = rng.normal(0, 1, (20_000, n))
    v1, v2 = x.var(1, ddof=1), y.var(1, ddof=1)
    df = (v1 / n + v2 / n)**2 / ((v1 / n)**2 / (n - 1) + (v2 / n)**2 / (n - 1))
    tw = stats.t.ppf(0.975, df).mean()
    tp = stats.t.ppf(0.975, 2 * n - 2)
    print(f"{n:5d} {2 * n - 2:8d} {df.mean():12.1f} {tw / tp:8.4f}")

    n    합동 df     웰치 평균 df      폭 비
    5        8          6.8   1.0373
   10       18         16.5   1.0074
   20       38         36.3   1.0017
   50       98         96.1   1.0003

import numpy as np
from scipy import stats

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

def welch_cover(x, y, n1, n2, true=0.0):
    v1, v2 = x.var(1, ddof=1), y.var(1, ddof=1)
    se = np.sqrt(v1 / n1 + v2 / n2)
    df = (v1 / n1 + v2 / n2)**2 / ((v1 / n1)**2 / (n1 - 1)
                                   + (v2 / n2)**2 / (n2 - 1))
    h = stats.t.ppf(0.975, df) * se
    d = x.mean(1) - y.mean(1)
    below = np.mean(d + h < true)
    above = np.mean(d - h > true)
    return 1 - below - above, below, above

cases = [
    ("정규 / 정규", lambda m: rng.normal(0, 1, m), lambda m: rng.normal(0, 1, m), 15, 15),
    ("지수 / 지수", lambda m: rng.exponential(1, m) - 1,
     lambda m: rng.exponential(1, m) - 1, 15, 15),
    ("지수 / 정규", lambda m: rng.exponential(1, m) - 1,
     lambda m: rng.normal(0, 1, m), 15, 15),
    ("지수 / 지수 (n=50)", lambda m: rng.exponential(1, m) - 1,
     lambda m: rng.exponential(1, m) - 1, 50, 50),
    ("t(3) / t(3)", lambda m: rng.standard_t(3, m),
     lambda m: rng.standard_t(3, m), 20, 20),
]
for name, g1, g2, n1, n2 in cases:
    c, lo, hi = welch_cover(g1((M, n1)), g2((M, n2)), n1, n2)
    print(f"{name:20s} 포함 {c:.4f}   아래 누락 {lo:.4f}   위 누락 {hi:.4f}")

정규 / 정규              포함 0.9513   아래 누락 0.0239   위 누락 0.0248
지수 / 지수              포함 0.9567   아래 누락 0.0209   위 누락 0.0224
지수 / 정규              포함 0.9456   아래 누락 0.0406   위 누락 0.0138
지수 / 지수 (n=50)       포함 0.9526   아래 누락 0.0242   위 누락 0.0233
t(3) / t(3)          포함 0.9573   아래 누락 0.0222   위 누락 0.0205

import numpy as np
from scipy import stats

rng = np.random.default_rng(88)
x = rng.lognormal(1.0, 0.5, 40)          # 집단 1
y = rng.lognormal(0.7, 0.5, 45)          # 집단 2
z = stats.norm.ppf(0.975)

# 1) 평균 차이 (웰치)
n1, n2 = len(x), len(y)
v1, v2 = x.var(ddof=1), y.var(ddof=1)
se = np.sqrt(v1 / n1 + v2 / n2)
df = (v1 / n1 + v2 / n2)**2 / ((v1 / n1)**2 / (n1 - 1) + (v2 / n2)**2 / (n2 - 1))
h = stats.t.ppf(0.975, df) * se
print(f"평균 차이   {x.mean() - y.mean():7.4f}  "
      f"({x.mean() - y.mean() - h:7.4f}, {x.mean() - y.mean() + h:7.4f})")

# 2) 로그 척도 차이 → 기하평균의 비
lx, ly = np.log(x), np.log(y)
sl = np.sqrt(lx.var(ddof=1) / n1 + ly.var(ddof=1) / n2)
dfl = (lx.var(ddof=1) / n1 + ly.var(ddof=1) / n2)**2 / (
    (lx.var(ddof=1) / n1)**2 / (n1 - 1) + (ly.var(ddof=1) / n2)**2 / (n2 - 1))
hl = stats.t.ppf(0.975, dfl) * sl
dl = lx.mean() - ly.mean()
print(f"기하평균 비 {np.exp(dl):7.4f}  "
      f"({np.exp(dl - hl):7.4f}, {np.exp(dl + hl):7.4f})")

# 3) 코헨의 d (합동 표준편차 기준) + 부트스트랩 구간
sp = np.sqrt(((n1 - 1) * v1 + (n2 - 1) * v2) / (n1 + n2 - 2))
dcoh = (x.mean() - y.mean()) / sp
B = 20_000
bx = x[rng.integers(0, n1, (B, n1))]
by = y[rng.integers(0, n2, (B, n2))]
spb = np.sqrt(((n1 - 1) * bx.var(1, ddof=1) + (n2 - 1) * by.var(1, ddof=1))
              / (n1 + n2 - 2))
db = (bx.mean(1) - by.mean(1)) / spb
print(f"코헨의 d    {dcoh:7.4f}  "
      f"({np.percentile(db, 2.5):7.4f}, {np.percentile(db, 97.5):7.4f})")

# 4) 중앙값 차이 (부트스트랩)
mb = np.median(bx, 1) - np.median(by, 1)
print(f"중앙값 차이 {np.median(x) - np.median(y):7.4f}  "
      f"({np.percentile(mb, 2.5):7.4f}, {np.percentile(mb, 97.5):7.4f})")

평균 차이    0.8968  ( 0.3329,  1.4608)
기하평균 비  1.3590  ( 1.1250,  1.6418)
코헨의 d     0.7099  ( 0.3063,  1.1280)
중앙값 차이  0.6539  (-0.0641,  1.6675)