검정력 분석과 표본크기¶

import numpy as np
from scipy import stats

def power_ttest(n, delta, sigma=1.0, alpha=0.05):
    """양측 이표본 t-검정의 검정력 (정규근사)."""
    # 두 집단의 차이라 분산이 두 번 들어간다. 그래서 sqrt(2/n)이다.
    se = sigma * np.sqrt(2 / n)
    z_crit = stats.norm.ppf(1 - alpha / 2)
    z_effect = delta / se
    # 두 꼬리를 모두 세는 것이 정확하다. 두 번째 항은 참 효과가 양수인데도
    # 통계량이 반대쪽 꼬리로 넘어가 기각되는 경우이며, 보통 무시할 만큼 작다.
    power = (1 - stats.norm.cdf(z_crit - z_effect)
             + stats.norm.cdf(-z_crit - z_effect))
    return power

for n in [20, 50, 100]:
    print(f"n={n:>4} per group: power = {power_ttest(n, delta=0.5):.4f}")

# (1) 두 항을 따로 본다.
from statsmodels.stats.power import TTestIndPower
zc = stats.norm.ppf(0.975)
print(f"\n{'n':>5} {'lambda':>9} {'첫째 항':>10} {'둘째 항':>11}"
      f" {'정규근사':>9} {'비중심 t':>9} {'차이':>8}")
for n in (20, 50, 100):
    lam = 0.5 * np.sqrt(n / 2)
    up, lo = stats.norm.sf(zc - lam), stats.norm.cdf(-zc - lam)
    nct = TTestIndPower().power(effect_size=0.5, nobs1=n, alpha=0.05,
                                ratio=1.0, alternative='two-sided')
    print(f"{n:>5} {lam:>9.6f} {up:>10.6f} {lo:>11.3e}"
          f" {up + lo:>9.4f} {nct:>9.4f} {up + lo - nct:>+8.4f}")

n=  20 per group: power = 0.3526
n=  50 per group: power = 0.7054
n= 100 per group: power = 0.9424

    n    lambda       첫째 항        둘째 항      정규근사     비중심 t       차이
   20  1.581139   0.352409   1.992e-04    0.3526    0.3379  +0.0147
   50  2.500000   0.705414   4.099e-06    0.7054    0.6969  +0.0085
  100  3.535534   0.942438   1.948e-08    0.9424    0.9404  +0.0020

def sample_size_ttest(delta, sigma=1.0, alpha=0.05, power=0.80):
    """이표본 t-검정의 집단당 최소 표본크기."""
    z_alpha = stats.norm.ppf(1 - alpha / 2)
    z_beta = stats.norm.ppf(power)
    # 앞의 계수 2가 "두 집단"의 대가다. 일표본 공식과 여기서 갈린다.
    n = 2 * ((z_alpha + z_beta) * sigma / delta) ** 2
    return int(np.ceil(n))

def sample_size_proportion(p1, p2, alpha=0.05, power=0.80):
    """이표본 비율검정의 집단당 최소 표본크기."""
    p_bar = (p1 + p2) / 2
    z_alpha = stats.norm.ppf(1 - alpha / 2)
    z_beta = stats.norm.ppf(power)
    # 두 항의 제곱근 안이 서로 다르다는 점이 핵심이다.
    #   z_alpha 쪽: H0("두 비율이 같다") 아래의 분산이므로 합동비율 p_bar를 쓴다.
    #   z_beta  쪽: H1 아래의 분산이므로 p1과 p2를 각각 쓴다.
    # 검정력 계산은 두 가설 아래의 분포를 동시에 다루므로 이렇게 섞인다.
    numer = (z_alpha * np.sqrt(2 * p_bar * (1 - p_bar))
             + z_beta * np.sqrt(p1 * (1 - p1) + p2 * (1 - p2))) ** 2
    n = numer / (p1 - p2) ** 2
    return int(np.ceil(n))

# (2) 두 제곱근을 하나로 통일하면 어떻게 되는가.
za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
for p1, p2 in [(0.0121, 0.011), (0.6, 0.3)]:
    pb = (p1 + p2) / 2
    s0 = np.sqrt(2 * pb * (1 - pb))                   # H0 아래
    s1 = np.sqrt(p1 * (1 - p1) + p2 * (1 - p2))       # H1 아래
    ok = (za * s0 + zb * s1) ** 2 / (p1 - p2) ** 2
    only0 = ((za + zb) * s0) ** 2 / (p1 - p2) ** 2
    only1 = ((za + zb) * s1) ** 2 / (p1 - p2) ** 2
    print(f"p1={p1}, p2={p2}:  s0={s0:.6f}  s1={s1:.6f}  s0/s1={s0 / s1:.4f}")
    print(f"   옳은 식 {ok:>12.3f}   s0 만 {only0:>12.3f}"
          f"   s1 만 {only1:>12.3f}")
    print(f"   옳은 n 에서의 검정력 "
          f"{stats.norm.sf((za * s0 - (p1 - p2) * np.sqrt(np.ceil(ok))) / s1):.4f}")

