분산에 대한 카이제곱 검정¶

import numpy as np
import scipy.stats as stats

def chi2_test_for_variance(data, sigma2_0=1.0):
    """모분산이 sigma2_0 인지 검정한다.

    카이제곱 분포는 좌우가 대칭이 아니다. 그래서 양측 p-값을 만들 때
    한쪽 꼬리 넓이의 두 배를 쓰되, 둘 중 작은 쪽을 골라 두 배 한다.
    """
    n = len(data)
    s2 = np.var(data, ddof=1)
    statistic = (n - 1) * s2 / sigma2_0
    p_value = 2 * min(
        stats.chi2(df=n - 1).cdf(statistic),
        stats.chi2(df=n - 1).sf(statistic)
    )
    return statistic, p_value

# (1) 두 기각규칙이 같은 결정을 주는가. 1000 개 표본에서 맞춰 본다.
rng = np.random.default_rng(0)
n, alpha = 100, 0.05
lo = stats.chi2(df=n - 1).ppf(alpha / 2)
hi = stats.chi2(df=n - 1).ppf(1 - alpha / 2)
same, pmax = True, 0.0
for _ in range(1000):
    x = rng.normal(size=n)
    T, p = chi2_test_for_variance(x)
    pmax = max(pmax, p)
    if (p < alpha) != (T < lo or T > hi):
        same = False
print(f"임계값 규칙 = p-값 규칙 ?  {same}")
print(f"관측된 p-값의 최대      = {pmax:.4f}  (1 을 넘지 않는다)")

# (2) 실제 1종오류율. H0: sigma^2 = 1 이 참인 두 모집단에서 센다.
B = 20_000
print(f"\n{'모집단':>10}{'beta2':>8}{'n=20':>9}{'n=100':>9}{'n=400':>9}{'극한(이론)':>12}")
for name, draw, beta2 in [
    ("정규", lambda m: rng.normal(size=m), 3.0),
    ("지수", lambda m: rng.exponential(size=m), 9.0),
]:
    row = []
    for n in (20, 100, 400):
        T = (n - 1) * np.var(draw((B, n)), axis=1, ddof=1)
        c = stats.chi2(df=n - 1)
        row.append(np.mean((T < c.ppf(0.025)) | (T > c.ppf(0.975))))
    lim = 2 * stats.norm.sf(1.96 * np.sqrt(2 / (beta2 - 1)))
    print(f"{name:>10}{beta2:>8.1f}" + "".join(f"{v:>9.3f}" for v in row) + f"{lim:>12.3f}")

임계값 규칙 = p-값 규칙 ?  True
관측된 p-값의 최대      = 0.9990  (1 을 넘지 않는다)

       모집단   beta2     n=20    n=100    n=400      극한(이론)
        정규     3.0    0.049    0.049    0.051       0.050
        지수     9.0    0.269    0.311    0.324       0.327

seed = 1
x = stats.norm(loc=0, scale=1).rvs(100, random_state=seed)

for scale in [1.00, 1.05, 1.10, 1.15, 1.20]:
    # random_state를 고정했으므로 scale만 바뀐 **같은** 난수열이 쓰인다.
    # 표본마다 새로 뽑았다면 표집 변동이 섞여 scale의 효과만 보기 어렵다.
    y = stats.norm(loc=1, scale=scale).rvs(100, random_state=seed)
    stat, pval = chi2_test_for_variance(y, sigma2_0=1.0)
    print(f"sigma={scale:.2f}: s2={np.var(y, ddof=1):.3f} T={stat:.2f} p={pval:.3f}")

sigma=1.00: s2=0.791 T=78.35 p=0.125
sigma=1.05: s2=0.873 T=86.38 p=0.373
sigma=1.10: s2=0.958 T=94.80 p=0.799
sigma=1.15: s2=1.047 T=103.62 p=0.711
sigma=1.20: s2=1.140 T=112.82 p=0.324

import numpy as np
from scipy import stats

def chi2_var_p(x, s0=1.0):
    n = len(x)
    T = (n - 1) * x.var(ddof=1) / s0
    return 2 * min(stats.chi2(n - 1).cdf(T), stats.chi2(n - 1).sf(T))

