대응 평균 신뢰구간의 포함확률 모의실험¶

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

n_simulations = 100
n = 12
mu_x, mu_y = 0.5, 0.0
sigma_x, sigma_y = 1.0, 1.2
rho = 0.6
alpha = 0.05
method = "t"  # 't' | 'z_known' | 'z_plugin'

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

delta_true = mu_x - mu_y
var_d_true = sigma_x**2 + sigma_y**2 - 2 * rho * sigma_x * sigma_y
sigma_d_true = np.sqrt(var_d_true)

# 상관이 있는 짝 (X, Y)를 만들어야 하므로 공분산행렬을 세우고
# Cholesky 분해 L을 쓴다. 독립인 표준정규 z에 L을 곱하면
# 공분산이 Sigma인 자료가 된다. Cov(Lz) = L L^T = Sigma 이기 때문이다.
cov = rho * sigma_x * sigma_y
Sigma = np.array([[sigma_x**2, cov], [cov, sigma_y**2]])
L = np.linalg.cholesky(Sigma)

df = n - 1
t_star = t.ppf(1 - alpha / 2, df=df)
z_star = norm.ppf(1 - alpha / 2)

lowers = np.empty(n_simulations)
uppers = np.empty(n_simulations)
centers = np.empty(n_simulations)

for i in range(n_simulations):
    z_vals = rng.standard_normal(size=(2, n))
    xy = (L @ z_vals).T
    x = xy[:, 0] + mu_x
    y = xy[:, 1] + mu_y
    d = x - y
    dbar = d.mean()
    s_d = d.std(ddof=1)

    if method == "t":
        se, crit = s_d / np.sqrt(n), t_star
    elif method == "z_known":
        se, crit = sigma_d_true / np.sqrt(n), z_star
    else:
        se, crit = s_d / np.sqrt(n), z_star

    lowers[i] = dbar - crit * se
    uppers[i] = dbar + crit * se
    centers[i] = dbar

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

# (2) rho 를 바꾸면 폭만 달라지고 포함률은 그대로다.
print(f"\n{'rho':>6}{'sigma_D':>10}{'기대폭':>10}{'모의 평균폭':>12}{'포함':>6}")
for r in (-0.3, 0.0, 0.6, 0.9):
    sd_r = np.sqrt(sigma_x**2 + sigma_y**2 - 2 * r * sigma_x * sigma_y)
    Lr = np.linalg.cholesky(np.array([[sigma_x**2, r * sigma_x * sigma_y],
                                      [r * sigma_x * sigma_y, sigma_y**2]]))
    rng_r = np.random.default_rng(42)
    hit, wsum = 0, 0.0
    for _ in range(n_simulations):
        d_r = (Lr @ rng_r.standard_normal((2, n))).T
        d_r = (d_r[:, 0] + mu_x) - (d_r[:, 1] + mu_y)
        half = t_star * d_r.std(ddof=1) / np.sqrt(n)
        hit += abs(d_r.mean() - delta_true) <= half
        wsum += 2 * half
    # 기대 폭 = 2 t* c4 sigma_D / sqrt(n),  c4 = E[S]/sigma
    c4 = np.exp(gammaln(n / 2) - gammaln((n - 1) / 2)) * np.sqrt(2 / (n - 1))
    exp_w = 2 * t_star * c4 * sd_r / np.sqrt(n)
    print(f"{r:>6.1f}{sd_r:>10.5f}{exp_w:>10.5f}{wsum / n_simulations:>12.5f}{hit:>6d}")

Paired t coverage: 95.0%

   rho   sigma_D       기대폭      모의 평균폭    포함
  -0.3   1.77764   2.20822     2.13020    95
   0.0   1.56205   1.94041     1.88188    94
   0.6   1.00000   1.24222     1.23565    95
   0.9   0.52915   0.65732     0.67048    98

# 구간 하나를 가로선 하나로 그린다. 참값을 담은 구간은 검정, 놓친 구간은
# 빨강이다. 세로 점선이 참값이고, 빨간 선이 몇 개인지 세는 것이 곧 포함확률을
# 재는 일이다. 구간마다 길이가 다른 까닭은 표본마다 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, color="r")
n_fail = int((~covered).sum())
ax.set_title(f"{n_simulations} Paired {method} CIs | n={n}, rho={rho}, CL=95%")
ax.set_yticks([])
ax.set_xlabel("Mean difference")
plt.tight_layout()
plt.show()

import numpy as np
from scipy import stats

