모분산의 신뢰구간¶

import numpy as np
from scipy import stats

n = 25
s_squared = 120
alpha = 0.05
df = n - 1

# 카이제곱은 좌우가 대칭이 아니므로 두 기각값을 따로 구한다.
chi2_lower = stats.chi2.ppf(alpha / 2, df)
chi2_upper = stats.chi2.ppf(1 - alpha / 2, df)

# 위아래가 뒤집혀 들어간다. (n-1)S^2/sigma^2 이 카이제곱을 따르므로
# sigma^2 로 풀면 분모에 기각값이 오고, 큰 기각값이 구간의 아래끝을 만든다.
# 부호를 헷갈리기 쉬운 대목이다.
ci_lower = df * s_squared / chi2_upper
ci_upper = df * s_squared / chi2_lower

# 표준편차의 구간은 분산 구간에 제곱근을 씌우면 된다. 제곱근이 단조함수라
# 순서가 그대로 보존되기 때문이다.
ci_sd_lower = np.sqrt(ci_lower)
ci_sd_upper = np.sqrt(ci_upper)

print(f"95% CI for variance: ({ci_lower:.2f}, {ci_upper:.2f})")
print(f"95% CI for std dev:  ({ci_sd_lower:.2f}, {ci_sd_upper:.2f})")

# (1) 양 끝이 S^2 의 상수배다. 그 상수는 자유도와 alpha 만으로 정해진다.
k_lo = df / chi2_upper
k_hi = df / chi2_lower
print(f"\n하한 배수 (n-1)/chi2_0.975 = {k_lo:.5f}")
print(f"상한 배수 (n-1)/chi2_0.025 = {k_hi:.5f}")
print(f"분산 구간의 상한/하한 비 = {k_hi / k_lo:.5f}  (= chi2_0.975/chi2_0.025 = {chi2_upper / chi2_lower:.5f})")
print(f"표준편차 구간의 비        = {np.sqrt(k_hi / k_lo):.5f}")

# S^2 을 바꾸어도 비가 그대로인가.
print(f"\n{'S^2':>10}{'하한':>12}{'상한':>12}{'비':>10}{'sd 비':>9}")
for s2 in (0.05, 7.0, 120.0, 9_999.0):
    lo, hi = df * s2 / chi2_upper, df * s2 / chi2_lower
    print(f"{s2:>10.2f}{lo:>12.5f}{hi:>12.5f}{hi / lo:>10.5f}{np.sqrt(hi / lo):>9.5f}")

# (2) 구간은 S^2 을 산술적으로도 기하적으로도 중심에 두지 않는다.
print(f"\n점추정값 S^2 = {s_squared}")
print(f"  아래로 {s_squared - ci_lower:.2f},  위로 {ci_upper - s_squared:.2f}  (위쪽이 {(ci_upper - s_squared) / (s_squared - ci_lower):.3f} 배)")
print(f"  산술 중점 = {(ci_lower + ci_upper) / 2:.2f}")
print(f"  기하 중점 = {np.sqrt(ci_lower * ci_upper):.2f}   (배수 {np.sqrt(k_lo * k_hi):.5f})")
print(f"  두 분위수의 기하평균 = {np.sqrt(chi2_lower * chi2_upper):.4f}  vs  n-1 = {df}")

95% CI for variance: (73.16, 232.24)
95% CI for std dev:  (8.55, 15.24)

하한 배수 (n-1)/chi2_0.975 = 0.60969
상한 배수 (n-1)/chi2_0.025 = 1.93530
분산 구간의 상한/하한 비 = 3.17423  (= chi2_0.975/chi2_0.025 = 3.17423)
표준편차 구간의 비        = 1.78164

       S^2          하한          상한         비     sd 비
      0.05     0.03048     0.09677   3.17423  1.78164
      7.00     4.26785    13.54713   3.17423  1.78164
    120.00    73.16315   232.23652   3.17423  1.78164
   9999.00  6096.31974 19351.10823   3.17423  1.78164

점추정값 S^2 = 120
  아래로 46.84,  위로 112.24  (위쪽이 2.396 배)
  산술 중점 = 152.70
  기하 중점 = 130.35   (배수 1.08625)
  두 분위수의 기하평균 = 22.0943  vs  n-1 = 24

from scipy import stats

print(f"{'n':>6} {'var ratio':>10} {'sd ratio':>10}")
for n in [10, 25, 50, 100, 500, 1000]:
    df = n - 1
    r = stats.chi2.ppf(0.975, df) / stats.chi2.ppf(0.025, df)
    print(f"{n:>6} {r:>10.3f} {r**0.5:>10.3f}")

     n  var ratio   sd ratio
    10      7.044      2.654
    25      3.174      1.782
    50      2.225      1.492
   100      1.751      1.323
   500      1.282      1.132
  1000      1.192      1.092

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n, df, R = 25, 24, 20000
lo_q = stats.chi2.ppf(0.025, df)
hi_q = stats.chi2.ppf(0.975, df)

cases = [
    ("Normal",      lambda: rng.normal(0, 1, n),      1.0),
    ("t(5)",        lambda: rng.standard_t(5, n),     5 / 3),
    ("Uniform",     lambda: rng.uniform(0, 1, n),     1 / 12),
    ("Exponential", lambda: rng.exponential(1, n),    1.0),
]

for name, gen, true_var in cases:
    hits = 0
    for _ in range(R):
        s2 = gen().var(ddof=1)
        lo, hi = df * s2 / hi_q, df * s2 / lo_q
        hits += (lo <= true_var <= hi)
    print(f"{name:>12}: coverage = {hits / R:.4f}")

      Normal: coverage = 0.9504
        t(5): coverage = 0.8292
     Uniform: coverage = 0.9961
 Exponential: coverage = 0.7218

from scipy import stats
import numpy as np

n, s = 30, 0.42
ub_var = (n - 1) * s**2 / stats.chi2.ppf(0.05, n - 1)
print(f"95% upper bound: sigma^2 <= {ub_var:.4f}, sigma <= {np.sqrt(ub_var):.4f}")

95% upper bound: sigma^2 <= 0.2889, sigma <= 0.5375

from scipy import stats

n, s2, df = 25, 0.05, 24
lo = df * s2 / stats.chi2.ppf(0.975, df)
hi = df * s2 / stats.chi2.ppf(0.025, df)
print(f"95% CI: ({lo:.5f}, {hi:.5f})\n")

for sigma0_sq in [0.025, 0.0305, 0.04, 0.07, 0.0968, 0.12]:
    c = df * s2 / sigma0_sq
    p = 2 * min(stats.chi2.cdf(c, df), stats.chi2.sf(c, df))
    inside = lo <= sigma0_sq <= hi
    print(f"sigma0^2 = {sigma0_sq:.4f}: chi2 = {c:7.3f}, "
          f"p = {p:.4f}, in CI = {inside}")

95% CI: (0.03048, 0.09677)

sigma0^2 = 0.0250: chi2 =  48.000, p = 0.0050, in CI = False
sigma0^2 = 0.0305: chi2 =  39.344, p = 0.0502, in CI = True
sigma0^2 = 0.0400: chi2 =  30.000, p = 0.3695, in CI = True
sigma0^2 = 0.0700: chi2 =  17.143, p = 0.3150, in CI = True
sigma0^2 = 0.0968: chi2 =  12.397, p = 0.0499, in CI = False
sigma0^2 = 0.1200: chi2 =  10.000, p = 0.0109, in CI = False