rng = np.random.default_rng(111)
B = 20_000
# 모두 분산이 정확히 1 이 되도록 척도를 맞춘 분포들
dists = {
    "정규 (초과첨도 0)": lambda n: rng.normal(0, 1, n),
    "균등 (−1.2)": lambda n: rng.uniform(-np.sqrt(3), np.sqrt(3), n),
    "라플라스 (+3)": lambda n: rng.laplace(0, 1 / np.sqrt(2), n),
    "t(8) (+1.5)": lambda n: rng.standard_t(8, n) / np.sqrt(8 / 6),
    "t(5) (+6)": lambda n: rng.standard_t(5, n) / np.sqrt(5 / 3),
    "지수 (+6)": lambda n: rng.exponential(1, n),
}
for n in [20, 50]:
    print(f"  n = {n}")
    for lab, f in dists.items():
        a = sum(chi2_var_p(f(n)) < 0.05 for _ in range(B))
        print(f"    {lab:>18s} {a / B:8.4f}")

  n = 20
           정규 (초과첨도 0)   0.0485
             균등 (−1.2)   0.0046
             라플라스 (+3)   0.1808
           t(8) (+1.5)   0.1033
             t(5) (+6)   0.1583
               지수 (+6)   0.2687
  n = 50
           정규 (초과첨도 0)   0.0496
             균등 (−1.2)   0.0024
             라플라스 (+3)   0.2006
           t(8) (+1.5)   0.1201
             t(5) (+6)   0.1880
               지수 (+6)   0.2969

import numpy as np
from scipy import stats

rng = np.random.default_rng(222)
B = 10_000

def chi2_ci(x, a=0.05):
    n = len(x)
    s2 = x.var(ddof=1)
    return ((n - 1) * s2 / stats.chi2.ppf(1 - a / 2, n - 1),
            (n - 1) * s2 / stats.chi2.ppf(a / 2, n - 1))

def boot_ci(x, a=0.05, R=800, rng=None):
    n = len(x)
    bs = np.array([x[rng.integers(0, n, n)].var(ddof=1) for _ in range(R)])
    return np.quantile(bs, [a / 2, 1 - a / 2])

dists = {
    "정규": lambda n: rng.normal(0, 1, n),
    "라플라스": lambda n: rng.laplace(0, 1 / np.sqrt(2), n),
    "t(5)": lambda n: rng.standard_t(5, n) / np.sqrt(5 / 3),
    "지수": lambda n: rng.exponential(1, n),
    "균등": lambda n: rng.uniform(-np.sqrt(3), np.sqrt(3), n),
}
print("참값 σ² = 1 을 덮는 비율 (명목 95%)")
print(f"{'분포':>10s} {'n=20 카이제곱':>13s} {'n=50 카이제곱':>13s} "
      f"{'n=50 부트스트랩':>15s}")
for lab, f in dists.items():
    row = []
    for n in [20, 50]:
        c = 0
        for _ in range(B):
            lo, hi = chi2_ci(f(n))
            c += lo <= 1.0 <= hi
        row.append(c / B)
    c = 0
    for _ in range(1_500):
        lo, hi = boot_ci(f(50), rng=rng)
        c += lo <= 1.0 <= hi
    row.append(c / 1_500)
    print(f"{lab:>10s} {row[0]:13.4f} {row[1]:13.4f} {row[2]:15.4f}")

참값 σ² = 1 을 덮는 비율 (명목 95%)
        분포     n=20 카이제곱     n=50 카이제곱      n=50 부트스트랩
        정규        0.9470        0.9477          0.9107
      라플라스        0.8140        0.8032          0.8660
      t(5)        0.8477        0.8090          0.8293
        지수        0.7413        0.6944          0.8300
        균등        0.9952        0.9964          0.9280

from scipy import stats

def n_needed(ratio, power=0.80, alpha=0.05):
    for n in range(3, 200_000):
        lo = stats.chi2.ppf(alpha / 2, n - 1)
        hi = stats.chi2.ppf(1 - alpha / 2, n - 1)
        pw = (stats.chi2.cdf(lo / ratio, n - 1)
              + stats.chi2.sf(hi / ratio, n - 1))
        if pw >= power:
            return n

