일표본 분산 검정¶

from scipy.stats import chi2

def test_variance_one_sample(n, s2, sigma0, alt="two-sided", alpha=0.05):
    """H0: sigma^2 = sigma0^2. **정규모집단**을 가정한다.

    이 가정은 형식적인 단서가 아니다. 평균 검정과 달리 여기서는
    중심극한정리가 도와주지 않아서 n을 키워도 비정규성이 상쇄되지 않는다.
    s2에는 ddof=1로 계산한 표본분산을 넣는다.
    """
    df = n - 1
    chi2_stat = df * s2 / (sigma0 ** 2)
    if alt == "two-sided":
        # 카이제곱분포는 비대칭이라 "양쪽 꼬리"를 나누는 방식이 여럿이다.
        # 여기서는 작은 쪽 꼬리를 두 배 하는 관례를 따랐다. 간단하고
        # 신뢰구간과 어긋나지 않지만, 확률이 반씩 나뉘지는 않는다.
        p = 2 * min(chi2.cdf(chi2_stat, df), 1 - chi2.cdf(chi2_stat, df))
    elif alt == "less":
        p = chi2.cdf(chi2_stat, df)
    else:
        p = 1 - chi2.cdf(chi2_stat, df)
    return chi2_stat, p, (p < alpha)

from scipy.optimize import brentq, minimize_scalar

nu, alpha = 11, 0.05
a = chi2.ppf(alpha / 2, nu)        # 왼쪽 임계값
b = chi2.ppf(1 - alpha / 2, nu)    # 오른쪽 임계값
print(f"등꼬리 임계값:  a = {a:.6f},  b = {b:.6f}")

# (1) 함수가 쓰는 p-값이 두 임계값에서 정확히 alpha 인가.
#     sigma0 = 1 로 두면 chi2 = (n-1) s2 이므로 s2 = x/(n-1) 을 넣으면 된다.
for x, name in ((a, "a"), (b, "b")):
    for eps in (-1e-6, 1e-6):
        _, p, rej = test_variance_one_sample(12, (x + eps) / nu, 1.0)
        print(f"  chi2 = {name}{eps:+.0e} : p = {p:.12f}  기각 {rej}")

# (2) 검정력 함수와 편향
def power(theta):
    """theta = sigma^2/sigma0^2 일 때 기각확률."""
    return chi2.cdf(a / theta, nu) + chi2.sf(b / theta, nu)

fa, fb = chi2.pdf(a, nu), chi2.pdf(b, nu)
print(f"\na f(a) = {a * fa:.6f},  b f(b) = {b * fb:.6f},  "
      f"기울기 pi'(1) = {b * fb - a * fa:+.6f}")
res = minimize_scalar(power, bounds=(0.3, 3.0), method="bounded",
                      options={"xatol": 1e-12})
left = brentq(lambda th: power(th) - alpha, 0.3, 0.9)
print(f"pi(1) = {power(1.0):.6f}   최소 검정력 {res.fun:.6f} at theta = {res.x:.6f}")
print(f"검정력이 명목 아래로 내려가는 구간 = ({left:.6f}, 1)")
print("\n  theta   검정력")
for th in (0.80, 0.85, 0.90, 0.9414, 0.97, 1.00, 1.05, 1.20):
    print(f"  {th:.4f}  {power(th):.6f}")

등꼬리 임계값:  a = 3.815748,  b = 21.920049
  chi2 = a-1e-06 : p = 0.049999948116  기각 True
  chi2 = a+1e-06 : p = 0.050000051885  기각 False
  chi2 = b-1e-06 : p = 0.050000015864  기각 False
  chi2 = b+1e-06 : p = 0.049999984136  기각 True

a f(a) = 0.098989,  b f(b) = 0.173872,  기울기 pi'(1) = +0.074882
pi(1) = 0.050000   최소 검정력 0.047778 at theta = 0.941417
검정력이 명목 아래로 내려가는 구간 = (0.884423, 1)

  theta   검정력
  0.8000  0.062191
  0.8500  0.053607
  0.9000  0.048942
  0.9414  0.047778
  0.9700  0.048314
  1.0000  0.050000
  1.0500  0.055249
  1.2000  0.087465

stat, p, reject = test_variance_one_sample(
    n=12, s2=2.1**2, sigma0=2.0, alt="greater"
)
print("chi2:", stat, "p:", p, "reject:", reject)

# 표본분산은 그대로 두고 표본크기만 키우면 어떻게 되는지 본다.
for n in [12, 50, 200]:
    st, pv, rj = test_variance_one_sample(n=n, s2=2.1**2, sigma0=2.0, alt="greater")
    print(f"n={n:>4}: chi2={st:8.2f}  p={pv:.4f}  reject={rj}")

