비율 신뢰구간의 포함확률 모의실험¶

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

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

n_simulations = 100
n = 20
p_true = 0.20
alpha = 0.05
method = "wilson"  # 'wald' | 'wilson' | 'ac' | 'cp'

# 표본을 만들 필요가 없다. 필요한 것은 성공 횟수 k 하나뿐이므로
# 0/1을 n개 뽑는 대신 이항분포에서 k를 바로 뽑는다.
k = np.random.binomial(n=n, p=p_true, size=n_simulations)
phat = k / n
z = norm.ppf(1 - alpha / 2)

lower = np.empty(n_simulations)
upper = np.empty(n_simulations)

for i, ki in enumerate(k):
    p = ki / n
    if method == "wald":
        se = np.sqrt(p * (1 - p) / n)
        lo, hi = p - z * se, p + z * se
    elif method == "wilson":
        denom = 1 + z**2 / n
        center = (p + z**2 / (2 * n)) / denom
        half = z * np.sqrt(p * (1 - p) / n + z**2 / (4 * n**2)) / denom
        lo, hi = center - half, center + half
    elif method == "ac":
        n_tilde = n + z**2
        p_tilde = (ki + 0.5 * z**2) / n_tilde
        se_tilde = np.sqrt(p_tilde * (1 - p_tilde) / n_tilde)
        lo, hi = p_tilde - z * se_tilde, p_tilde + z * se_tilde
    else:  # cp
        lo = 0.0 if ki == 0 else beta.ppf(alpha / 2, ki, n - ki + 1)
        hi = 1.0 if ki == n else beta.ppf(1 - alpha / 2, ki + 1, n - ki)

    lower[i] = max(0.0, lo)
    upper[i] = min(1.0, hi)

covered = (lower <= p_true) & (p_true <= upper)
coverage_pct = 100.0 * covered.mean()
print(f"{method} coverage: {coverage_pct:.1f}%")


# (1) 포함률은 유한합이다. k 가 21가지뿐이므로 담는 k 를 모아 더하면 끝이다.
from scipy.stats import binom


def interval(ki, which):
    """성공 횟수 ki 하나에서 네 방법의 구간 끝점을 돌려준다."""
    p = ki / n
    if which == "wald":
        se = np.sqrt(p * (1 - p) / n)
        lo, hi = p - z * se, p + z * se
    elif which == "wilson":
        denom = 1 + z**2 / n
        center = (p + z**2 / (2 * n)) / denom
        half = z * np.sqrt(p * (1 - p) / n + z**2 / (4 * n**2)) / denom
        lo, hi = center - half, center + half
    elif which == "ac":
        n_tilde = n + z**2
        p_tilde = (ki + 0.5 * z**2) / n_tilde
        se_tilde = np.sqrt(p_tilde * (1 - p_tilde) / n_tilde)
        lo, hi = p_tilde - z * se_tilde, p_tilde + z * se_tilde
    else:
        lo = 0.0 if ki == 0 else beta.ppf(alpha / 2, ki, n - ki + 1)
        hi = 1.0 if ki == n else beta.ppf(1 - alpha / 2, ki + 1, n - ki)
    return max(0.0, lo), min(1.0, hi)


def covers(ki, which):
    lo, hi = interval(ki, which)
    return lo <= p_true <= hi


print()
print(f"{'방법':<8}{'담는 k':<12}{'정확':>9}{'모의':>8}{'차이/MC오차':>10}")
for m in ("wald", "wilson", "ac", "cp"):
    S = [ki for ki in range(n + 1) if covers(ki, m)]
    exact = sum(binom.pmf(ki, n, p_true) for ki in S)
    sim = np.mean([covers(ki, m) for ki in k])
    se_mc = np.sqrt(exact * (1 - exact) / n_simulations)
    print(f"{m:<8}{f'{S[0]}..{S[-1]}':<12}{exact:>9.5f}{sim:>8.2f}{abs(sim - exact) / se_mc:>10.2f}")

wilson coverage: 96.0%

방법      담는 k               정확      모의   차이/MC오차
wald    2..8          0.92084    0.91      0.40
wilson  1..7          0.95633    0.96      0.18
ac      1..7          0.95633    0.96      0.18
cp      1..8          0.97849    0.99      0.79

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

ax.axvline(p_true, linestyle="--", linewidth=1.5)
ax.set_title(f"{n_simulations} {method.upper()} CIs | n={n}, p={p_true}, CL=95%")
ax.set_yticks([])
ax.set_xlabel("Proportion value")
plt.tight_layout()
plt.show()

import numpy as np
from scipy.stats import norm, binom

n, p, alpha = 20, 0.05, 0.05
z = norm.ppf(1 - alpha / 2)
coverage = 0.0
for k in range(n + 1):
    phat = k / n
    se = np.sqrt(phat * (1 - phat) / n)
    lo, hi = max(0.0, phat - z * se), min(1.0, phat + z * se)
    if lo <= p <= hi:
        coverage += binom.pmf(k, n, p)
