F 검정 검정력 모의실험¶

import numpy as np
from scipy.stats import f

rng = np.random.default_rng(0)


def f_test_two_sided(x1, x2, alpha=0.05):
    """양측 F 검정이 기각하면 True 를 돌려준다."""
    n1, n2 = x1.size, x2.size
    df1, df2 = n1 - 1, n2 - 1
    F_obs = x1.var(ddof=1) / x2.var(ddof=1)
    p_left = f.cdf(F_obs, df1, df2)
    p_right = f.sf(F_obs, df1, df2)
    p_two = 2 * min(p_left, p_right)
    return p_two < alpha


def estimate_power(n1=12, n2=12, sigma1=1.0, sigma2=1.5,
                   n_sims=2000, alpha=0.05):
    """모의실험으로 양측 F 검정의 검정력을 추정한다.

    두 모집단의 표준편차를 실제로 다르게 두고, 검정이 그 차이를 몇 번이나
    잡아내는지 센다.
    """
    hits = 0
    for _ in range(n_sims):
        x1 = rng.normal(0, sigma1, size=n1)
        x2 = rng.normal(0, sigma2, size=n2)
        if f_test_two_sided(x1, x2, alpha=alpha):
            hits += 1
    return hits / n_sims


def exact_power(n1, n2, sigma1, sigma2, alpha=0.05):
    """(1) 에서 유도한 정확 검정력. k = sigma1^2/sigma2^2 이다."""
    df1, df2 = n1 - 1, n2 - 1
    lo = f.ppf(alpha / 2, df1, df2)
    hi = f.ppf(1 - alpha / 2, df1, df2)
    k = (sigma1 / sigma2) ** 2
    return f.cdf(lo / k, df1, df2) + f.sf(hi / k, df1, df2)


# 표준편차가 두 배나 차이 나는데도 n=10 이면 검정력이 절반에 못 미친다.
# 분산에 대한 검정이 평균에 대한 검정보다 훨씬 무딘 까닭은, 분산의
# 표집분포가 그만큼 넓게 퍼져 있기 때문이다.
power = estimate_power(n1=10, n2=10, sigma1=1.0, sigma2=2.0,
                       n_sims=5000, alpha=0.05)
se = np.sqrt(power * (1 - power) / 5000)
print(f"Estimated power: {power:.3f} (SE: {se:.3f})")

ex = exact_power(10, 10, 1.0, 2.0)
print(f"Exact power:     {ex:.4f}")
print(f"차이 = {abs(power - ex):.4f}   ( = {abs(power - ex) / se:.2f} 표준오차)")

Estimated power: 0.506 (SE: 0.007)
Exact power:     0.4964
차이 = 0.0094   ( = 1.33 표준오차)

# 이번에는 차이를 1.5 배로 줄이고 n 을 키워 간다. 쓸 만한 검정력을 얻으려면
# 표본이 얼마나 필요한지가 드러난다.
for n in [10, 20, 30, 50, 100]:
    pw = estimate_power(n1=n, n2=n, sigma1=1.0, sigma2=1.5, n_sims=3000)
    print(f"n1=n2={n:3d}: power = {pw:.3f}")

# (1) 의 정확 검정력과 나란히 놓는다. SE 는 sqrt(p(1-p)/3000) 이다.
import numpy as np
print(f"\n{'n':>5}{'모의':>9}{'정확':>9}{'차':>9}{'SE':>8}{'차/SE':>8}")
for n, pw in zip([10, 20, 30, 50, 100], [0.209, 0.395, 0.573, 0.809, 0.977]):
    ex = exact_power(n, n, 1.0, 1.5)
    se = np.sqrt(ex * (1 - ex) / 3000)
    print(f"{n:>5}{pw:>9.3f}{ex:>9.4f}{pw - ex:>9.4f}{se:>8.4f}{(pw - ex) / se:>8.2f}")

# 검정력 0.8 을 넘는 가장 작은 n
n = 5
while exact_power(n, n, 1.0, 1.5) < 0.8:
    n += 1
print(f"\n검정력 0.8 을 넘는 가장 작은 n = {n}  (검정력 {exact_power(n, n, 1.0, 1.5):.4f})")
print(f"  바로 앞 n = {n - 1} 의 검정력 = {exact_power(n - 1, n - 1, 1.0, 1.5):.4f}")

n1=n2= 10: power = 0.209
n1=n2= 20: power = 0.395
n1=n2= 30: power = 0.573
n1=n2= 50: power = 0.809
n1=n2=100: power = 0.977

    n       모의       정확        차      SE    차/SE
   10    0.209   0.2011   0.0079  0.0073    1.08
   20    0.395   0.4018  -0.0068  0.0090   -0.76
   30    0.573   0.5726   0.0004  0.0090    0.04
   50    0.809   0.8023   0.0067  0.0073    0.92
  100    0.977   0.9798  -0.0028  0.0026   -1.11