rng = np.random.default_rng(77)
n, M = 20, 20_000
print(f"{'ρ':>6s} {'대응':>8s} {'독립(오분석)':>13s} "
      f"{'대응 폭':>9s} {'독립 폭':>9s}")
for rho in [-0.5, -0.2, 0.0, 0.3, 0.6, 0.9]:
    L = np.linalg.cholesky(np.array([[1, rho], [rho, 1]]))
    z = rng.standard_normal((M, n, 2)) @ L.T
    x, y = 10 + z[:, :, 0], 12 + z[:, :, 1]     # 참 차이 -2
    d = x - y

    h = stats.t.ppf(0.975, n - 1) * d.std(1, ddof=1) / np.sqrt(n)
    cov_p = np.mean((d.mean(1) - h <= -2) & (-2 <= d.mean(1) + h))

    sp = np.sqrt((x.var(1, ddof=1) + y.var(1, ddof=1)) / 2)
    hi = stats.t.ppf(0.975, 2 * n - 2) * sp * np.sqrt(2 / n)
    dd = x.mean(1) - y.mean(1)
    cov_i = np.mean((dd - hi <= -2) & (-2 <= dd + hi))

    print(f"{rho:6.1f} {cov_p:8.4f} {cov_i:13.4f} "
          f"{2 * h.mean():9.4f} {2 * hi.mean():9.4f}")

     ρ       대응       독립(오분석)      대응 폭      독립 폭
  -0.5   0.9500        0.8906    1.6001    1.2694
  -0.2   0.9485        0.9258    1.4353    1.2739
   0.0   0.9484        0.9493    1.3068    1.2708
   0.3   0.9495        0.9788    1.0913    1.2714
   0.6   0.9533        0.9971    0.8266    1.2673
   0.9   0.9506        1.0000    0.4132    1.2648

import numpy as np
from scipy import stats

def hl_ci(d, alpha=0.05):
    n = len(d)
    w = np.add.outer(d, d)[np.triu_indices(n)] / 2     # 왈시 평균
    w.sort()
    mu = n * (n + 1) / 4
    sd = np.sqrt(n * (n + 1) * (2 * n + 1) / 24)
    k = max(int(np.floor(mu - stats.norm.ppf(1 - alpha / 2) * sd)), 0)
    return w[k], w[len(w) - 1 - k]

# 지수분포(1)의 세 중심
pm = stats.gamma.ppf(0.5, 2) / 2                       # 유사중앙값
print(f"Exp(1)  평균 1.0000  중앙값 {np.log(2):.4f}  유사중앙값 {pm:.4f}")

rng = np.random.default_rng(5)
n, M = 15, 4_000
ht = hm = hp = 0
for _ in range(M):
    d = rng.exponential(1.0, n)
    h = stats.t.ppf(0.975, n - 1) * d.std(ddof=1) / np.sqrt(n)
    ht += (d.mean() - h <= 1.0 <= d.mean() + h)
    lo, hi = hl_ci(d)
    hm += (lo <= 1.0 <= hi)
    hp += (lo <= pm <= hi)
print(f"t 구간이 평균을 담을 확률        {ht / M:.4f}")
print(f"HL 구간이 평균을 담을 확률       {hm / M:.4f}")
print(f"HL 구간이 유사중앙값을 담을 확률  {hp / M:.4f}")

Exp(1)  평균 1.0000  중앙값 0.6931  유사중앙값 0.8392
t 구간이 평균을 담을 확률        0.9120
HL 구간이 평균을 담을 확률       0.8852
HL 구간이 유사중앙값을 담을 확률  0.9430

import numpy as np
from scipy import stats

rng = np.random.default_rng(99)
d = np.round(rng.normal(3.0, 1.5, 20), 2)
d2 = d.copy()
d2[0] = 25.0                                   # 기록 오류 한 건

def t_ci(v):
    n = len(v)
    h = stats.t.ppf(0.975, n - 1) * v.std(ddof=1) / np.sqrt(n)
    return v.mean() - h, v.mean() + h

def hl_ci(v, alpha=0.05):
    n = len(v)
    w = np.add.outer(v, v)[np.triu_indices(n)] / 2
    w.sort()
    mu = n * (n + 1) / 4
    sd = np.sqrt(n * (n + 1) * (2 * n + 1) / 24)
    k = max(int(np.floor(mu - stats.norm.ppf(1 - alpha / 2) * sd)), 0)
    return w[k], w[len(w) - 1 - k]