chi2: 12.127500000000001 p: 0.35413619761553705 reject: False
n=  12: chi2=   12.13  p=0.3541  reject=False
n=  50: chi2=   54.02  p=0.2885  reject=False
n= 200: chi2=  219.40  p=0.1532  reject=False

import math

import numpy as np
from scipy.stats import chi2, norm
from scipy.optimize import brentq

r = 2.1**2 / 2.0**2
z95 = norm.ppf(0.95)

print(f"r = {r},  r - 1 = {r - 1:.4f}")
print("    n      정확 p     단순근사 p    WH 근사 p")
for n in (12, 50, 200, 538, 1000):
    nu = n - 1
    z1 = (r - 1) * math.sqrt(nu / 2)                             # (1) 의 식
    z2 = (r**(1 / 3) - 1 + 2 / (9 * nu)) * math.sqrt(9 * nu / 2)  # 윌슨-힐퍼티
    print(f"{n:5d}   {chi2.sf(nu * r, nu):.6f}    {norm.sf(z1):.6f}"
          f"    {norm.sf(z2):.6f}")

# p < 0.05 인 n 을 2 부터 20000 까지 모두 훑는다. "첫 칸" 만 보지 않는다.
ns = np.arange(2, 20_001)
pv = chi2.sf((ns - 1) * r, ns - 1)
ok = ns[pv < 0.05]
print(f"\np < 0.05 인 n 의 개수 = {len(ok)},  가장 작은 n = {ok.min()},"
      f"  그 구간이 연속인가 = {bool(np.all(np.diff(ok) == 1))}")
print(f"p 가 n 에 대해 단조감소인가 = {bool(np.all(np.diff(pv) < 0))}"
      f"   (최대 p = {pv.max():.6f} at n = {ns[pv.argmax()]})")

print(f"\n단순근사가 주는 n      = {1 + 2 * (z95 / (r - 1))**2:.2f}")
wh = brentq(lambda nu: (r**(1 / 3) - 1 + 2 / (9 * nu)) * math.sqrt(9 * nu / 2) - z95,
            10, 1e6)
print(f"윌슨-힐퍼티가 주는 n   = {wh + 1:.2f}")

r = 1.1025,  r - 1 = 0.1025
    n      정확 p     단순근사 p    WH 근사 p
   12   0.354136    0.405016    0.353926
   50   0.288478    0.305955    0.288326
  200   0.153248    0.153288    0.153206
  538   0.049927    0.046521    0.049925
 1000   0.012817    0.010987    0.012819

p < 0.05 인 n 의 개수 = 19463,  가장 작은 n = 538,  그 구간이 연속인가 = True
p 가 n 에 대해 단조감소인가 = False   (최대 p = 0.358162 at n = 8)

단순근사가 주는 n      = 516.04
윌슨-힐퍼티가 주는 n   = 537.51

import numpy as np
from scipy import stats

def power_var(n, r, alpha=0.05):
    nu = n - 1
    lo = stats.chi2.ppf(alpha / 2, nu)
    hi = stats.chi2.ppf(1 - alpha / 2, nu)
    return stats.chi2.cdf(lo / r, nu) + stats.chi2.sf(hi / r, nu)

rs = [0.5, 0.67, 1.5, 2.0, 3.0]
print(f"{'n':>5s} " + " ".join(f"{'r='+str(r):>9s}" for r in rs))
for n in [10, 20, 30, 50, 100, 200]:
    print(f"{n:5d} " + " ".join(f"{power_var(n, r):9.4f}" for r in rs))

print()
for r in [1.5, 2.0, 3.0, 0.5]:
    n = next(m for m in range(5, 5000) if power_var(m, r) >= 0.80)
    print(f"σ²/σ0² = {r}: 검정력 80%에 필요한 n = {n}")

    n     r=0.5    r=0.67     r=1.5     r=2.0     r=3.0
   10    0.2020    0.0914    0.1833    0.3934    0.7057
   20    0.4650    0.1770    0.2911    0.6289    0.9255
   30    0.6842    0.2687    0.3910    0.7829    0.9830
   50    0.9152    0.4494    0.5623    0.9324    0.9993
  100    0.9987    0.7787    0.8289    0.9974    1.0000
  200    1.0000    0.9788    0.9806    1.0000    1.0000

σ²/σ0² = 1.5: 검정력 80%에 필요한 n = 93
σ²/σ0² = 2.0: 검정력 80%에 필요한 n = 32
σ²/σ0² = 3.0: 검정력 80%에 필요한 n = 13
σ²/σ0² = 0.5: 검정력 80%에 필요한 n = 38

za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
print(f"평균이 0.5σ 이동: n = {np.ceil((za + zb)**2 / 0.5**2):.0f}")
print(f"평균이 1.0σ 이동: n = {np.ceil((za + zb)**2 / 1.0**2):.0f}")