print(f"exact Wald coverage: {100 * coverage:.1f}%")

exact Wald coverage: 63.9%

import numpy as np
from scipy import stats

n, alpha = 40, 0.05
z = stats.norm.ppf(1 - alpha / 2)
k = np.arange(n + 1)
ph = k / n

def bounds(m):
    if m == "wald":
        se = np.sqrt(ph * (1 - ph) / n)
        return ph - z * se, ph + z * se
    if m == "wilson":
        d = 1 + z**2 / n
        c = (ph + z**2 / (2 * n)) / d
        h = z / d * np.sqrt(ph * (1 - ph) / n + z**2 / (4 * n**2))
        return c - h, c + h
    if m == "ac":
        nt = n + z**2
        pt = (k + z**2 / 2) / nt
        se = np.sqrt(pt * (1 - pt) / nt)
        return pt - z * se, pt + z * se
    if m == "jeffreys":
        return (stats.beta.ppf(alpha / 2, k + 0.5, n - k + 0.5),
                stats.beta.ppf(1 - alpha / 2, k + 0.5, n - k + 0.5))
    if m == "cp":
        return (np.where(k == 0, 0.0, stats.beta.ppf(alpha / 2, k, n - k + 1)),
                np.where(k == n, 1.0, stats.beta.ppf(1 - alpha / 2, k + 1, n - k)))

methods = ["wald", "wilson", "ac", "jeffreys", "cp"]
names = ["왈드", "윌슨", "AC", "제프리스", "CP"]
for p in [0.05, 0.15, 0.30, 0.50]:
    w = stats.binom.pmf(k, n, p)
    out = []
    for m in methods:
        lo, hi = bounds(m)
        cov = w[(lo <= p) & (p <= hi)].sum()
        wid = ((hi - lo) * w).sum()
        out.append(f"{cov:.3f}/{wid:.3f}")
    print(f"p={p:.2f}  " + "  ".join(f"{nm} {o}" for nm, o in zip(names, out)))

p=0.05  왈드 0.868/0.121  윌슨 0.952/0.146  AC 0.986/0.166  제프리스 0.986/0.133  CP 0.986/0.157
p=0.15  왈드 0.939/0.215  윌슨 0.958/0.216  AC 0.958/0.224  제프리스 0.958/0.213  CP 0.976/0.236
p=0.30  왈드 0.930/0.280  윌슨 0.944/0.270  AC 0.944/0.272  제프리스 0.944/0.273  CP 0.961/0.296
p=0.50  왈드 0.919/0.306  윌슨 0.962/0.293  AC 0.962/0.293  제프리스 0.962/0.297  CP 0.962/0.320

import numpy as np
from scipy import stats

def min_cov(n, method, alpha=0.05, grid=2001):
    z = stats.norm.ppf(1 - alpha / 2)
    k = np.arange(n + 1)
    ph = k / n
    if method == "wald":
        se = np.sqrt(ph * (1 - ph) / n)
        lo, hi = ph - z * se, ph + z * se
    elif method == "wilson":
        d = 1 + z**2 / n
        c = (ph + z**2 / (2 * n)) / d
        h = z / d * np.sqrt(ph * (1 - ph) / n + z**2 / (4 * n**2))
        lo, hi = c - h, c + h
    elif method == "jeffreys":
        lo = stats.beta.ppf(alpha / 2, k + 0.5, n - k + 0.5)
        hi = stats.beta.ppf(1 - alpha / 2, k + 0.5, n - k + 0.5)
    elif method == "cp":
        lo = np.where(k == 0, 0.0, stats.beta.ppf(alpha / 2, k, n - k + 1))
        hi = np.where(k == n, 1.0, stats.beta.ppf(1 - alpha / 2, k + 1, n - k))
    ps = np.linspace(0.001, 0.999, grid)
    pmf = stats.binom.pmf(k[:, None], n, ps[None, :])
    inside = (lo[:, None] <= ps[None, :]) & (ps[None, :] <= hi[:, None])
    return (pmf * inside).sum(axis=0).min()

print(f"{'n':>5s} {'왈드':>8s} {'윌슨':>8s} {'제프리스':>9s} {'CP':>8s}")
for n in [10, 25, 50, 100, 200]:
    vals = [min_cov(n, m) for m in ["wald", "wilson", "jeffreys", "cp"]]
    print(f"{n:5d} " + " ".join(f"{v:8.4f}" for v in vals))

    n       왈드       윌슨      제프리스       CP
   10   0.0100   0.8384   0.8682   0.9611
   25   0.0247   0.8392   0.8906   0.9505
   50   0.0488   0.8394   0.8840   0.9509
  100   0.0952   0.8607   0.8805   0.9503
  200   0.1813   0.9102   0.8781   0.9504

import numpy as np
from scipy import stats

n, p_true, alpha = 1000, 0.002, 0.05
k = np.arange(0, 25)

