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

import numpy as np
import scipy.stats as stats


def chi2_test_for_variance(data, sigma2_0=1.0):
    """모분산이 sigma2_0 인지 검정하는 일표본 카이제곱 검정.

    H0: sigma^2 = sigma2_0
    H1: sigma^2 != 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


# 한 쌍으로 손풀기. H0 이 참인 표본이므로 p 가 커야 한다.
rng = np.random.default_rng(5)
stat, pval = chi2_test_for_variance(rng.normal(size=25), sigma2_0=1.0)
print(f"표본 하나:  T = {stat:.4f},  p = {pval:.4f}")

n, alpha = 25, 0.05
dist = stats.chi2(df=n - 1)
lo, hi = dist.ppf([alpha / 2, 1 - alpha / 2])
print(f"\nchi2(24) 의 2.5%/97.5% 임계값 = {lo:.4f}  {hi:.4f}")
print(f"  경계에서의 p:  p({lo:.4f}) = {2 * min(dist.cdf(lo), dist.sf(lo)):.6f}"
      f"   p({hi:.4f}) = {2 * min(dist.cdf(hi), dist.sf(hi)):.6f}")

# (1) 두 판정규칙이 정말 같은 기각역을 주는가. 임계값 바로 안팎까지 격자에 끼운다.
grid = np.unique(np.concatenate([
    np.linspace(1e-6, 120, 200_001),
    lo + np.array([-1e-9, 0.0, 1e-9]),
    hi + np.array([-1e-9, 0.0, 1e-9]),
]))
U = dist.cdf(grid)
p_grid = 2 * np.minimum(U, 1 - U)
rule_crit = (grid < lo) | (grid > hi)      # 임계값으로 자르는 규칙
rule_pval = p_grid < alpha                 # p < alpha 규칙
print(f"\n격자점 {grid.size:,} 개에서 두 규칙이 어긋난 횟수 = {(rule_crit != rule_pval).sum()}")
print(f"격자에서 p 의 최대값 = {p_grid.max():.6f}  (1 을 넘지 않는다)")

# 비대칭이라는 것이 눈에 보이게. 평균은 n-1 = 24 지만 p = 1 이 되는 곳은 중앙값이다.
print(f"\nchi2(24) 의 평균 = {n - 1},  중앙값 = {dist.median():.4f}")
print(f"  기각역 하단은 평균에서 {(n - 1) - lo:.4f} 아래,  상단은 {hi - (n - 1):.4f} 위")
print(f"  p(평균 24) = {2 * min(dist.cdf(24), dist.sf(24)):.6f}")
print(f"  p(중앙값)  = {2 * min(dist.cdf(dist.median()), dist.sf(dist.median())):.6f}")

표본 하나:  T = 18.8045,  p = 0.4753

chi2(24) 의 2.5%/97.5% 임계값 = 12.4012  39.3641
  경계에서의 p:  p(12.4012) = 0.050000   p(39.3641) = 0.050000

격자점 200,007 개에서 두 규칙이 어긋난 횟수 = 2
격자에서 p 의 최대값 = 0.999968  (1 을 넘지 않는다)

chi2(24) 의 평균 = 24,  중앙값 = 23.3367
  기각역 하단은 평균에서 11.5988 아래,  상단은 15.3641 위
  p(평균 24) = 0.923195
  p(중앙값)  = 1.000000

# 참 표준편차를 1 에서 1.2 까지 올려 가며 검정이 언제부터 잡아내는지 본다.
size, seed = 100, 0

null = stats.chi2(df=size - 1)
lo, hi = null.ppf([0.025, 0.975])
T1, _ = chi2_test_for_variance(
    stats.norm(loc=1, scale=1.0).rvs(size, random_state=seed), sigma2_0=1.0)
print(f"chi2(99) 의 2.5%/97.5% 임계값 = {lo:.4f}  {hi:.4f}")
print(f"sigma = 1 일 때의 T(1) = {T1:.4f}")
print(f"기각 문턱:  sigma > sqrt({hi:.4f}/{T1:.4f}) = {np.sqrt(hi / T1):.4f}\n")

print(f"{'sigma':>6}{'s':>9}{'T':>10}{'sigma^2*T(1)':>14}{'차':>10}{'p':>8}  판정")
for scale in [1.00, 1.05, 1.10, 1.15, 1.20]:
    y = stats.norm(loc=1, scale=scale).rvs(size, random_state=seed)
    stat, pval = chi2_test_for_variance(y, sigma2_0=1.0)
    pred = scale**2 * T1
    verdict = "기각" if (stat < lo or stat > hi) else "기각 못함"
    print(f"{scale:>6.2f}{y.std(ddof=1):>9.4f}{stat:>10.4f}{pred:>14.4f}"
          f"{abs(stat - pred):>10.1e}{pval:>8.3f}  {verdict}")

# 다섯 표본이 같은 z 를 척도만 바꾸어 쓴 것인지 직접 확인한다.
z = stats.norm(loc=0, scale=1).rvs(size, random_state=seed)
y12 = stats.norm(loc=1, scale=1.2).rvs(size, random_state=seed)
print(f"\ny(1.2) - (1 + 1.2*z) 의 최대 절대오차 = {np.abs(y12 - (1 + 1.2 * z)).max():.3e}")

# 독립 표본이면 기각률이 얼마인가. T_obs = sigma^2 * W,  W ~ chi2(99) 이므로
print("\n독립 표본일 때 정확한 기각률 (n=100, alpha=0.05)")
print(f"{'sigma':>6}{'일표본 chi2':>14}{'두표본 F':>12}")
f_null = stats.f(size - 1, size - 1)
f_lo, f_hi = f_null.ppf([0.025, 0.975])
for scale in [1.00, 1.05, 1.10, 1.15, 1.20, 1.30]:
    r = scale**2
    chi_pow = null.cdf(lo / r) + null.sf(hi / r)
    f_pow = f_null.cdf(f_lo * r) + f_null.sf(f_hi * r)
    print(f"{scale:>6.2f}{chi_pow:>14.4f}{f_pow:>12.4f}")

chi2(99) 의 2.5%/97.5% 임계값 = 73.3611  128.4220
sigma = 1 일 때의 T(1) = 101.5827
기각 문턱:  sigma > sqrt(128.4220/101.5827) = 1.1244

 sigma        s         T  sigma^2*T(1)         차       p  판정
  1.00   1.0130  101.5827      101.5827   0.0e+00   0.819  기각 못함
  1.05   1.0636  111.9949      111.9949   2.8e-14   0.351  기각 못함
  1.10   1.1143  122.9150      122.9150   4.3e-14   0.104  기각 못함
  1.15   1.1649  134.3431      134.3431   5.7e-14   0.021  기각
  1.20   1.2156  146.2790      146.2790   5.7e-14   0.003  기각

y(1.2) - (1 + 1.2*z) 의 최대 절대오차 = 0.000e+00

독립 표본일 때 정확한 기각률 (n=100, alpha=0.05)
 sigma      일표본 chi2       두표본 F
  1.00        0.0500      0.0500
  1.05        0.1158      0.0769
  1.10        0.2946      0.1559
  1.15        0.5352      0.2817
  1.20        0.7500      0.4378
  1.30        0.9586      0.7381

import scipy.stats as stats

T = 24 * 0.06 / 0.04  # 36.0
p = 2 * min(stats.chi2(df=24).cdf(T), stats.chi2(df=24).sf(T))
print(f"T = {T:.1f}, right tail = {stats.chi2(df=24).sf(T):.4f}, "
      f"two-sided p = {p:.4f}")

T = 36.0, right tail = 0.0549, two-sided p = 0.1098

import numpy as np
import scipy.stats as stats

rng = np.random.default_rng(42)
n, sigma2_0, alpha, n_sims = 30, 1.0, 0.05, 10000
rejections = 0

for _ in range(n_sims):
    x = rng.normal(0, 1, size=n)
    s2 = np.var(x, ddof=1)
    T = (n - 1) * s2 / sigma2_0
    p = 2 * min(stats.chi2(df=n - 1).cdf(T), stats.chi2(df=n - 1).sf(T))
    if p < alpha:
        rejections += 1

print(f"Empirical Type I error: {rejections / n_sims:.4f}")

Empirical Type I error: 0.0500

import numpy as np
import scipy.stats as stats

rng = np.random.default_rng(42)
n, sigma2_0, alpha, n_sims = 30, 5 / 3, 0.05, 10000
# Var(t(5)) = 5/(5-2) = 5/3
rejections = 0

for _ in range(n_sims):
    x = stats.t(df=5).rvs(size=n, random_state=rng)
    s2 = np.var(x, ddof=1)
    T = (n - 1) * s2 / sigma2_0
    p = 2 * min(stats.chi2(df=n - 1).cdf(T), stats.chi2(df=n - 1).sf(T))
    if p < alpha:
        rejections += 1

print(f"Empirical rejection rate (t(5)): {rejections / n_sims:.4f}")

Empirical rejection rate (t(5)): 0.1725

import numpy as np
import scipy.stats as stats

def chi2_test_variance(data, sigma2_0=1.0, alternative='two-sided'):
    """모분산에 대한 일표본 카이제곱 검정.

    alternative: 'two-sided', 'greater', or 'less'
    """
    n = len(data)
    s2 = np.var(data, ddof=1)
    T = (n - 1) * s2 / sigma2_0
    dist = stats.chi2(df=n - 1)
    if alternative == 'greater':
        p = dist.sf(T)
    elif alternative == 'less':
        p = dist.cdf(T)
    elif alternative == 'two-sided':
        p = 2 * min(dist.cdf(T), dist.sf(T))
    else:
        raise ValueError(f"unknown alternative: {alternative}")
    return T, p