Bayes 분산 검정¶

import numpy as np
from scipy.stats import invgamma

def posterior_params(x, m0=0.0, k0=1e-6, a0=1e-2, b0=1e-2):
    """정규-역감마 켤레모형의 사후 모수를 구한다.

    기본값은 거의 정보를 주지 않는 사전분포다(k0, a0, b0 가 모두 작다).
    이러면 사후분포가 자료에 거의 전적으로 맡겨진다.
    """
    x = np.asarray(x, dtype=float)
    n = x.size
    xbar = x.mean()
    S = np.sum((x - xbar)**2)
    k_n = k0 + n
    m_n = (k0 * m0 + n * xbar) / k_n
    a_n = a0 + n / 2.0
    b_n = b0 + 0.5 * (S + (k0 * n / k_n) * (xbar - m0)**2)
    return m_n, k_n, a_n, b_n

def draw_posterior_sigma2(x, n_draws=10000, rng=None):
    """분산의 사후분포에서 표본을 뽑는다.

    평균을 적분해 없앤 sigma^2 의 주변 사후분포가 역감마가 된다.
    그래서 MCMC 없이 바로 뽑을 수 있다.
    """
    if rng is None:
        rng = np.random.default_rng()
    m_n, k_n, a_n, b_n = posterior_params(x)
    sig2 = invgamma(a=a_n, scale=b_n).rvs(size=n_draws, random_state=rng)
    return sig2

# === 사전분포의 기여가 얼마나 작은가 ===
from scipy.stats import chi2

x1 = np.array([12, 15, 14, 10, 13, 14, 12, 11], dtype=float)
n = x1.size
S = np.sum((x1 - x1.mean()) ** 2)
m_n, k_n, a_n, b_n = posterior_params(x1)

print(f"n = {n},  xbar = {x1.mean():.4f},  S = {S:.4f},  s^2 = {x1.var(ddof=1):.4f}")
print(f"a_n = {a_n:.4f}   (무정보 극한 n/2 = {n / 2:.4f})")
print(f"b_n = {b_n:.6f}   (무정보 극한 S/2 = {S / 2:.6f})")
print(f"  b0 가 더한 것        = {1e-2:.6f}")
print(f"  평균 항이 더한 것    = {0.5 * (1e-6 * n / k_n) * x1.mean() ** 2:.3e}")
print(f"  둘을 합쳐 b_n 의 {100 * (b_n - S / 2) / b_n:.4f} %")

# === 추축량: 2 b_n / sigma^2 ~ chi2(2 a_n) 인가 ===
lo, hi = invgamma(a=a_n, scale=b_n).ppf([0.025, 0.975])
piv = 2 * b_n / chi2(2 * a_n).ppf([0.975, 0.025])
print(f"\n역감마 분위수로  : ({lo:.6f}, {hi:.6f})")
print(f"2b_n/chi2(2a_n) 로: ({piv[0]:.6f}, {piv[1]:.6f})")

# === 신용구간 대 빈도주의 신뢰구간 ===
f_lo, f_hi = S / chi2(n - 1).ppf([0.975, 0.025])
print(f"\n95% 신용구간      = ({lo:.3f}, {hi:.3f})      상한/하한 = {hi / lo:.4f}")
print(f"95% 신뢰구간      = ({f_lo:.3f}, {f_hi:.3f})     상한/하한 = {f_hi / f_lo:.4f}")
print(f"로그 폭           = {np.log(hi / lo):.4f} 대 {np.log(f_hi / f_lo):.4f}"
      f"   -> {100 * (1 - np.log(hi / lo) / np.log(f_hi / f_lo)):.1f} % 좁다")
print(f"자유도 비교: 2a_n = {2 * a_n:.2f}  대  n-1 = {n - 1}")
print(f"chi2 분위수 비 (자료와 무관): {chi2(2 * a_n).ppf(0.975) / chi2(2 * a_n).ppf(0.025):.4f}"
      f"  대  {chi2(n - 1).ppf(0.975) / chi2(n - 1).ppf(0.025):.4f}")

