신뢰수준과 포함확률¶

#!/usr/bin/env python3
"""평균 신뢰구간을 세 방법으로 만들어 포함확률을 비교한다.

sigma 를 아는 z, sigma 자리에 s 를 꽂아 넣은 z, 그리고 t 세 가지다.
가운데 방법이 왜 명목수준에 못 미치는지가 이 모의실험의 요점이다.

Usage:
    python mean_ci_simulation.py --method t --n-sim 100 --n 10 --alpha 0.05
    python mean_ci_simulation.py --method z_known --sigma 1.0
    python mean_ci_simulation.py --method z_plugin --N 500  # with FPC
"""

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


def finite_population_correction(n: int, N: int | None) -> float:
    """유한모집단 수정 인자 sqrt((N-n)/(N-1)). N 을 주지 않으면 1.0 을 돌려준다.

        표본이 모집단의 상당 부분을 차지하면 표준오차가 공식보다 작아진다.
        그 몫을 되돌리는 인자다.
        """
    if N is None:
        return 1.0
    if N <= 1 or n >= N:
        raise ValueError("FPC requires N > 1 and n < N.")
    return float(np.sqrt((N - n) / (N - 1)))


def simulate_data(n_sim: int, n: int, mu: float, sigma: float, rng) -> np.ndarray:
    """모양 (n_sim, n) 인 배열을 돌려준다. 행 하나가 표본 하나다."""
    return rng.normal(loc=mu, scale=sigma, size=(n_sim, n))


def compute_intervals(xbar, s, n, alpha, method, sigma_known=None, N=None):
    fpc = finite_population_correction(n, N)
    if method == "z_known":
        if sigma_known is None:
            raise ValueError("z_known requires sigma_known.")
        z_star = norm.ppf(1 - alpha / 2.0)
        se = sigma_known / np.sqrt(n) * fpc
        moe = z_star * se
    elif method == "z_plugin":
        z_star = norm.ppf(1 - alpha / 2.0)
        se = (s / np.sqrt(n)) * fpc
        moe = z_star * se
    elif method == "t":
        df = n - 1
        t_star = t.ppf(1 - alpha / 2.0, df=df)
        se = (s / np.sqrt(n)) * fpc
        moe = t_star * se
    else:
        raise ValueError(f"Unknown method: {method}")
    return xbar - moe, xbar + moe


def plot_intervals(ax, lower, upper, xbar, covered, mu, title):
    for i in range(len(xbar)):
        color = "k" if covered[i] else "r"
        ax.plot([lower[i], upper[i]], [i, i], lw=2, color=color)
        ax.plot(xbar[i], i, marker="o", ms=3, color=color)
    ax.axvline(mu, linestyle="--", linewidth=1.5, color="r")
    ax.set_title(title, fontsize=12)
    ax.set_yticks([])
    for sp in ["left", "right", "top"]:
        ax.spines[sp].set_visible(False)
    ax.set_xlabel("Mean value")


def main():
    p = argparse.ArgumentParser()
    p.add_argument("--method", choices=["z_known", "z_plugin", "t"], default="t")
    p.add_argument("--rng-seed", type=int, default=42)   # 아래 그림을 재현하려면 고정한다
    p.add_argument("--n-sim", type=int, default=100)
    p.add_argument("--n", type=int, default=10)
    p.add_argument("--mu", type=float, default=0.0)
    p.add_argument("--sigma", type=float, default=1.0)
    p.add_argument("--alpha", type=float, default=0.05)
    p.add_argument("--N", type=int, default=None)
    args, _ = p.parse_known_args()

    rng = np.random.default_rng(args.rng_seed)
    X = simulate_data(args.n_sim, args.n, args.mu, args.sigma, rng)
    xbar = X.mean(axis=1)
    s = X.std(axis=1, ddof=1)

    lower, upper = compute_intervals(
        xbar, s, args.n, args.alpha, args.method,
        sigma_known=args.sigma if args.method == "z_known" else None, N=args.N
    )
    covered = (lower <= args.mu) & (args.mu <= upper)
    n_fail = int((~covered).sum())
    coverage_pct = 100.0 * covered.mean()

    fig, ax = plt.subplots(figsize=(12, 12))
    title = (f"{args.n_sim} {args.method} CIs | n={args.n}, "
             f"CL={int((1 - args.alpha) * 100)}% | "
             f"Fail={n_fail} (Coverage ≈ {coverage_pct:.1f}%)")
    plot_intervals(ax, lower, upper, xbar, covered, args.mu, title)
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""비율 신뢰구간을 네 방법으로 만들어 포함확률을 비교한다.

Wald 는 교과서에 가장 먼저 나오지만 실제 포함확률이 가장 나쁘다.
Wilson 과 Agresti-Coull 이 그 대안이고, Clopper-Pearson 은 보수적이다.

Usage:
    python proportion_ci_simulation.py  # defaults to Wald