p1=0.0121, p2=0.011:  s0=0.151107  s1=0.151105  s0/s1=1.0000
   옳은 식   148110.393   s0 만   148111.571   s1 만   148107.647
   옳은 n 에서의 검정력 0.8000
p1=0.6, p2=0.3:  s0=0.703562  s1=0.670820  s0/s1=1.0488
   옳은 식       41.970   s0 만       43.169   s1 만       39.244
   옳은 n 에서의 검정력 0.8003

# 이표본 t-검정: 중간 크기 효과 (Cohen's d = 0.5)
delta = 0.5
n_req = sample_size_ttest(delta, sigma=1.0, alpha=0.05, power=0.80)
print(f"Required n per group: {n_req}")

# 비율 검정 (A/B 검정): 전환율 1.10% -> 1.21%, 즉 10% 상대 개선
p1, p2 = 0.0121, 0.011
n_prop = sample_size_proportion(p1, p2, alpha=0.05, power=0.80)
print(f"Required n per group: {n_prop:,}")

# (2) 유효 효과크기로 읽어 본다.
pb = (p1 + p2) / 2
d_eff = (p1 - p2) / np.sqrt(pb * (1 - pb))
print(f"\n유효 효과크기 d_eff = {d_eff:.6f}"
      f"   (t 검정의 d = 0.5 의 1/{0.5 / d_eff:.1f})")
print(f"d_eff 로 예측한 n   = {2 * ((za + zb) / d_eff) ** 2:>12.1f}")
print(f"공식이 준 n         = {n_prop:>12,}")
print(f"비: 예측 {(0.5 / d_eff) ** 2:.1f}배   실제 {n_prop / n_req:.1f}배")
print(f"\n필요한 방문자: 집단당 {n_prop:,}명,  합계 {2 * n_prop:,}명")

Required n per group: 63
Required n per group: 148,111

유효 효과크기 d_eff = 0.010295   (t 검정의 d = 0.5 의 1/48.6)
d_eff 로 예측한 n   =     148111.6
공식이 준 n         =      148,111
비: 예측 2358.8배   실제 2351.0배

필요한 방문자: 집단당 148,111명,  합계 296,222명

import matplotlib.pyplot as plt

# 효과가 작을수록 같은 검정력에 필요한 표본이 가파르게 늘어난다.
# d=0.2 곡선이 0.8 에 닿는 자리를 d=0.8 곡선의 그것과 견주어 보면 된다.
ns = np.arange(10, 500)
fig, ax = plt.subplots(figsize=(10, 5))
for d, ls in [(0.2, '--'), (0.5, '-'), (0.8, ':')]:
    powers = [power_ttest(n, d) for n in ns]
    ax.plot(ns, powers, ls, label=f'd = {d}')

ax.axhline(0.80, color='grey', linestyle='-.', alpha=0.5, label='Power = 0.80')
ax.set_xlabel('Sample size per group (n)')
ax.set_ylabel('Power')
ax.set_title('Power Curves for Two-Sample t-Test')
ax.legend()
plt.tight_layout()
plt.show()

# (1) 의 예측과 곡선 위에서 찾은 자리를 맞춘다.
print(f"{'d':>5} {'공식':>10} {'곡선':>6} {'그 n 의 검정력':>14} {'직전 n':>7}")
for d in (0.2, 0.5, 0.8):
    pw = np.array([power_ttest(n, d) for n in ns])
    hit = np.flatnonzero(pw >= 0.80)
    # **조건을 만족하는 칸이 있는지 먼저 확인한다.** 없으면 아래 인덱싱이
    # 조용히 엉뚱한 답을 주는 대신 여기서 멈춘다.
    assert hit.size > 0, f"d={d} 는 격자 안에서 0.80 에 닿지 않는다"
    k = hit[0]
    print(f"{d:>5.1f} {2 * ((za + zb) / d) ** 2:>10.3f} {ns[k]:>6}"
          f" {pw[k]:>14.4f} {pw[k - 1]:>7.4f}")

# 평평해지는 구간: d = 0.5 에서 검정력을 더 올리는 값.
print()
for target, zq in ((0.80, zb), (0.90, stats.norm.ppf(0.90)),
                   (0.95, stats.norm.ppf(0.95))):
    n = 2 * ((za + zq) / 0.5) ** 2
    print(f"  d=0.5, 검정력 {target:.2f}: n = {n:7.2f}"
          f" -> {int(np.ceil(n)):>3}  (실제 {power_ttest(int(np.ceil(n)), 0.5):.4f})")