# === a0 = -1/2 로 두면 빈도주의와 정확히 같아진다 ===
g_lo, g_hi = invgamma(a=(n - 1) / 2, scale=S / 2).ppf([0.025, 0.975])
print(f"\na0 = -1/2 (a_n = (n-1)/2) 신용구간 = ({g_lo:.6f}, {g_hi:.6f})")
print(f"빈도주의 신뢰구간                  = ({f_lo:.6f}, {f_hi:.6f})")
print(f"차 = {abs(g_lo - f_lo):.3e}, {abs(g_hi - f_hi):.3e}")

# === 함수가 제대로 도는지 ===
rng = np.random.default_rng(0)
d = draw_posterior_sigma2(x1, n_draws=200_000, rng=rng)
print(f"\n사후표본 20 만 개:  평균 {d.mean():.4f} (닫힌 꼴 {b_n / (a_n - 1):.4f}),  "
      f"중앙값 {np.median(d):.4f} (닫힌 꼴 {invgamma(a=a_n, scale=b_n).ppf(0.5):.4f})")

n = 8,  xbar = 12.6250,  S = 19.8750,  s^2 = 2.8393
a_n = 4.0100   (무정보 극한 n/2 = 4.0000)
b_n = 9.947580   (무정보 극한 S/2 = 9.937500)
  b0 가 더한 것        = 0.010000
  평균 항이 더한 것    = 7.970e-05
  둘을 합쳐 b_n 의 0.1013 %

역감마 분위수로  : (1.132684, 9.085126)
2b_n/chi2(2a_n) 로: (1.132684, 9.085126)

95% 신용구간      = (1.133, 9.085)      상한/하한 = 8.0209
95% 신뢰구간      = (1.241, 11.761)     상한/하한 = 9.4757
로그 폭           = 2.0820 대 2.2487   -> 7.4 % 좁다
자유도 비교: 2a_n = 8.02  대  n-1 = 7
chi2 분위수 비 (자료와 무관): 8.0209  대  9.4757

a0 = -1/2 (a_n = (n-1)/2) 신용구간 = (1.241197, 11.761265)
빈도주의 신뢰구간                  = (1.241197, 11.761265)
차 = 0.000e+00, 0.000e+00

사후표본 20 만 개:  평균 3.2993 (닫힌 꼴 3.3048),  중앙값 2.6953 (닫힌 꼴 2.7016)

x1 = np.array([12, 15, 14, 10, 13, 14, 12, 11], dtype=float)
x2 = np.array([22, 25, 20, 18, 24, 23, 19, 21], dtype=float)

# 두 집단의 사후표본을 각각 뽑아 나눈다. 이 비의 분포가 곧 분산비의
# 사후분포다. 빈도주의 F 검정이 p-값 하나를 주는 자리에서, 베이즈 쪽은
# "비가 1 보다 클 확률"을 그대로 셈할 수 있다.
rng = np.random.default_rng(0)
s1 = draw_posterior_sigma2(x1, n_draws=20000, rng=rng)
s2 = draw_posterior_sigma2(x2, n_draws=20000, rng=rng)
ratio = s1 / s2

print(f"Sample variances:           {np.var(x1, ddof=1):.4f}, "
      f"{np.var(x2, ddof=1):.4f}")
print(f"Posterior mean of sigma1^2: {s1.mean():.3f}")
print(f"Posterior mean of sigma2^2: {s2.mean():.3f}")
print(f"Posterior mean of ratio:    {ratio.mean():.3f}")
print(f"95% credible interval:      ({np.percentile(ratio, 2.5):.3f}, "
      f"{np.percentile(ratio, 97.5):.3f})")
print(f"P(sigma1^2 > sigma2^2):     {np.mean(ratio > 1.0):.4f}")

# === 닫힌 꼴과 맞춰 본다 ===
from scipy.stats import f as f_dist

_, _, a1, b1 = posterior_params(x1)
_, _, a2, b2 = posterior_params(x2)
assert a1 == a2, "두 집단의 a_n 이 같아야 F 로 적힌다"
Fd = f_dist(2 * a1, 2 * a2)
sc = b1 / b2

