σ₁² / σ₂²의 신뢰구간¶

import numpy as np
from scipy.stats import f

n1, n2 = 15, 12
alpha = 0.05

# 참 분산비는 1/1.5² = 0.444다. 답을 알고 있는 상태에서 구간을 본다.
rng = np.random.default_rng(42)
x = rng.normal(loc=0, scale=1.0, size=n1)
y = rng.normal(loc=0, scale=1.5, size=n2)

s1_sq = x.var(ddof=1)
s2_sq = y.var(ddof=1)
rhat = s1_sq / s2_sq

# F분포는 두 자유도의 **순서**가 중요하다. dfn이 분자, dfd가 분모다.
# 두 표본의 크기가 다르므로 여기서 뒤바꾸면 조용히 틀린 답이 나온다.
df1, df2 = n1 - 1, n2 - 1
F_lo = f(dfn=df1, dfd=df2).ppf(alpha / 2.0)
F_hi = f(dfn=df1, dfd=df2).ppf(1 - alpha / 2.0)

# 여기서도 큰 임계값이 아래끝의 분모로 간다.
# (s1²/s2²)/(σ1²/σ2²) ~ F 를 σ1²/σ2² 에 대해 풀면 대소가 뒤집히기 때문이다.
ci_lower = rhat / F_hi
ci_upper = rhat / F_lo

print(f"95% CI for σ₁²/σ₂²: ({ci_lower:.4f}, {ci_upper:.4f})")

95% CI for σ₁²/σ₂²: (0.2396, 2.4903)