print(f"  d=0.2 는 격자 끝 n=499 에서도 {power_ttest(499, 0.2):.4f} 까지만 간다")

    d         공식     곡선      그 n 의 검정력    직전 n
  0.2    392.444    393         0.8006  0.7996
  0.5     62.791     63         0.8013  0.7950
  0.8     24.528     25         0.8074  0.7914

  d=0.5, 검정력 0.80: n =   62.79 ->  63  (실제 0.8013)
  d=0.5, 검정력 0.90: n =   84.06 ->  85  (실제 0.9031)
  d=0.5, 검정력 0.95: n =  103.96 -> 104  (실제 0.9501)
  d=0.2 는 격자 끝 n=499 에서도 0.8848 까지만 간다

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.stats.power import TTestPower

analysis = TTestPower()
ns = np.arange(5, 201)

fig, ax = plt.subplots(figsize=(9, 5))
for d in [0.2, 0.5, 0.8]:
    powers = [analysis.power(effect_size=d, nobs=n, alpha=0.05) for n in ns]
    ax.plot(ns, powers, label=f'd = {d}')

ax.axhline(0.80, color='grey', linestyle='--', label='80% power')
ax.set_xlabel('Sample size n')
ax.set_ylabel('Power')
ax.set_title('Power Curves (One-Sample t-Test)')
ax.legend()
plt.tight_layout()
plt.show()

import numpy as np
from scipy import stats
from scipy.optimize import brentq

rng = np.random.default_rng(6)
M = 400_000
za = stats.norm.ppf(0.975)

print(f"{'검정력':>8s} {'θ/SE':>7s} {'과장 비 M':>10s} {'부호 오류 S':>12s}")
for pw in [0.06, 0.10, 0.20, 0.50, 0.80, 0.95]:
    # 목표 검정력을 주는 참 효과(표준오차 단위)
    nc = brentq(lambda c: stats.norm.sf(za - c)
                + stats.norm.cdf(-za - c) - pw, 1e-9, 10)
    z = rng.normal(nc, 1, M)
    est = z[np.abs(z) > za]                      # 유의한 결과만
    print(f"{pw:8.0%} {nc:7.3f} {np.mean(np.abs(est)) / nc:10.3f} "
          f"{np.mean(est < 0):12.4f}")

     검정력    θ/SE     과장 비 M      부호 오류 S
      6%   0.295      8.015       0.2030
     10%   0.652      3.710       0.0449
     20%   1.115      2.258       0.0051
     50%   1.960      1.407       0.0001
     80%   2.802      1.125       0.0000
     95%   3.605      1.030       0.0000

import numpy as np
from scipy import stats

za = stats.norm.ppf(0.975)
deltas = np.array([0, 0.2, 0.5, 0.8, 1.0, 1.5])

def power(d, n, alpha=0.05):
    z = stats.norm.ppf(1 - alpha / 2)
    nc = d * np.sqrt(n)
    return stats.norm.sf(z - nc) + stats.norm.cdf(-z - nc)

print(f"{'δ':>5s} " + " ".join(f"{'n='+str(n):>9s}"
                               for n in [10, 25, 50, 100]))
for d in deltas:
    print(f"{d:5.1f} " + " ".join(f"{power(d, n):9.4f}"
                                  for n in [10, 25, 50, 100]))

print("\nα 를 바꾸면 (n = 25)")
print(f"{'δ':>5s} " + " ".join(f"{'α='+str(a):>9s}"
                               for a in [0.10, 0.05, 0.01, 0.001]))
for d in deltas:
    print(f"{d:5.1f} " + " ".join(f"{power(d, 25, a):9.4f}"
                                  for a in [0.10, 0.05, 0.01, 0.001]))

    δ      n=10      n=25      n=50     n=100
  0.0    0.0500    0.0500    0.0500    0.0500
  0.2    0.0969    0.1701    0.2930    0.5160
  0.5    0.3526    0.7054    0.9424    0.9988
  0.8    0.7156    0.9793    0.9999    1.0000
  1.0    0.8854    0.9988    1.0000    1.0000
  1.5    0.9973    1.0000    1.0000    1.0000

α 를 바꾸면 (n = 25)
    δ     α=0.1    α=0.05    α=0.01   α=0.001
  0.0    0.1000    0.0500    0.0100    0.0010
  0.2    0.2636    0.1701    0.0577    0.0110
  0.5    0.8038    0.7054    0.4698    0.2146
  0.8    0.9907    0.9793    0.9228    0.7610
  1.0    0.9996    0.9988    0.9923    0.9563
  1.5    1.0000    1.0000    1.0000    1.0000