print(f"{'σ²/σ0²':>8s} {'필요한 n':>9s}")
for r in [1.25, 1.5, 2.0, 3.0, 0.8, 0.667, 0.5]:
    print(f"{r:8.3f} {n_needed(r):9d}")

  σ²/σ0²     필요한 n
   1.250       308
   1.500        93
   2.000        32
   3.000        13
   0.800       327
   0.667       103
   0.500        38

import numpy as np
from scipy import stats, optimize

def shortest(n, conf=0.95):
    """길이 1/a - 1/b 를 최소화하는 (a, b) 를 찾는다."""
    df = n - 1

    def length(alo):
        a = stats.chi2.ppf(alo, df)
        b = stats.chi2.ppf(alo + conf, df)
        return 1 / a - 1 / b

    r = optimize.minimize_scalar(length, bounds=(1e-6, 1 - conf - 1e-6),
                                 method="bounded")
    return (stats.chi2.ppf(r.x, df), stats.chi2.ppf(r.x + conf, df), r.fun)

print("동일꼬리 구간 대 최단 구간 (95%, (n-1)S² = 1 로 둔 길이)")
print(f"{'n':>5s} {'동일꼬리 [a,b]':>22s} {'길이':>8s} "
      f"{'최단 [a,b]':>22s} {'길이':>8s} {'절감':>7s}")
for n in [5, 10, 20, 50, 100]:
    df = n - 1
    a1, b1 = stats.chi2.ppf(0.025, df), stats.chi2.ppf(0.975, df)
    L1 = 1 / a1 - 1 / b1
    a2, b2, L2 = shortest(n)
    print(f"{n:5d} [{a1:8.3f},{b1:9.3f}] {L1:8.4f} "
          f"[{a2:8.3f},{b2:9.3f}] {L2:8.4f} {100 * (L1 - L2) / L1:6.2f}%")

print("\n최단 구간의 꼬리 확률 배분")
for n in [5, 20, 100]:
    a2, b2, _ = shortest(n)
    print(f"  n={n:3d}: 아래 꼬리 {stats.chi2.cdf(a2, n - 1):.4f}, "
          f"위 꼬리 {stats.chi2.sf(b2, n - 1):.4f}")

동일꼬리 구간 대 최단 구간 (95%, (n-1)S² = 1 로 둔 길이)
    n             동일꼬리 [a,b]       길이               최단 [a,b]       길이      절감
    5 [   0.484,   11.143]   1.9746 [   0.708,   21.052]   1.3645  30.90%
   10 [   2.700,   19.023]   0.3177 [   3.284,   26.076]   0.2662  16.22%
   20 [   8.907,   32.852]   0.0818 [   9.899,   38.327]   0.0749   8.44%
   50 [  31.555,   70.222]   0.0175 [  33.044,   74.517]   0.0168   3.48%
  100 [  73.361,  128.422]   0.0058 [  75.151,  132.189]   0.0057   1.76%

최단 구간의 꼬리 확률 배분
  n=  5: 아래 꼬리 0.0497, 위 꼬리 0.0003
  n= 20: 아래 꼬리 0.0446, 위 꼬리 0.0054
  n=100: 아래 꼬리 0.0355, 위 꼬리 0.0145

분산에 대한 추론이 필요하다
    │
    ├─ 정규성을 신뢰할 수 있는가?
    │     │
    │     ├─ 예 (물리 측정, 공정 관리 등)
    │     │     └─→ 카이제곱 검정 / 구간
    │     │
    │     └─ 아니오
    │           ├─→ BCa 부트스트랩 구간
    │           ├─→ 변환 후 검정 (로그 등)
    │           └─→ 로버스트 산포 측도 (IQR, MAD)
    │
    └─ 표본이 충분한가?  (비 1.25 → n ≈ 300)

n = 50, S² = 0.0612 mm²
H0: σ² = 0.04 mm²
T = 49 × 0.0612 / 0.04 = 74.97,  df = 49
p = 2 × min(0.9910, 0.0090) = 0.0180
95% 신뢰구간: σ² ∈ [0.0427, 0.0951]
→ 목표 분산을 넘어섰다고 판단 (단, 정규성 가정 아래)