# 포아송 정확 구간 (λ 척도 → p = λ/n)
lam_lo = np.where(k == 0, 0.0, stats.chi2.ppf(alpha / 2, 2 * k) / 2)
lam_hi = stats.chi2.ppf(1 - alpha / 2, 2 * k + 2) / 2
po_lo, po_hi = lam_lo / n, lam_hi / n

# 클로퍼-피어슨
cp_lo = np.where(k == 0, 0.0, stats.beta.ppf(alpha / 2, k, n - k + 1))
cp_hi = stats.beta.ppf(1 - alpha / 2, k + 1, n - k)

print(f"{'k':>3s} {'포아송 구간':>24s} {'클로퍼-피어슨':>24s}")
for i in [0, 1, 2, 5]:
    print(f"{i:3d}  ({po_lo[i]:.6f}, {po_hi[i]:.6f})   "
          f"({cp_lo[i]:.6f}, {cp_hi[i]:.6f})")

kk = np.arange(n + 1)
for name, lo, hi in [("포아송", po_lo, po_hi), ("CP", cp_lo, cp_hi)]:
    full_lo = np.concatenate([lo, np.full(n + 1 - len(lo), 1.0)])
    full_hi = np.concatenate([hi, np.full(n + 1 - len(hi), 1.0)])
    inside = (full_lo <= p_true) & (p_true <= full_hi)
    print(f"{name} 포함확률 {stats.binom.pmf(kk, n, p_true)[inside].sum():.4f}")

  k                   포아송 구간                  클로퍼-피어슨
  0  (0.000000, 0.003689)   (0.000000, 0.003682)
  1  (0.000025, 0.005572)   (0.000025, 0.005559)
  2  (0.000242, 0.007225)   (0.000242, 0.007206)
  5  (0.001623, 0.011668)   (0.001625, 0.011629)
포아송 포함확률 0.9835
CP 포함확률 0.9835

import numpy as np
from scipy import stats

rng = np.random.default_rng(77)
n, M, B = 30, 4_000, 999
z = stats.norm.ppf(0.975)

print(f"{'p':>5s} {'백분위':>8s} {'왈드':>8s} {'윌슨':>8s} {'퇴화':>7s}")
for p_true in [0.03, 0.05, 0.10, 0.20, 0.40]:
    hb = hw = hs = degenerate = 0
    for _ in range(M):
        k = rng.binomial(n, p_true)
        ph = k / n
        star = rng.binomial(n, ph, B) / n          # 비모수 부트스트랩과 동치
        lo, hi = np.percentile(star, [2.5, 97.5])
        hb += (lo <= p_true <= hi)
        degenerate += (hi - lo == 0)
        se = np.sqrt(ph * (1 - ph) / n)
        hw += (ph - z * se <= p_true <= ph + z * se)
        d = 1 + z**2 / n
        c = (ph + z**2 / (2 * n)) / d
        h = z / d * np.sqrt(ph * (1 - ph) / n + z**2 / (4 * n**2))
        hs += (c - h <= p_true <= c + h)
    print(f"{p_true:5.2f} {hb / M:8.4f} {hw / M:8.4f} "
          f"{hs / M:8.4f} {degenerate / M:7.4f}")

    p      백분위       왈드       윌슨      퇴화
 0.03   0.5965   0.6062   0.9345  0.3927
 0.05   0.7792   0.7823   0.9433  0.2152
 0.10   0.9555   0.8100   0.9742  0.0395
 0.20   0.9440   0.9440   0.9627  0.0013
 0.40   0.9547   0.9265   0.9560  0.0000

import numpy as np
from scipy import stats

def exact_coverage(n, p, alpha=0.05):
    z = stats.norm.ppf(1 - alpha / 2)
    k = np.arange(n + 1)
    ph = k / n
    se = np.sqrt(ph * (1 - ph) / n)
    inside = (ph - z * se <= p) & (p <= ph + z * se)
    return stats.binom.pmf(k, n, p)[inside].sum()

def mc_coverage(n, p, M, seed, alpha=0.05):
    rng = np.random.default_rng(seed)
    z = stats.norm.ppf(1 - alpha / 2)
    ph = rng.binomial(n, p, M) / n
    se = np.sqrt(ph * (1 - ph) / n)
    return np.mean((ph - z * se <= p) & (p <= ph + z * se))

n, p = 20, 0.05
exact = exact_coverage(n, p)
print(f"정확 계산      {exact:.6f}")
for M in [1_000, 10_000, 100_000]:
    est = mc_coverage(n, p, M, 1)
    print(f"모의실험 M={M:>7,d}  {est:.6f}   오차 {est - exact:+.6f}   "
          f"MCSE {np.sqrt(exact * (1 - exact) / M):.6f}")

정확 계산      0.638940
모의실험 M=  1,000  0.636000   오차 -0.002940   MCSE 0.015189
모의실험 M= 10,000  0.643000   오차 +0.004060   MCSE 0.004803
모의실험 M=100,000  0.637680   오차 -0.001260   MCSE 0.001519