두 분산 비교를 위한 F 검정¶

from scipy.stats import f

df1, df2, alpha = 14, 19, 0.05
print(f"lower = {f.ppf(alpha / 2, df1, df2):.4f}")
print(f"upper = {f.ppf(1 - alpha / 2, df1, df2):.4f}")

# 관측값과 양측 p 값
lo, hi = f.ppf(alpha / 2, df1, df2), f.ppf(1 - alpha / 2, df1, df2)
F = 25 / 16
p = 2 * min(f.cdf(F, df1, df2), f.sf(F, df1, df2))
print(f"F = 25/16 = {F:.4f}   ->  {'기각' if (F < lo or F > hi) else '기각 못함'},  p = {p:.4f}")

# 기각하려면 분산비가 얼마나 되어야 하는가. 두 방향이 다르다.
print(f"\n집단 1 이 더 퍼져 있어야 할 때:  s1^2/s2^2 > {hi:.4f}")
print(f"집단 2 가 더 퍼져 있어야 할 때:  s2^2/s1^2 > {1 / lo:.4f}")

# 두 집단이 같은 크기 n 일 때의 문턱
print("\n두 집단이 모두 크기 n 일 때 기각에 필요한 분산비")
for n in (5, 15, 30, 100):
    print(f"  n = {n:>3}  ->  {f.ppf(0.975, n - 1, n - 1):.4f} 배")

lower = 0.3496
upper = 2.6469
F = 25/16 = 1.5625   ->  기각 못함,  p = 0.3607

집단 1 이 더 퍼져 있어야 할 때:  s1^2/s2^2 > 2.6469
집단 2 가 더 퍼져 있어야 할 때:  s2^2/s1^2 > 2.8607

두 집단이 모두 크기 n 일 때 기각에 필요한 분산비
  n =   5  ->  9.6045 배
  n =  15  ->  2.9786 배
  n =  30  ->  2.1010 배
  n = 100  ->  1.4862 배

import numpy as np
from scipy.stats import f

sample1 = [12, 15, 14, 10, 13, 14, 12, 11]
sample2 = [22, 25, 20, 18, 24, 23, 19, 21]

# ddof=1 로 표본분산을 쓴다. F 검정의 이론이 그렇게 세워져 있다.
var1 = np.var(sample1, ddof=1)
var2 = np.var(sample2, ddof=1)

# 두 표본분산의 비. 귀무가설이 참이면 이 값이 1 근처에 놓인다.
f_statistic = var1 / var2

# 자유도
df1 = len(sample1) - 1
df2 = len(sample2) - 1

# F 분포는 대칭이 아니므로 작은 쪽 꼬리를 두 배 한다.
p_value = 2 * min(f.cdf(f_statistic, df1, df2),
                  f.sf(f_statistic, df1, df2))

print(f"Sample variances: {var1:.4f}, {var2:.4f}")
print(f"F-statistic: {f_statistic:.4f}")
print(f"Degrees of freedom: {df1}, {df2}")
print(f"Two-sided p-value: {p_value:.4f}")

# 결과 해석
alpha = 0.05
if p_value < alpha:
    print("Reject H0: variances are significantly different.")
else:
    print("Fail to reject H0: no significant difference in variances.")

# 세 절차의 실제 1종오류율을 잰다. H0 가 참이면 cdf(F) 는 U(0,1) 이므로
# F 를 직접 뽑지 않고 귀무분포에서 바로 뽑아도 같다.
rng = np.random.default_rng(11)
M, a = 400_000, 0.05
d = f(14, 19)                      # 보기 1 의 자유도로 잰다
U = d.cdf(d.rvs(M, random_state=rng))

print(f"\n명목 alpha = {a} 일 때 실제 1종오류율 (반복 {M:,}회, 자유도 14, 19)")
print(f"  (가) 방향을 미리 정한 단측  p = cdf(F)        : {np.mean(U < a):.4f}")
print(f"  (나) 자료를 보고 작은 꼬리  p = min(cdf, sf)  : {np.mean(np.minimum(U, 1 - U) < a):.4f}")
print(f"  (다) 등꼬리 양측           p = 2*min(cdf, sf): {np.mean(2 * np.minimum(U, 1 - U) < a):.4f}")

Sample variances: 2.8393, 6.0000
F-statistic: 0.4732
Degrees of freedom: 7, 7
Two-sided p-value: 0.3448
Fail to reject H0: no significant difference in variances.

명목 alpha = 0.05 일 때 실제 1종오류율 (반복 400,000회, 자유도 14, 19)
  (가) 방향을 미리 정한 단측  p = cdf(F)        : 0.0497
  (나) 자료를 보고 작은 꼬리  p = min(cdf, sf)  : 0.1001
  (다) 등꼬리 양측           p = 2*min(cdf, sf): 0.0498

import numpy as np
from scipy.stats import f

line1 = np.array([14.2, 13.8, 15.1, 14.7, 14.5, 14.6, 15.0, 14.8])
line2 = np.array([15.2, 14.8, 15.6, 15.0, 14.9, 15.3, 15.5, 14.9])

var1 = np.var(line1, ddof=1)
var2 = np.var(line2, ddof=1)
F_stat = var1 / var2
df1, df2 = len(line1) - 1, len(line2) - 1

alpha = 0.05
lower = f.ppf(alpha / 2, df1, df2)
upper = f.ppf(1 - alpha / 2, df1, df2)
p = 2 * min(f.cdf(F_stat, df1, df2), f.sf(F_stat, df1, df2))

print(f"variances: {var1:.5f}, {var2:.5f}")
print(f"F = {F_stat:.4f}, df = ({df1}, {df2})")
print(f"critical values: [{lower:.4f}, {upper:.4f}]")
print(f"two-sided p = {p:.4f}")

variances: 0.18125, 0.08857
F = 2.0464, df = (7, 7)
critical values: [0.2002, 4.9949]
two-sided p = 0.3655

from scipy.stats import f

print(f"{'n':>5} {'F_upper':>10} {'sd ratio':>10}")
for n in [5, 10, 20, 50, 100]:
    up = f.ppf(0.975, n - 1, n - 1)
    print(f"{n:>5} {up:>10.3f} {up**0.5:>10.3f}")

    n    F_upper   sd ratio
    5      9.605      3.099
   10      4.026      2.006
   20      2.526      1.589
   50      1.762      1.327
  100      1.486      1.219

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n, R, alpha = 30, 5000, 0.05

cases = [
    ("Normal",      lambda: rng.normal(0, 1, n)),
    ("t(10)",       lambda: rng.standard_t(10, n)),
    ("t(5)",        lambda: rng.standard_t(5, n)),
    ("Exponential", lambda: rng.exponential(1, n)),
]

for name, gen in cases:
    rej = 0
    for _ in range(R):
        F = gen().var(ddof=1) / gen().var(ddof=1)
        p = 2 * min(stats.f.cdf(F, n - 1, n - 1),
                    stats.f.sf(F, n - 1, n - 1))
        rej += (p < alpha)
    print(f"{name:>12}: empirical size = {rej / R:.4f}")

      Normal: empirical size = 0.0482
       t(10): empirical size = 0.0970
        t(5): empirical size = 0.1714
 Exponential: empirical size = 0.2768