print(f"\na_n = {a1:.4f} (두 집단 공통),  b_n = {b1:.6f}, {b2:.6f}")
print(f"rho ~ (b1/b2) x F(2a,2a) = {sc:.6f} x F({2 * a1:.2f},{2 * a2:.2f})")
print("                닫힌 꼴    모의 2만개")
for name, closed, mc in (
        ("median    ", sc * Fd.ppf(0.5), np.median(ratio)),
        ("mean      ", (b1 / (a1 - 1)) * (a2 / b2), ratio.mean()),
        ("q(0.025)  ", sc * Fd.ppf(0.025), np.percentile(ratio, 2.5)),
        ("q(0.975)  ", sc * Fd.ppf(0.975), np.percentile(ratio, 97.5)),
        ("P(rho>1)  ", Fd.sf(b2 / b1), np.mean(ratio > 1))):
    print(f"  {name}  {closed:10.6f}  {mc:12.6f}")

print(f"\n표본분산의 비 s1^2/s2^2 = {np.var(x1, ddof=1) / np.var(x2, ddof=1):.6f}")
print(f"사후중앙값 b1/b2        = {sc:.6f}   (차 {abs(sc - np.var(x1, ddof=1) / np.var(x2, ddof=1)):.2e})")
print(f"평균/중앙값 = {((b1 / (a1 - 1)) * (a2 / b2)) / sc:.4f} 배")

# === 빈도주의 F 검정과 견준다 ===
n = x1.size
Fs = np.var(x1, ddof=1) / np.var(x2, ddof=1)
p_two = 2 * min(f_dist(n - 1, n - 1).cdf(Fs), f_dist(n - 1, n - 1).sf(Fs))
fl, fh = Fs / f_dist(n - 1, n - 1).ppf([0.975, 0.025])
print(f"\nF 검정: F = {Fs:.4f}, 양측 p = {p_two:.4f}, 단측(하단) p = {p_two / 2:.4f}")
print(f"  P(rho > 1) = {Fd.sf(b2 / b1):.4f}  <- 단측 p 와 견주라")
print(f"빈도주의 신뢰구간 = ({fl:.6f}, {fh:.6f})   로그 폭 {np.log(fh / fl):.4f}")
print(f"베이즈 신용구간   = ({sc * Fd.ppf(0.025):.6f}, {sc * Fd.ppf(0.975):.6f})"
      f"   로그 폭 {np.log(Fd.ppf(0.975) / Fd.ppf(0.025)):.4f}")
print(f"  -> {100 * (1 - np.log(Fd.ppf(0.975) / Fd.ppf(0.025)) / np.log(fh / fl)):.1f} % 좁다")

Sample variances:           2.8393, 6.0000
Posterior mean of sigma1^2: 3.317
Posterior mean of sigma2^2: 6.956
Posterior mean of ratio:    0.634
95% credible interval:      (0.106, 2.074)
P(sigma1^2 > sigma2^2):     0.1573

a_n = 4.0100 (두 집단 공통),  b_n = 9.947580, 21.010231
rho ~ (b1/b2) x F(2a,2a) = 0.473464 x F(8.02,8.02)
                닫힌 꼴    모의 2만개
  median        0.473464      0.477285
  mean          0.630760      0.633695
  q(0.025)      0.107025      0.106483
  q(0.975)      2.094531      2.073990
  P(rho>1)      0.155034      0.157350

표본분산의 비 s1^2/s2^2 = 0.473214
사후중앙값 b1/b2        = 0.473464   (차 2.49e-04)
평균/중앙값 = 1.3322 배

F 검정: F = 0.4732, 양측 p = 0.3448, 단측(하단) p = 0.1724
  P(rho > 1) = 0.1550  <- 단측 p 와 견주라
빈도주의 신뢰구간 = (0.094739, 2.363662)   로그 폭 3.2168
베이즈 신용구간   = (0.107025, 2.094531)   로그 폭 2.9740
  -> 7.5 % 좁다

import numpy as np
from scipy.stats import invgamma

def posterior_params(x, m0=0.0, k0=1e-6, a0=1e-2, b0=1e-2):
    x = np.asarray(x, dtype=float)
    n = x.size
    xbar = x.mean()
    S = np.sum((x - xbar)**2)
    k_n = k0 + n
    a_n = a0 + n / 2.0
    b_n = b0 + 0.5 * (S + (k0 * n / k_n) * (xbar - m0)**2)
    return a_n, b_n

x1 = np.array([12, 15, 14, 10, 13, 14, 12, 11], dtype=float)
rng = np.random.default_rng(0)