"""

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

rng_seed = 42        # 아래 그림을 재현하려면 고정한다
n_simulations = 100
n = 20
p_true = 0.20
alpha = 0.05
method = "wald"  # 'wald' | 'wilson' | 'ac' | 'cp'


def main():
    if rng_seed is not None:
        np.random.seed(rng_seed)

    k = np.random.binomial(n=n, p=p_true, size=n_simulations)
    phat = k / n
    lower = np.empty(n_simulations)
    upper = np.empty(n_simulations)
    z = norm.ppf(1 - alpha / 2.0)

    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
        elif method == "cp":
            lo = 0.0 if ki == 0 else beta.ppf(alpha / 2.0, ki, n - ki + 1)
            hi = 1.0 if ki == n else beta.ppf(1 - alpha / 2.0, ki + 1, n - ki)
        lower[i] = max(0.0, lo)
        upper[i] = min(1.0, hi)

    covered = (lower <= p_true) & (p_true <= upper)
    n_fail = int((~covered).sum())
    coverage_pct = 100.0 * covered.mean()

    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()} Proportion CIs | n={n}, p={p_true:.3f}, "
        f"CL={int((1 - alpha) * 100)}% | Fail={n_fail} (Coverage ≈ {coverage_pct:.1f}%)")
    ax.set_yticks([])
    for sp in ["left", "right", "top"]:
        ax.spines[sp].set_visible(False)
    ax.set_xlabel("Proportion value")
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

import numpy as np
from scipy import stats

def coverage(n, p, method, alpha=0.05):
    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 == "ac":                      # 애그레스티-콜: 성공·실패에 2씩 더함
        nt = n + z**2
        pt = (k + z**2 / 2) / nt
        se = np.sqrt(pt * (1 - pt) / nt)
        lo, hi = pt - z * se, pt + z * se
    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))
    inside = (lo <= p) & (p <= hi)
    return stats.binom.pmf(k, n, p)[inside].sum()

for m in ["wald", "wilson", "ac", "cp"]:
    print(f"{m:7s} {coverage(40, 0.1, m):.4f}")

wald    0.9145
wilson  0.9433
ac      0.9581
cp      0.9697

import numpy as np
from scipy import stats

n, z = 40, stats.norm.ppf(0.975)
k = np.arange(n + 1)
ph = k / n
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

grid = np.linspace(0.05, 0.30, 6)
for p in grid:
    inside = (lo <= p) & (p <= hi)
    print(f"p={p:.2f}  포함확률 {stats.binom.pmf(k, n, p)[inside].sum():.4f}")

p=0.05  포함확률 0.9520
p=0.10  포함확률 0.9433
p=0.15  포함확률 0.9580
p=0.20  포함확률 0.9283
p=0.25  포함확률 0.9577
p=0.30  포함확률 0.9443

import numpy as np

for M in [1_000, 10_000, 100_000]:
    se = np.sqrt(0.95 * 0.05 / M)
    print(f"M={M:>7,d}  SE={se:.5f}  ±{1.96 * se:.4f}")

M=  1,000  SE=0.00689  ±0.0135
M= 10,000  SE=0.00218  ±0.0043
M=100,000  SE=0.00069  ±0.0014

import numpy as np
from scipy import stats

rng = np.random.default_rng(20250908)
n, M = 25, 20_000

def dists():
    yield "정규", lambda size: rng.normal(0, 1, size), 0.0, 1.0
    yield "균등", lambda size: rng.uniform(-np.sqrt(3), np.sqrt(3), size), 0.0, 1.0
    yield "지수", lambda size: rng.exponential(1, size) - 1, 0.0, 1.0
    yield "t(5)", lambda size: rng.standard_t(5, size) / np.sqrt(5 / 3), 0.0, 1.0

print(f"{'분포':6s} {'초과첨도':>8s} {'t구간':>8s} {'χ²구간':>8s}")
for name, gen, mu, var in dists():
    x = gen((M, n))
    xb, s2 = x.mean(1), x.var(1, ddof=1)
    # 평균에 대한 t 구간
    h = stats.t.ppf(0.975, n - 1) * np.sqrt(s2 / n)
    cov_t = np.mean((xb - h <= mu) & (mu <= xb + h))
    # 분산에 대한 카이제곱 구간
    lo = (n - 1) * s2 / stats.chi2.ppf(0.975, n - 1)
    hi = (n - 1) * s2 / stats.chi2.ppf(0.025, n - 1)
    cov_c = np.mean((lo <= var) & (var <= hi))
    kurt = stats.kurtosis(x.ravel())
    print(f"{name:6s} {kurt:8.2f} {cov_t:8.3f} {cov_c:8.3f}")

분포         초과첨도      t구간     χ²구간
정규         0.01    0.951    0.951
균등        -1.20    0.949    0.997
지수         5.95    0.923    0.720
t(5)       4.73    0.953    0.834