검정력 0.8 을 넘는 가장 작은 n = 50  (검정력 0.8023)
  바로 앞 n = 49 의 검정력 = 0.7940

import numpy as np
from scipy.stats import f

rng = np.random.default_rng(0)

def f_test_two_sided(x1, x2, alpha=0.05):
    F = x1.var(ddof=1) / x2.var(ddof=1)
    df1, df2 = x1.size - 1, x2.size - 1
    p = 2 * min(f.cdf(F, df1, df2), f.sf(F, df1, df2))
    return p < alpha

print(f"{'n':>5s}", end="")
for ratio in [1.25, 1.5, 2.0, 3.0]:
    print(f"  ratio={ratio:.2f}", end="")
print()

for n in [10, 25, 50, 100]:
    print(f"{n:5d}", end="")
    for ratio in [1.25, 1.5, 2.0, 3.0]:
        hits = 0
        for _ in range(5000):
            x1 = rng.normal(0, 1.0, n)
            x2 = rng.normal(0, ratio, n)
            if f_test_two_sided(x1, x2):
                hits += 1
        print(f"  {hits/5000:10.3f}", end="")
    print()

    n  ratio=1.25  ratio=1.50  ratio=2.00  ratio=3.00
   10       0.096       0.205       0.500       0.873
   25       0.184       0.487       0.910       1.000
   50       0.338       0.801       0.998       1.000
  100       0.598       0.981       1.000       1.000

import numpy as np
from scipy.stats import levene, f as fdist

rng = np.random.default_rng(0)

for n in [10, 25, 50, 100]:
    for ratio in [1.5, 2.0]:
        hits_f, hits_l = 0, 0
        for _ in range(3000):
            x1 = rng.normal(0, 1.0, n)
            x2 = rng.normal(0, ratio, n)
            F = np.var(x1, ddof=1) / np.var(x2, ddof=1)
            p = 2 * min(fdist(n-1, n-1).cdf(F), fdist(n-1, n-1).sf(F))
            if p < 0.05:
                hits_f += 1
            _, p_l = levene(x1, x2, center='median')
            if p_l < 0.05:
                hits_l += 1
        print(f"n={n:3d}, ratio={ratio}: "
              f"F-test={hits_f/3000:.3f}, Levene={hits_l/3000:.3f}")

n= 10, ratio=1.5: F-test=0.218, Levene=0.131
n= 10, ratio=2.0: F-test=0.508, Levene=0.311
n= 25, ratio=1.5: F-test=0.487, Levene=0.384
n= 25, ratio=2.0: F-test=0.903, Levene=0.812
n= 50, ratio=1.5: F-test=0.812, Levene=0.725
n= 50, ratio=2.0: F-test=0.997, Levene=0.992
n=100, ratio=1.5: F-test=0.976, Levene=0.958
n=100, ratio=2.0: F-test=1.000, Levene=1.000

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import f as fdist

rng = np.random.default_rng(0)
n, n_sims, alpha = 20, 3000, 0.05
ratios = np.arange(1.0, 3.05, 0.1)
powers = []

for r in ratios:
    hits = 0
    for _ in range(n_sims):
        x1 = rng.normal(0, 1.0, n)
        x2 = rng.normal(0, r, n)
        F = np.var(x1, ddof=1) / np.var(x2, ddof=1)
        p = 2 * min(fdist(n-1, n-1).cdf(F), fdist(n-1, n-1).sf(F))
        if p < alpha:
            hits += 1
    powers.append(hits / n_sims)

plt.figure(figsize=(8, 4))
plt.plot(ratios, powers, "o-", markersize=4)
plt.axhline(alpha, ls="--", color="red", label=f"alpha = {alpha}")
plt.axhline(0.8, ls=":", color="gray", label="power = 0.8")
plt.xlabel("sigma2 / sigma1")
plt.ylabel("Estimated power")
plt.title("F-test power curve (n1 = n2 = 20)")
plt.legend()
plt.tight_layout()
plt.show()

import numpy as np
from scipy.stats import f as fdist

rng = np.random.default_rng(0)
n_sims = 5000

for n1, n2 in [(5, 100), (50, 50)]:
    hits = 0
    for _ in range(n_sims):
        x1 = rng.normal(0, 1.0, n1)
        x2 = rng.normal(0, 2.0, n2)
        F = np.var(x1, ddof=1) / np.var(x2, ddof=1)
        p = 2 * min(fdist(n1-1, n2-1).cdf(F), fdist(n1-1, n2-1).sf(F))
        if p < 0.05:
            hits += 1
    print(f"n1={n1:3d}, n2={n2:3d}: power = {hits/n_sims:.3f}")

n1=  5, n2=100: power = 0.257
n1= 50, n2= 50: power = 0.997