검정력 분석¶

from scipy import stats
import numpy as np

def sample_size_z_test(effect_size, alpha=0.05, power=0.80, alternative='two-sided'):
    """일표본 z-검정에 필요한 표본크기."""
    if alternative == 'two-sided':
        z_alpha = stats.norm.ppf(1 - alpha / 2)
    else:
        z_alpha = stats.norm.ppf(1 - alpha)
    # z_beta는 ppf(power)이지 ppf(1 - power)가 아니다.
    # 검정력이 클수록 z_beta가 커져 n이 늘어나야 하기 때문이다.
    z_beta = stats.norm.ppf(power)
    # 두 임계값을 **더한다**. alpha는 H0 분포의 오른쪽 꼬리에서,
    # beta는 Ha 분포의 왼쪽 꼬리에서 재므로 둘 사이의 거리가 두 값의 합이다.
    n = ((z_alpha + z_beta) / effect_size) ** 2
    return int(np.ceil(n))

# 효과크기 0.5를 검정력 80%로 탐지하려면?
n = sample_size_z_test(effect_size=0.5)
print(f"Required sample size: {n}")

# (2) 올림하기 전의 값과 1/d^2 비례를 확인한다.
za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
print(f"za + zb = {za + zb:.6f}")
for d in (0.5, 0.25):
    print(f"  d = {d:.2f}: n = {((za + zb) / d) ** 2:.4f}"
          f"  -> {sample_size_z_test(effect_size=d)}")

Required sample size: 32
za + zb = 2.801585
  d = 0.50: n = 31.3955  -> 32
  d = 0.25: n = 125.5821  -> 126

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

def plot_power(mu_0, mu_a, sigma, n, alpha=0.05):
    """단측 z 검정의 검정력을 그림으로 보인다.

        귀무분포와 대립분포를 겹쳐 그리고, 기각역에 해당하는 두 넓이를 칠한다.
        귀무 쪽 넓이가 alpha, 대립 쪽 넓이가 검정력이다.
        """
    se = sigma / np.sqrt(n)
    z_crit = stats.norm.ppf(1 - alpha)
    x_crit = mu_0 + z_crit * se

    x = np.linspace(mu_0 - 4*se, mu_a + 4*se, 300)

    # 귀무가설이 참일 때의 분포
    y_h0 = stats.norm(mu_0, se).pdf(x)
    # 대립가설이 참일 때의 분포. 중심만 오른쪽으로 옮겨 간다.
    y_ha = stats.norm(mu_a, se).pdf(x)

    fig, ax = plt.subplots(figsize=(12, 4))
    ax.plot(x, y_h0, 'b-', label=f'$H_0$: $\\mu = {mu_0}$')
    ax.plot(x, y_ha, 'r-', label=f'$H_a$: $\\mu = {mu_a}$')

    # 귀무 분포에서 기각역에 해당하는 꼬리. 그 넓이가 유의수준 alpha 다.
    x_reject = x[x >= x_crit]
    ax.fill_between(x_reject, stats.norm(mu_0, se).pdf(x_reject), alpha=0.3, color='blue', label=f'$\\alpha$ = {alpha}')

    # 대립 분포에서 같은 기각역의 넓이. 그것이 검정력이다.
    ax.fill_between(x_reject, stats.norm(mu_a, se).pdf(x_reject), alpha=0.3, color='red', label=f'Power = {1 - stats.norm(mu_a, se).cdf(x_crit):.3f}')

    ax.axvline(x_crit, color='k', linestyle='--', label=f'Critical value = {x_crit:.2f}')
    ax.legend()
    ax.set_xlabel('Sample Mean')
    ax.set_ylabel('Density')
    ax.set_title('Power of a Hypothesis Test')
    plt.show()

plot_power(mu_0=50, mu_a=52, sigma=10, n=25)

# (1) 과 (2) 의 닫힌 꼴을 확인한다.
mu0, mua, sigma, n_now, al = 50, 52, 10, 25, 0.05
se = sigma / np.sqrt(n_now)
z_al, z_be = stats.norm.ppf(1 - al), stats.norm.ppf(0.80)
x_c = mu0 + z_al * se
print(f"SE = {se}, x_c = {x_c:.6f}")
print(f"검정력  닫힌 꼴 {stats.norm.sf(z_al - (mua - mu0) / se):.6f}"
      f"   직접 {stats.norm(mua, se).sf(x_c):.6f}")