a_v, b_v = posterior_params(x1, a0=1e-2, b0=1e-2)      # vague
draws_v = invgamma(a=a_v, scale=b_v).rvs(20000, random_state=rng)

a_i, b_i = posterior_params(x1, a0=3, b0=10)           # informative
draws_i = invgamma(a=a_i, scale=b_i).rvs(20000, random_state=rng)

print(f"Vague:       a={a_v:.2f}, b={b_v:.2f}, mean={draws_v.mean():.2f}, "
      f"CI=({np.percentile(draws_v,2.5):.2f}, {np.percentile(draws_v,97.5):.2f})")
print(f"Informative: a={a_i:.2f}, b={b_i:.2f}, mean={draws_i.mean():.2f}, "
      f"CI=({np.percentile(draws_i,2.5):.2f}, {np.percentile(draws_i,97.5):.2f})")

Vague:       a=4.01, b=9.95, mean=3.32, CI=(1.13, 9.11)
Informative: a=7.00, b=19.94, mean=3.32, CI=(1.53, 7.09)

import numpy as np
from scipy.stats import invgamma, f as f_dist

def posterior_params(x, m0=0.0, k0=1e-6, a0=1e-2, b0=1e-2):
    n = x.size
    xbar = x.mean()
    S = np.sum((x - xbar)**2)
    k_n = k0 + n
    a_n = a0 + n / 2.0
    b_n = b0 + 0.5 * (S + (k0 * n / k_n) * (xbar - m0)**2)
    return a_n, b_n

rng = np.random.default_rng(42)
x1 = rng.normal(0, 2, 30)     # sd = 2, so variance = 4
x2 = rng.normal(0, 3, 30)     # sd = 3, so variance = 9

print(f"sample variances: {np.var(x1, ddof=1):.3f}, "
      f"{np.var(x2, ddof=1):.3f} (true: 4, 9)")

a1, b1 = posterior_params(x1)
a2, b2 = posterior_params(x2)
s1 = invgamma(a=a1, scale=b1).rvs(50000, random_state=rng)
s2 = invgamma(a=a2, scale=b2).rvs(50000, random_state=rng)
ratio = s1 / s2

print(f"P(rho < 1 | data) = {np.mean(ratio < 1):.4f}")

F_stat = np.var(x1, ddof=1) / np.var(x2, ddof=1)
p_val = 2 * min(f_dist.cdf(F_stat, 29, 29), f_dist.sf(F_stat, 29, 29))
print(f"F-test: F = {F_stat:.4f}, p = {p_val:.4f}")

sample variances: 2.413, 5.828 (true: 4, 9)
P(rho < 1 | data) = 0.9908
F-test: F = 0.4141, p = 0.0205

import numpy as np
from scipy.stats import invgamma, skew
import matplotlib.pyplot as plt

def posterior_params(x, m0=0.0, k0=1e-6, a0=1e-2, b0=1e-2):
    n = x.size; xbar = x.mean(); S = np.sum((x - xbar)**2)
    k_n = k0 + n
    a_n = a0 + n / 2.0
    b_n = b0 + 0.5 * (S + (k0 * n / k_n) * (xbar - m0)**2)
    return a_n, b_n

x1 = np.array([12, 15, 14, 10, 13, 14, 12, 11], dtype=float)
x2 = np.array([22, 25, 20, 18, 24, 23, 19, 21], dtype=float)
rng = np.random.default_rng(0)

a1, b1 = posterior_params(x1)
a2, b2 = posterior_params(x2)
s1 = invgamma(a=a1, scale=b1).rvs(50000, random_state=rng)
s2 = invgamma(a=a2, scale=b2).rvs(50000, random_state=rng)
ratio = s1 / s2

print(f"skewness of ratio:     {skew(ratio):.3f}")
print(f"skewness of log ratio: {skew(np.log(ratio)):.3f}")

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].hist(ratio, bins=60, edgecolor='k', alpha=0.7)
axes[0].set_title("Posterior of ratio")
axes[1].hist(np.log(ratio), bins=60, edgecolor='k', alpha=0.7)
axes[1].set_title("Posterior of log(ratio)")
plt.tight_layout()
plt.show()

skewness of ratio:     9.518
skewness of log ratio: -0.017