import numpy as np
from scipy import stats

def power(delta, n, sigma=1.0, alpha=0.05, rho=0.0, paired=False):
    se = sigma * np.sqrt(2 * (1 - rho) / n) if paired else sigma * np.sqrt(2 / n)
    z = stats.norm.ppf(1 - alpha / 2)
    nc = delta / se
    return stats.norm.sf(z - nc) + stats.norm.cdf(-z - nc)

base = power(0.5, 30)
print(f"기준: δ=0.5, σ=1, 집단당 n=30, α=0.05 → 검정력 {base:.4f}\n")
cases = [
    ("n 을 2배(60)로", power(0.5, 60)),
    ("n 을 4배(120)로", power(0.5, 120)),
    ("σ 를 20% 줄임", power(0.5, 30, sigma=0.8)),
    ("σ 를 절반으로", power(0.5, 30, sigma=0.5)),
    ("δ 를 1.5배로(설계 강화)", power(0.75, 30)),
    ("α 를 0.10으로", power(0.5, 30, alpha=0.10)),
    ("대응설계 ρ=0.5", power(0.5, 30, paired=True, rho=0.5)),
    ("대응설계 ρ=0.8", power(0.5, 30, paired=True, rho=0.8)),
]
for name, v in cases:
    print(f"{name:24s} {v:.4f}   (+{v - base:.4f})")

기준: δ=0.5, σ=1, 집단당 n=30, α=0.05 → 검정력 0.4907

n 을 2배(60)로              0.7819   (+0.2912)
n 을 4배(120)로             0.9721   (+0.4814)
σ 를 20% 줄임               0.6775   (+0.1868)
σ 를 절반으로                 0.9721   (+0.4814)
δ 를 1.5배로(설계 강화)         0.8276   (+0.3369)
α 를 0.10으로               0.6149   (+0.1242)
대응설계 ρ=0.5               0.7819   (+0.2912)
대응설계 ρ=0.8               0.9911   (+0.5004)

import numpy as np
from scipy import stats

za = stats.norm.ppf(0.975)
t = 0.5                                  # 정보 비율(절반 모음)

def cp(z_now, delta_over_se_final, t=0.5, za=za):
    """z_now: 현재 시점의 z 통계량 (정보비율 t)"""
    num = za - np.sqrt(t) * z_now - (1 - t) * delta_over_se_final
    return stats.norm.sf(num / np.sqrt(1 - t))

print(f"{'현재 z':>8s} {'δ=계획값':>10s} {'δ=현재추정':>12s} {'δ=0':>8s}")
for z_now in [0.0, 0.5, 1.0, 1.5, 2.0]:
    planned = za + stats.norm.ppf(0.90)          # 검정력 90% 설계
    print(f"{z_now:8.1f} {cp(z_now, planned):10.4f} "
          f"{cp(z_now, z_now / np.sqrt(t)):12.4f} {cp(z_now, 0.0):8.4f}")

    현재 z      δ=계획값       δ=현재추정      δ=0
     0.0     0.3157       0.0028   0.0028
     0.5     0.5081       0.0382   0.0115
     1.0     0.6986       0.2201   0.0382
     1.5     0.8462       0.5903   0.1017
     2.0     0.9358       0.8903   0.2201

import numpy as np
from scipy import stats

rng = np.random.default_rng(101)
M = 20_000

def sim_power(gen, analyze, M=M, alpha=0.05):
    return np.mean([analyze(*gen()) < alpha for _ in range(M)])

# 1) 공식이 있는 경우로 검증
n, d = 30, 0.5
gen1 = lambda: (rng.normal(d, 1, n), rng.normal(0, 1, n))
ana1 = lambda x, y: stats.ttest_ind(x, y, equal_var=False).pvalue
theo = stats.norm.sf(stats.norm.ppf(0.975) - d / np.sqrt(2 / n))
print(f"t 검정  모의 {sim_power(gen1, ana1):.4f}   이론(정규근사) {theo:.4f}")

# 2) 공식이 마땅치 않은 경우: 치우친 자료에 만-휘트니
def gen2():
    x = rng.lognormal(0.3, 1.0, n)
    y = rng.lognormal(0.0, 1.0, n)
    return x, y
ana2 = lambda x, y: stats.mannwhitneyu(x, y, alternative="two-sided").pvalue
ana3 = lambda x, y: stats.ttest_ind(x, y, equal_var=False).pvalue
print(f"로그정규 — 만-휘트니 {sim_power(gen2, ana2):.4f}   "
      f"웰치 t {sim_power(gen2, ana3):.4f}")

t 검정  모의 0.4768   이론(정규근사) 0.4906
로그정규 — 만-휘트니 0.1961   웰치 t 0.1422