d = (mua - mu0) / sigma
n_need = ((z_al + z_be) / d) ** 2
print(f"검정력 0.80 에 필요한 n = {n_need:.2f} -> {int(np.ceil(n_need))}")
se2 = sigma / np.sqrt(np.ceil(n_need))
print(f"  그때 SE = {se2:.4f} (지금의 1/{se / se2:.2f}),"
      f"  x_c = {mu0 + z_al * se2:.2f},"
      f"  검정력 = {stats.norm.sf(z_al - (mua - mu0) / se2):.4f}")

SE = 2.0, x_c = 53.289707
검정력  닫힌 꼴 0.259511   직접 0.259511
검정력 0.80 에 필요한 n = 154.56 -> 155
  그때 SE = 0.8032 (지금의 1/2.49),  x_c = 51.32,  검정력 = 0.8010

from statsmodels.stats.power import TTestPower
import numpy as np

analysis = TTestPower()

# 5점 차이를 탐지하려면 몇 명이 필요한가? (mu_0 = 0, mu_a = 5, sigma = 15)
# 검정력 계산에 들어가는 것은 delta도 sigma도 아니라 그 비뿐이다.
# 5점/15와 1점/3은 검정력 관점에서 완전히 같은 문제다.
effect_size = 5 / 15  # Cohen's d = 0.333

n_needed = analysis.solve_power(effect_size=effect_size, alpha=0.05,
                                power=0.80, alternative='two-sided')
print(f"One-sample t-test:")
print(f"  Effect size (Cohen's d): {effect_size:.3f}")
print(f"  Sample size needed for 80% power: {int(np.ceil(n_needed))}")

# 반대 방향의 질문: n이 정해져 있을 때 검정력은 얼마인가?
power = analysis.power(effect_size=effect_size, nobs=50, alpha=0.05,
                      alternative='two-sided')
print(f"  Power with n=50: {power:.3f}")

# (1) 의 z 어림과 견준다.
za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
n_z = ((za + zb) / effect_size) ** 2
pw_z = stats.norm.sf(za - effect_size * np.sqrt(50))
print(f"\n  z 어림   n = {n_z:.4f},  n=50 의 검정력 = {pw_z:.6f}")
print(f"  비중심 t n = {n_needed:.4f},  n=50 의 검정력 = {power:.6f}")
print(f"  차이        {n_needed - n_z:+.4f}명,"
      f"             {power - pw_z:+.6f}")