#!/usr/bin/env python3
"""두 분산의 비에 대한 F 신뢰구간을 100번 만들어 포함확률을 확인한다.

두 표본분산의 비가 F 분포를 따른다는 사실에서 구간이 나온다. 이 방법은
정규성에 특히 민감해서, 모집단이 정규가 아니면 포함확률이 크게 어긋난다.
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import f

rng_seed = 42        # 아래 그림을 재현하려면 고정한다
n_simulations = 100
n1, n2 = 15, 12
mu1, mu2 = 0.0, 0.0
sigma1, sigma2 = 1.0, 1.5
alpha = 0.05


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

    theta_true = (sigma1**2) / (sigma2**2)
    df1, df2 = n1 - 1, n2 - 1

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

    F_lo = f(dfn=df1, dfd=df2).ppf(alpha / 2.0)
    F_hi = f(dfn=df1, dfd=df2).ppf(1 - alpha / 2.0)

    for i in range(n_simulations):
        x = np.random.normal(loc=mu1, scale=sigma1, size=n1)
        y = np.random.normal(loc=mu2, scale=sigma2, size=n2)
        s1_sq = x.var(ddof=1)
        s2_sq = y.var(ddof=1)
        rhat = s1_sq / s2_sq
        lowers[i] = rhat / F_hi
        uppers[i] = rhat / F_lo
        centers[i] = rhat

    covered = (lowers <= theta_true) & (theta_true <= uppers)
    n_fail = (~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([lowers[i], uppers[i]], [i, i], lw=2, color=color)
        ax.plot(centers[i], i, marker="o", ms=3, color=color)

    ax.axvline(theta_true, linestyle="--", linewidth=1.5, color="r")
    ax.set_title(
        f"{n_simulations} F-intervals for σ₁²/σ₂² | n1={n1}, n2={n2}, "
        f"df=({df1},{df2}), CL={int((1 - alpha) * 100)}% | "
        f"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("θ = σ₁² / σ₂²")
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

import numpy as np
from scipy import stats

rng = np.random.default_rng(23)
M = 20_000
cases = [("정규", lambda s: rng.normal(0, 1, s)),
         ("균등", lambda s: rng.uniform(-1.732, 1.732, s)),
         ("t(10)", lambda s: rng.standard_t(10, s)),
         ("t(5)", lambda s: rng.standard_t(5, s)),
         ("지수", lambda s: rng.exponential(1, s))]
print(f"{'분포':>7s} {'n=20':>9s} {'n=50':>9s}")
for name, gen in cases:
    row = []
    for n in [20, 50]:
        x, y = gen((M, n)), gen((M, n))
        r = x.var(1, ddof=1) / y.var(1, ddof=1)
        lo = r / stats.f.ppf(0.975, n - 1, n - 1)
        hi = r / stats.f.ppf(0.025, n - 1, n - 1)
        row.append(np.mean((lo <= 1) & (1 <= hi)))
    print(f"{name:>7s} {row[0]:9.4f} {row[1]:9.4f}")

     분포      n=20      n=50
     정규    0.9473    0.9504
     균등    0.9944    0.9967
  t(10)    0.9075    0.9000
   t(5)    0.8373    0.8025
     지수    0.7323    0.7046

import numpy as np
from scipy import stats

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

cases = [("정규", lambda s: rng.normal(0, 1, s)),
         ("t(5)", lambda s: rng.standard_t(5, s)),
         ("지수", lambda s: rng.exponential(1, s))]
for name, gen in cases:
    hf = hb = hk = 0
    for _ in range(M):
        x, y = gen(n1), gen(n2)
        r = x.var(ddof=1) / y.var(ddof=1)

        lo = r / stats.f.ppf(0.975, n1 - 1, n2 - 1)
        hi = r / stats.f.ppf(0.025, n1 - 1, n2 - 1)
        hf += (lo <= 1 <= hi)

        bx = x[rng.integers(0, n1, (B, n1))]
        by = y[rng.integers(0, n2, (B, n2))]
        br = np.log(bx.var(1, ddof=1) / by.var(1, ddof=1))
        l, h = np.percentile(br, [2.5, 97.5])
        hb += (l <= 0 <= h)

        kk = (stats.kurtosis(x, bias=False) + stats.kurtosis(y, bias=False)) / 2
        sd = np.sqrt((kk + 2) * (1 / (n1 - 1) + 1 / (n2 - 1)))
        hk += (abs(np.log(r)) <= z * sd)

    print(f"{name:5s} F {hf / M:.4f}   부트스트랩(로그) {hb / M:.4f}   "
          f"첨도 보정 {hk / M:.4f}")

정규    F 0.9485   부트스트랩(로그) 0.9340   첨도 보정 0.9400
t(5)  F 0.8100   부트스트랩(로그) 0.9110   첨도 보정 0.9015
지수    F 0.7065   부트스트랩(로그) 0.9040   첨도 보정 0.8710

import numpy as np
from scipy import stats

s1sq, n1 = 25.0, 16
s2sq, n2 = 10.0, 21
r = s1sq / s2sq
lo = r / stats.f.ppf(0.975, n1 - 1, n2 - 1)
hi = r / stats.f.ppf(0.025, n1 - 1, n2 - 1)
print(f"분산비 {r:.4f}      구간 ({lo:.4f}, {hi:.4f})")
print(f"표준편차비 {np.sqrt(r):.4f}  구간 ({np.sqrt(lo):.4f}, {np.sqrt(hi):.4f})")

# 순서를 바꾸면
r2 = s2sq / s1sq
lo2 = r2 / stats.f.ppf(0.975, n2 - 1, n1 - 1)
hi2 = r2 / stats.f.ppf(0.025, n2 - 1, n1 - 1)
print(f"\n역순 분산비 {r2:.4f}  구간 ({lo2:.4f}, {hi2:.4f})")
print(f"원 구간의 역수      ({1 / hi:.4f}, {1 / lo:.4f})")

분산비 2.5000      구간 (0.9716, 6.8898)
표준편차비 1.5811  구간 (0.9857, 2.6248)

역순 분산비 0.4000  구간 (0.1452, 1.0293)
원 구간의 역수      (0.1452, 1.0293)

import numpy as np
from scipy import stats

print(f"{'n(양군)':>8s} {'하한/r':>9s} {'상한/r':>9s} {'상한/하한':>10s}")
for n in [10, 20, 30, 50, 100, 200, 500]:
    lo = 1 / stats.f.ppf(0.975, n - 1, n - 1)
    hi = 1 / stats.f.ppf(0.025, n - 1, n - 1)
    print(f"{n:8d} {lo:9.4f} {hi:9.4f} {hi / lo:10.2f}")

print()
for target in [10, 4, 2, 1.5]:
    n = next(m for m in range(5, 5000)
             if (1 / stats.f.ppf(0.025, m - 1, m - 1))
             / (1 / stats.f.ppf(0.975, m - 1, m - 1)) <= target)
    print(f"상한/하한 {target}배 이하 → 집단당 n = {n}")

   n(양군)      하한/r      상한/r      상한/하한
      10    0.2484    4.0260      16.21
      20    0.3958    2.5265       6.38
      30    0.4760    2.1010       4.41
      50    0.5675    1.7622       3.11
     100    0.6728    1.4862       2.21
     200    0.7568    1.3214       1.75
     500    0.8389    1.1921       1.42

상한/하한 10배 이하 → 집단당 n = 14
상한/하한 4배 이하 → 집단당 n = 35
상한/하한 2배 이하 → 집단당 n = 131
상한/하한 1.5배 이하 → 집단당 n = 376

for lam in [0.5, 1, 2, 4, 9, 19]:
    print(f"σ_b²/σ_w² = {lam:4.1f}  →  ICC = {lam / (1 + lam):.3f}")

σ_b²/σ_w² =  0.5  →  ICC = 0.333
σ_b²/σ_w² =  1.0  →  ICC = 0.500
σ_b²/σ_w² =  2.0  →  ICC = 0.667
σ_b²/σ_w² =  4.0  →  ICC = 0.800
σ_b²/σ_w² =  9.0  →  ICC = 0.900
σ_b²/σ_w² = 19.0  →  ICC = 0.950

import numpy as np
from scipy import stats

rng = np.random.default_rng(404)
n, M = 30, 3_000
print(f"{'분포':>6s} {'F':>8s} {'바틀렛':>8s} {'레빈':>8s} {'BF':>8s}")
for name, gen in [("정규", lambda s: rng.normal(0, 1, s)),
                  ("t(5)", lambda s: rng.standard_t(5, s)),
                  ("지수", lambda s: rng.exponential(1, s))]:
    rej = np.zeros(4)
    for _ in range(M):
        x, y = gen(n), gen(n)
        f = x.var(ddof=1) / y.var(ddof=1)
        pf = 2 * min(stats.f.cdf(f, n - 1, n - 1), stats.f.sf(f, n - 1, n - 1))
        pb = stats.bartlett(x, y).pvalue
        pl = stats.levene(x, y, center="mean").pvalue
        pv = stats.levene(x, y, center="median").pvalue
        rej += np.array([pf, pb, pl, pv]) < 0.05
    print(f"{name:>6s} " + " ".join(f"{v / M:8.4f}" for v in rej))

    분포        F      바틀렛       레빈       BF
    정규   0.0527   0.0527   0.0517   0.0440
  t(5)   0.1667   0.1667   0.0470   0.0370
    지수   0.2790   0.2790   0.1330   0.0577