평균이 0.5σ 이동: n = 32
평균이 1.0σ 이동: n = 8

import numpy as np
from scipy import stats

rng = np.random.default_rng(5)
n, M = 30, 20_000
nu = n - 1
lo = stats.chi2.ppf(0.025, nu)
hi = stats.chi2.ppf(0.975, nu)

cases = [("정규", lambda s: rng.normal(0, 1, s), 1.0),
         ("t(5)", lambda s: rng.standard_t(5, s), 5 / 3),
         ("지수", lambda s: rng.exponential(1, s), 1.0)]
print(f"{'분포':>7s} {'카이제곱':>10s} {'첨도 보정':>10s}")
for name, gen, var in cases:
    x = gen((M, n))
    s2 = x.var(1, ddof=1)
    chi = nu * s2 / var
    rej_chi = (chi < lo) | (chi > hi)
    kur = stats.kurtosis(x, axis=1, bias=False)
    z = (np.log(s2) - np.log(var)) / np.sqrt((kur + 2) / nu)
    print(f"{name:>7s} {rej_chi.mean():10.4f} "
          f"{np.mean(np.abs(z) > 1.96):10.4f}")

     분포       카이제곱      첨도 보정
     정규     0.0492     0.0733
   t(5)     0.1766     0.1323
     지수     0.2830     0.1711

import numpy as np
from scipy import stats

n, s2, sigma0sq = 15, 0.0025, 0.0020
nu = n - 1
W = nu * s2 / sigma0sq
print(f"W = {W:.4f}  (자유도 {nu})")
print(f"양측 p = {2 * min(stats.chi2.cdf(W, nu), stats.chi2.sf(W, nu)):.4f}")
print(f"상단 단측 p (σ² > σ0²) = {stats.chi2.sf(W, nu):.4f}")
print(f"하단 단측 p (σ² < σ0²) = {stats.chi2.cdf(W, nu):.4f}")

print(f"\n임계값: 상단 {stats.chi2.ppf(0.95, nu):.3f}, "
      f"하단 {stats.chi2.ppf(0.05, nu):.3f}")
print(f"→ 상단 단측이면 s² 가 "
      f"{stats.chi2.ppf(0.95, nu) * sigma0sq / nu:.6f} 를 넘으면 기각")

W = 17.5000  (자유도 14)
양측 p = 0.4610
상단 단측 p (σ² > σ0²) = 0.2305
하단 단측 p (σ² < σ0²) = 0.7695

임계값: 상단 23.685, 하단 6.571
→ 상단 단측이면 s² 가 0.003384 를 넘으면 기각

import numpy as np
from scipy import stats

rng = np.random.default_rng(88)
x = rng.exponential(2.0, 40)          # 참 분산 4
sigma0sq = 2.5
n, B = len(x), 9_999
s2 = x.var(ddof=1)

# ① 카이제곱 검정
W = (n - 1) * s2 / sigma0sq
p_chi = 2 * min(stats.chi2.cdf(W, n - 1), stats.chi2.sf(W, n - 1))

# ② 부트스트랩: H0 가 참이 되도록 자료를 척도 조정한 뒤 재표본
xc = (x - x.mean()) * np.sqrt(sigma0sq / s2) + x.mean()
idx = rng.integers(0, n, (B, n))
s2b = xc[idx].var(1, ddof=1)
t_obs = np.log(s2 / sigma0sq)
tb = np.log(s2b / sigma0sq)
p_boot = (np.sum(np.abs(tb) >= abs(t_obs)) + 1) / (B + 1)

print(f"s² = {s2:.4f}  (σ0² = {sigma0sq})")
print(f"카이제곱 p = {p_chi:.4f}")
print(f"부트스트랩 p = {p_boot:.4f}")

s² = 3.0517  (σ0² = 2.5)
카이제곱 p = 0.3245
부트스트랩 p = 0.5945

rng = np.random.default_rng(3)
M, n, B = 1_000, 40, 999
c_chi = c_boot = 0
for _ in range(M):
    y = rng.exponential(1.0, n)        # 참 분산 1
    s2 = y.var(ddof=1)
    W = (n - 1) * s2 / 1.0
    c_chi += 2 * min(stats.chi2.cdf(W, n - 1),
                     stats.chi2.sf(W, n - 1)) < 0.05
    yc = (y - y.mean()) / np.sqrt(s2) + y.mean()
    ib = rng.integers(0, n, (B, n))
    s2b = yc[ib].var(1, ddof=1)
    c_boot += (np.sum(np.abs(np.log(s2b)) >= abs(np.log(s2))) + 1) \
        / (B + 1) < 0.05
print(f"지수분포, n=40: 카이제곱 {c_chi / M:.4f}, "
      f"부트스트랩 {c_boot / M:.4f}")

지수분포, n=40: 카이제곱 0.2900, 부트스트랩 0.1310