for name, v in [("원자료", d), ("이상치 1개", d2)]:
    lo, hi = t_ci(v)
    l2, h2 = hl_ci(v)
    w = np.add.outer(v, v)[np.triu_indices(len(v))] / 2
    print(f"{name}")
    print(f"  평균 {v.mean():6.3f}   t  ({lo:6.3f}, {hi:6.3f})  폭 {hi - lo:6.3f}")
    print(f"  HL   {np.median(w):6.3f}   HL ({l2:6.3f}, {h2:6.3f})  폭 {h2 - l2:6.3f}")

원자료
  평균  3.183   t  ( 2.584,  3.783)  폭  1.199
  HL    3.230   HL ( 2.615,  3.820)  폭  1.205
이상치 1개
  평균  4.278   t  ( 1.917,  6.638)  폭  4.720
  HL    3.340   HL ( 2.625,  4.130)  폭  1.505

import numpy as np
from scipy import stats

rng = np.random.default_rng(21)
n, M, rho = 20, 10_000, 0.7
L = np.linalg.cholesky(np.array([[1, rho], [rho, 1]]))
hit_cc = hit_br = 0
w_cc = w_br = 0.0

for _ in range(M):
    z = rng.standard_normal((n, 2)) @ L.T
    x, y = 10 + z[:, 0], 12 + z[:, 1]          # 참 차이 -2
    miss = rng.random(n) < 0.25                # y가 완전 무작위로 25% 결측

    d = (x - y)[~miss]                         # 완전사례 대응 분석
    m = len(d)
    h = stats.t.ppf(0.975, m - 1) * d.std(ddof=1) / np.sqrt(m)
    hit_cc += (d.mean() - h <= -2 <= d.mean() + h)
    w_cc += 2 * h

    yy = y[~miss]                              # 짝을 깨고 독립표본으로
    n1, n2 = len(x), len(yy)
    sp = np.sqrt(((n1 - 1) * x.var(ddof=1) + (n2 - 1) * yy.var(ddof=1))
                 / (n1 + n2 - 2))
    hi = stats.t.ppf(0.975, n1 + n2 - 2) * sp * np.sqrt(1 / n1 + 1 / n2)
    dd = x.mean() - yy.mean()
    hit_br += (dd - hi <= -2 <= dd + hi)
    w_br += 2 * hi

print(f"완전사례 대응   포함확률 {hit_cc / M:.4f}  평균 폭 {w_cc / M:.4f}")
print(f"짝 깨고 독립    포함확률 {hit_br / M:.4f}  평균 폭 {w_br / M:.4f}")

완전사례 대응   포함확률 0.9494  평균 폭 0.8498
짝 깨고 독립    포함확률 0.9968  평균 폭 1.3792

import numpy as np
from scipy import stats

rng = np.random.default_rng(3)
n, M, rho, tau = 30, 5_000, 0.7, 1.0
res = {"변화점수": [0, 0.0], "사후만": [0, 0.0], "ANCOVA": [0, 0.0]}

def two_sample(v, g):
    a, b = v[g == 0], v[g == 1]
    sp = np.sqrt(((len(a) - 1) * a.var(ddof=1) + (len(b) - 1) * b.var(ddof=1))
                 / (len(a) + len(b) - 2))
    h = stats.t.ppf(0.975, len(v) - 2) * sp * np.sqrt(1 / len(a) + 1 / len(b))
    return b.mean() - a.mean(), h

for _ in range(M):
    g = rng.integers(0, 2, 2 * n)                      # 무작위 배정
    pre = rng.normal(0, 1, 2 * n)
    post = rho * pre + np.sqrt(1 - rho**2) * rng.normal(0, 1, 2 * n) + tau * g

    for key, v in [("변화점수", post - pre), ("사후만", post)]:
        est, h = two_sample(v, g)
        res[key][0] += (est - h <= tau <= est + h)
        res[key][1] += 2 * h

    X = np.column_stack([np.ones(2 * n), g, pre])      # ANCOVA
    beta, *_ = np.linalg.lstsq(X, post, rcond=None)
    r = post - X @ beta
    df = 2 * n - 3
    se = np.sqrt((r @ r / df) * np.linalg.inv(X.T @ X)[1, 1])
    h = stats.t.ppf(0.975, df) * se
    res["ANCOVA"][0] += (beta[1] - h <= tau <= beta[1] + h)
    res["ANCOVA"][1] += 2 * h

for k, (hit, w) in res.items():
    print(f"{k:8s} 포함확률 {hit / M:.4f}   평균 폭 {w / M:.4f}")

변화점수     포함확률 0.9448   평균 폭 0.8037
사후만      포함확률 0.9488   평균 폭 1.0405
ANCOVA   포함확률 0.9442   평균 폭 0.7489