One-sample t-test:
  Effect size (Cohen's d): 0.333
  Sample size needed for 80% power: 73
  Power with n=50: 0.637

  z 어림   n = 70.6399,  n=50 의 검정력 = 0.654338
  비중심 t n = 72.5839,  n=50 의 검정력 = 0.637094
  차이        +1.9440명,             -0.017244

from statsmodels.stats.power import TTestIndPower

# 두 집단을 견주는 t-검정의 검정력 분석
analysis = TTestIndPower()

# 두 집단의 크기를 같게 두는 설계다. 같은 총인원이면 이때 검정력이 가장 높다.
# 효과크기 Cohen 의 d = 0.5. 두 평균이 표준편차의 절반만큼 떨어진 경우다.
effect_size = 0.5

n_per_group = analysis.solve_power(effect_size=effect_size, alpha=0.05,
                                   power=0.80, ratio=1.0,
                                   alternative='two-sided')
print(f"\nTwo-sample t-test (equal n):")
print(f"  Effect size (Cohen's d): {effect_size:.3f}")
print(f"  Sample size per group for 80% power: {int(np.ceil(n_per_group))}")
print(f"  Total sample size: {2 * int(np.ceil(n_per_group))}")

# 집단 크기가 다른 경우 (예: 2:1 배분)
# ratio는 nobs2/nobs1이고, solve_power가 돌려주는 것은 **nobs1**이다.
# 따라서 두 번째 집단은 ratio를 곱해서 얻어야 한다. 나누면 거꾸로다.
ratio = 2.0
n_group1 = analysis.solve_power(effect_size=effect_size, alpha=0.05,
                                power=0.80, ratio=ratio,
                                alternative='two-sided')
n_group2 = n_group1 * ratio
n1, n2 = int(np.ceil(n_group1)), int(np.ceil(n_group2))
print(f"\nTwo-sample t-test (2:1 ratio):")
print(f"  Group 1 n: {n1}")
print(f"  Group 2 n: {n2}")
print(f"  Total sample size: {n1 + n2}")

# (1) 의 식을 여러 r 에서 확인한다.
print(f"\n{'r':>5} {'n1':>8} {'n2':>8} {'총':>6} {'관측 비':>8} {'식':>8}")
base = None
for r in (1.0, 1.5, 2.0, 3.0, 4.0):
    a = analysis.solve_power(effect_size=0.5, alpha=0.05, power=0.80,
                             ratio=r, alternative='two-sided')
    total = a * (1 + r)
    base = total if base is None else base
    print(f"{r:>5.1f} {a:>8.2f} {a * r:>8.2f} {total:>6.1f}"
          f" {total / base:>8.4f} {(2 + r + 1 / r) / 4:>8.4f}")

Two-sample t-test (equal n):
  Effect size (Cohen's d): 0.500
  Sample size per group for 80% power: 64
  Total sample size: 128

Two-sample t-test (2:1 ratio):
  Group 1 n: 48
  Group 2 n: 96
  Total sample size: 144

    r       n1       n2      총     관측 비        식
  1.0    63.77    63.77  127.5   1.0000   1.0000
  1.5    53.11    79.66  132.8   1.0410   1.0417
  2.0    47.74    95.48  143.2   1.1231   1.1250
  3.0    42.35   127.04  169.4   1.3282   1.3333
  4.0    39.63   158.53  198.2   1.5538   1.5625

import numpy as np
import statsmodels.stats.api as sms
from statsmodels.stats.power import NormalIndPower

# 전환율을 견주는 A/B 검정
# 대조군 전환율 1.1%
# 실험군 전환율 1.65%. 상대적으로는 50% 개선이지만 절대차는 0.55%p 다.
p0 = 0.011   # Control baseline
p1 = 0.0165  # Treatment goal

# 비율에서는 두 값의 차이 대신 arcsin 변환 후의 차이를 효과크기로 쓴다.
#   h = 2*arcsin(sqrt(p1)) - 2*arcsin(sqrt(p0)) — 비율을 각도로 바꿔 재는 효과크기다.
# 이렇게 하면 분산이 p에 의존하는 문제가 사라져 하나의 공식으로 처리된다.
# 0.011 대 0.0165는 절대차로 0.55%p뿐이지만 상대적으로는 50% 증가다.
effect_size = sms.proportion_effectsize(p1, p0)

# alternative='larger'는 단측이다. A/B 검정에서는 "개선되었는가"만 보는 일이 많다.
analysis = NormalIndPower()
n_ab = analysis.solve_power(effect_size=effect_size, alpha=0.05,
                            power=0.80, ratio=1.0, alternative='larger')

print(f"\nA/B Test (Proportions):")
print(f"  Control rate: {p0:.2%}")
print(f"  Treatment goal: {p1:.2%}")
print(f"  Effect size: {effect_size:.4f}")
print(f"  Sample size per group for 80% power: {int(np.ceil(n_ab))}")

# (1) 의 손계산과 맞춰 본다.
h = 2 * np.arcsin(np.sqrt(p1)) - 2 * np.arcsin(np.sqrt(p0))
z_a, z_b = stats.norm.ppf(0.95), stats.norm.ppf(0.80)
n_formula = 2 * ((z_a + z_b) / h) ** 2
print(f"\n  h 직접 계산   {h:.8f}   statsmodels {effect_size:.8f}")
print(f"  n 공식        {n_formula:.4f}   statsmodels {n_ab:.4f}")

# (2) 그 표본에서 전환은 몇 건인가.
n_int = int(np.ceil(n_ab))
print(f"\n  집단당 {n_int:,}명 x 2 = {2 * n_int:,}명")
print(f"  기대 전환 건수: 대조군 {n_int * p0:.1f}건,"
      f" 실험군 {n_int * p1:.1f}건,  차이 {n_int * (p1 - p0):.1f}건")

A/B Test (Proportions):
  Control rate: 1.10%
  Treatment goal: 1.65%
  Effect size: 0.0475
  Sample size per group for 80% power: 5488

  h 직접 계산   0.04746819   statsmodels 0.04746819
  n 공식        5487.7312   statsmodels 5487.7312

  집단당 5,488명 x 2 = 10,976명
  기대 전환 건수: 대조군 60.4건, 실험군 90.6건,  차이 30.2건

from statsmodels.stats.power import FTestAnovaPower

# 분산분석의 효과크기는 Cohen 의 f 다. t-검정의 d 와는 눈금이 다르므로
# 0.25 를 d 로 읽으면 안 된다. f 는 0.10 작음, 0.25 중간, 0.40 큼으로 본다.
analysis = FTestAnovaPower()

effect_size = 0.25   # 중간 크기 효과
k_groups = 4         # 비교할 집단 수

# **nobs 는 총 관측수다.** 집단당이 아니다. 아래 k=2 검산이 그것을 보인다.
n_total = analysis.solve_power(effect_size=effect_size, alpha=0.05,
                               power=0.80, k_groups=k_groups)
n_each = int(np.ceil(n_total / k_groups))
print(f"\nOne-way ANOVA (4 groups):")
print(f"  Effect size (Cohen's f): {effect_size:.3f}")
print(f"  Total sample size for 80% power: {n_total:.3f}"
      f" -> {k_groups * n_each}")
print(f"  Sample size per group: {n_total / k_groups:.3f} -> {n_each}")
print(f"  실제 검정력: "
      f"{analysis.power(effect_size=effect_size, nobs=k_groups * n_each, alpha=0.05, k_groups=k_groups):.4f}")

# (1) 의 검산: k=2 는 d=0.5 인 이표본 t 검정과 같은 문제여야 한다.
n_k2 = analysis.solve_power(effect_size=0.25, alpha=0.05, power=0.80,
                            k_groups=2)
n_t = TTestIndPower().solve_power(effect_size=0.5, alpha=0.05, power=0.80,
                                  ratio=1.0, alternative='two-sided')
print(f"\n  ANOVA k=2 가 준 수      {n_k2:.5f}")
print(f"  이표본 t 의 집단당       {n_t:.5f}   총 {2 * n_t:.5f}")
print(f"  -> 돌려준 수는 '총'이다 (차이 {abs(n_k2 - 2 * n_t):.2e})")

# 집단 수를 늘리면 어떻게 되는가.
print(f"\n{'k':>3} {'총 N':>9} {'집단당':>8} {'올림 후 총':>10}")
for k in (2, 3, 4, 5, 6):
    N = analysis.solve_power(effect_size=0.25, alpha=0.05, power=0.80,
                             k_groups=k)
    each = int(np.ceil(N / k))
    print(f"{k:>3} {N:>9.3f} {N / k:>8.3f} {k * each:>10}")

One-way ANOVA (4 groups):
  Effect size (Cohen's f): 0.250
  Total sample size for 80% power: 178.396 -> 180
  Sample size per group: 44.599 -> 45
  실제 검정력: 0.8040

  ANOVA k=2 가 준 수      127.53062
  이표본 t 의 집단당       63.76561   총 127.53122
  -> 돌려준 수는 '총'이다 (차이 6.00e-04)

  k       총 N      집단당     올림 후 총
  2   127.531   63.765        128
  3   157.189   52.396        159
  4   178.396   44.599        180
  5   195.766   39.153        200
  6   210.845   35.141        216

import matplotlib.pyplot as plt

# 효과크기를 세 가지로 두고 표본크기에 따른 검정력을 그린다.
# 세 곡선이 모두 위로 볼록하다는 점이 중요하다. 표본을 늘릴수록 얻는 것이
# 줄어들므로, 검정력 0.8 을 0.9 로 올리는 비용이 0.5 를 0.8 로 올리는 비용보다 크다.
fig, ax = plt.subplots(figsize=(10, 6))

analysis = TTestPower()
sample_sizes = np.arange(10, 200, 5)

for d in [0.2, 0.5, 0.8]:
    power_values = [analysis.power(effect_size=d, nobs=n, alpha=0.05)
                    for n in sample_sizes]
    ax.plot(sample_sizes, power_values, linewidth=2, label=f"d = {d:.1f}")

# 관례로 쓰는 두 기준선. 0.80 이 가장 흔하다.
ax.axhline(0.80, color='red', linestyle='--', linewidth=1, label='Power = 0.80')
ax.axhline(0.90, color='orange', linestyle='--', linewidth=1, label='Power = 0.90')

ax.set_xlabel('Sample Size (n)', fontsize=12)
ax.set_ylabel('Power', fontsize=12)
ax.set_title('Power Curves: Sample Size vs. Power\n(One-Sample t-Test, α = 0.05)')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_ylim([0, 1])

plt.tight_layout()
plt.show()

# (1) 의 z 어림과 (2) 의 정확한 답을 나란히 둔다.
za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
print(f"{'d':>5} {'z 어림':>10} {'비중심 t':>10} {'차이':>8}")
exact = {}
for d in (0.2, 0.5, 0.8):
    nz = ((za + zb) / d) ** 2
    nt = analysis.solve_power(effect_size=d, alpha=0.05, power=0.80,
                              alternative='two-sided')
    exact[d] = nt
    print(f"{d:>5.1f} {nz:>10.3f} {nt:>10.3f} {nt - nz:>+8.3f}")

print(f"\n비 n(0.2)/n(0.8):  z 어림 {16:.2f}"
      f"   정확 {exact[0.2] / exact[0.8]:.2f}")
print(f"비 n(0.2)/n(0.5):  z 어림 {6.25:.2f}"
      f"   정확 {exact[0.2] / exact[0.5]:.2f}")

# 검정력을 더 올리는 비용 (d = 0.5)
print()
for pw in (0.80, 0.90, 0.95):
    n = analysis.solve_power(effect_size=0.5, alpha=0.05, power=pw,
                             alternative='two-sided')
    print(f"  d=0.5, 검정력 {pw:.2f}: n = {n:.3f} -> {int(np.ceil(n))}")

    d       z 어림      비중심 t       차이
  0.2    196.222    198.151   +1.929
  0.5     31.396     33.367   +1.972
  0.8     12.264     14.303   +2.039

비 n(0.2)/n(0.8):  z 어림 16.00   정확 13.85
비 n(0.2)/n(0.5):  z 어림 6.25   정확 5.94

  d=0.5, 검정력 0.80: n = 33.367 -> 34
  d=0.5, 검정력 0.90: n = 43.995 -> 44
  d=0.5, 검정력 0.95: n = 53.941 -> 54

def design_study(test_type, effect_size, alpha=0.05, power=0.80,
                 **kwargs):
    """연구 설계를 위한 검정력 분석을 한자리에 모은 함수.

    검정 종류마다 statsmodels의 클래스가 다르고 인자 이름도 다르다.
    그 차이를 여기서 흡수해 두면 설계 단계에서 시나리오를 바꿔 가며
    표본크기를 비교하기가 쉬워진다.

    test_type : 'one-sample', 'two-sample', 'anova', 'proportions'
    effect_size : 표준화된 효과크기 (종류마다 척도가 다르다: d, d, f)
    alpha : 유의수준
    power : 목표 검정력
    **kwargs : 검정별 추가 인자 (예: ANOVA의 k_groups)
    """
    results = {'test_type': test_type, 'alpha': alpha, 'power': power}

    if test_type == 'one-sample':
        analysis = TTestPower()
        n = analysis.solve_power(effect_size=effect_size, alpha=alpha,
                                power=power, alternative='two-sided')
        results['sample_size'] = int(np.ceil(n))

    elif test_type == 'two-sample':
        analysis = TTestIndPower()
        n = analysis.solve_power(effect_size=effect_size, alpha=alpha,
                                power=power, ratio=1.0, alternative='two-sided')
        results['sample_size_per_group'] = int(np.ceil(n))
        results['total_sample_size'] = 2 * int(np.ceil(n))

    elif test_type == 'anova':
        analysis = FTestAnovaPower()
        k = kwargs.get('k_groups', 3)
        # **FTestAnovaPower 의 nobs 는 총 관측수다**(보기 6).
        # 집단당으로 읽으면 표본을 k 배로 부풀리게 된다.
        n_total = analysis.solve_power(effect_size=effect_size, alpha=alpha,
                                       power=power, k_groups=k)
        each = int(np.ceil(n_total / k))
        results['k_groups'] = k
        results['sample_size_per_group'] = each
        results['total_sample_size'] = k * each

    return results

# 사용 예
print("\n" + "="*60)
print("STUDY DESIGN: Two-Sample Comparison")
print("="*60)
design = design_study('two-sample', effect_size=0.5)
for key, value in design.items():
    print(f"{key:.<30} {value}")

# 세 가지를 나란히 돌려 본다.
print("\n" + "="*60)
for args in [('one-sample', 0.5, {}),
             ('two-sample', 0.5, {}),
             ('anova', 0.25, {'k_groups': 4})]:
    kind, es, kw = args
    print(f"{kind:>12}  es={es}  ->  {design_study(kind, es, **kw)}")

# 말없이 실패하는 자리.
print("\n" + "="*60)
print("proportions :", design_study('proportions', 0.0475))
print("오타 'twosample':", design_study('twosample', 0.5))

============================================================
STUDY DESIGN: Two-Sample Comparison
============================================================
test_type..................... two-sample
alpha......................... 0.05
power......................... 0.8
sample_size_per_group......... 64
total_sample_size............. 128

============================================================
  one-sample  es=0.5  ->  {'test_type': 'one-sample', 'alpha': 0.05, 'power': 0.8, 'sample_size': 34}
  two-sample  es=0.5  ->  {'test_type': 'two-sample', 'alpha': 0.05, 'power': 0.8, 'sample_size_per_group': 64, 'total_sample_size': 128}
       anova  es=0.25  ->  {'test_type': 'anova', 'alpha': 0.05, 'power': 0.8, 'k_groups': 4, 'sample_size_per_group': 45, 'total_sample_size': 180}

============================================================
proportions : {'test_type': 'proportions', 'alpha': 0.05, 'power': 0.8}
오타 'twosample': {'test_type': 'twosample', 'alpha': 0.05, 'power': 0.8}

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

def power_t2(n, d, alpha=0.05):
    """이표본 t 검정(집단당 n)의 정확한 검정력"""
    nu = 2 * n - 2
    lam = d * np.sqrt(n / 2)
    tc = stats.t.ppf(1 - alpha / 2, nu)
    return stats.nct.sf(tc, nu, lam) + stats.nct.cdf(-tc, nu, lam)

# ① n 구하기
n = next(m for m in range(4, 5000) if power_t2(m, 0.5) >= 0.80)
print(f"① d=0.5, α=0.05, 검정력 0.80  →  집단당 n = {n}")

# ② 검정력 구하기
print(f"② d=0.5, α=0.05, n=40        →  검정력 = {power_t2(40, 0.5):.4f}")

# ③ 탐지 가능한 최소 효과크기
d = brentq(lambda x: power_t2(40, x) - 0.80, 0.01, 3.0)
print(f"③ α=0.05, n=40, 검정력 0.80  →  d = {d:.4f}")

# ④ 필요한 α
a = brentq(lambda x: power_t2(40, 0.5, x) - 0.80, 1e-6, 0.5)
print(f"④ d=0.5, n=40, 검정력 0.80   →  α = {a:.4f}")

① d=0.5, α=0.05, 검정력 0.80  →  집단당 n = 64
② d=0.5, α=0.05, n=40        →  검정력 = 0.5981
③ α=0.05, n=40, 검정력 0.80  →  d = 0.6343
④ d=0.5, n=40, 검정력 0.80   →  α = 0.1672

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
n, M = 25, 10_000

def run(gen, shift):
    cnt = np.zeros(3)
    for _ in range(M):
        x = gen(n)
        y = gen(n) + shift
        cnt[0] += stats.ttest_ind(x, y, equal_var=False).pvalue < 0.05
        cnt[1] += stats.mannwhitneyu(x, y).pvalue < 0.05
        # 부호검정 (중앙값 기준, 두 표본을 짝지어 비교하는 대신 대략적 대안)
        med = np.median(np.concatenate([x, y]))
        tab = [[np.sum(x > med), np.sum(x <= med)],
               [np.sum(y > med), np.sum(y <= med)]]
        cnt[2] += stats.chi2_contingency(tab).pvalue < 0.05
    return cnt / M

cases = [("정규", lambda m: rng.normal(0, 1, m), 0.8),
         ("t(3)", lambda m: rng.standard_t(3, m), 0.8),
         ("지수", lambda m: rng.exponential(1, m), 0.8),
         ("로그정규", lambda m: rng.lognormal(0, 1, m), 1.2)]
print(f"{'분포':>8s} {'웰치 t':>9s} {'만-휘트니':>11s} {'중앙값검정':>11s}")
for name, gen, sh in cases:
    r = run(gen, sh)
    print(f"{name:>8s} {r[0]:9.4f} {r[1]:11.4f} {r[2]:11.4f}")

      분포      웰치 t       만-휘트니       중앙값검정
      정규    0.7854      0.7585      0.4788
    t(3)    0.4438      0.5716      0.3980
      지수    0.7939      0.9406      0.6787
    로그정규    0.6173      0.9597      0.7971

import numpy as np
from scipy import stats

rng = np.random.default_rng(8)
M, k, n = 40_000, 3, 64
d = 0.5

print(f"{'ρ':>5s} {'개별':>8s} {'적어도 하나':>12s} {'모두':>8s} "
      f"{'본페로니 후 하나':>16s}")
for rho in [0.0, 0.3, 0.6, 0.9]:
    S = rho * np.ones((k, k)) + (1 - rho) * np.eye(k)
    L = np.linalg.cholesky(S)
    z = rng.standard_normal((M, k)) @ L.T + d * np.sqrt(n / 2)
    p = 2 * stats.norm.sf(np.abs(z))
    sig = p < 0.05
    sig_b = p < 0.05 / k
    print(f"{rho:5.1f} {sig[:, 0].mean():8.4f} {sig.any(1).mean():12.4f} "
          f"{sig.all(1).mean():8.4f} {sig_b.any(1).mean():16.4f}")

    ρ       개별       적어도 하나       모두        본페로니 후 하나
  0.0   0.8105       0.9932   0.5231           0.9631
  0.3   0.8115       0.9725   0.5853           0.9170
  0.6   0.8065       0.9382   0.6451           0.8563
  0.9   0.8092       0.8797   0.7304           0.7667

import numpy as np
from scipy import stats

def power_t2(n, d=0.5, alpha=0.05):
    nu = max(2 * n - 2, 1)
    lam = d * np.sqrt(n / 2)
    tc = stats.t.ppf(1 - alpha / 2, nu)
    return stats.nct.sf(tc, nu, lam) + stats.nct.cdf(-tc, nu, lam)

n0 = 64                                   # 계획 표본(검정력 0.80)
print(f"{'탈락률':>8s} {'유효 n':>8s} {'실제 검정력':>12s} "
      f"{'보정 모집 n':>12s}")
for q in [0.0, 0.05, 0.10, 0.20, 0.30, 0.50]:
    neff = int(np.floor(n0 * (1 - q)))
    nadj = int(np.ceil(n0 / (1 - q)))
    print(f"{q:8.0%} {neff:8d} {power_t2(neff):12.4f} {nadj:12d}")

     탈락률     유효 n       실제 검정력      보정 모집 n
      0%       64       0.8015           64
      5%       60       0.7753           68
     10%       57       0.7538           72
     20%       51       0.7056           80
     30%       44       0.6402           92
     50%       32       0.5036          128

for q in [0.10, 0.20, 0.30]:
    print(f"개체 탈락 {q:.0%} → 쌍 손실 {1 - (1 - q)**2:.1%}, "
          f"필요 모집 배율 {1 / (1 - q)**2:.2f}배")

개체 탈락 10% → 쌍 손실 19.0%, 필요 모집 배율 1.23배
개체 탈락 20% → 쌍 손실 36.0%, 필요 모집 배율 1.56배
개체 탈락 30% → 쌍 손실 51.0%, 필요 모집 배율 2.04배