대응 평균 검정¶

import math
from scipy.stats import t as tdist

def test_paired_mean(n, dbar, sd_d, mu_d0=0.0,
                     alt="two-sided", alpha=0.05):
    """대응 t-검정. 차이 D = X - Y에 대한 일표본 검정과 같다.

    받는 것은 짝의 개수 n, 차이의 평균, 차이의 표준편차뿐이다.
    두 집단의 산포나 상관을 따로 알 필요가 없다. 짝을 지으며 이미 흡수했기 때문이다.
    """
    df = n - 1               # 짝의 개수 - 1. 관측값 2n개가 아니다.
    se = sd_d / math.sqrt(n)
    t = (dbar - mu_d0) / se
    if alt == "two-sided":
        p = 2 * min(tdist.cdf(t, df), 1 - tdist.cdf(t, df))
    elif alt == "less":
        p = tdist.cdf(t, df)
    else:
        p = 1 - tdist.cdf(t, df)
    return t, p, (p < alpha)

from scipy import stats
from scipy.optimize import brentq

print("  n   t_{n-1}   t_{2n-2}    rho*")
for n in (5, 10, 30):
    t1 = stats.t.ppf(0.975, n - 1)
    t2 = stats.t.ppf(0.975, 2 * n - 2)
    print(f"{n:4d}   {t1:.4f}    {t2:.4f}    {1 - (t2 / t1)**2:.4f}")


def pow_paired(n, rho, mu):
    """쌍 n 개, 상관 rho, 참 평균차 mu 일 때 대응 검정의 검정력."""
    nc = mu * math.sqrt(n) / math.sqrt(2 * (1 - rho))
    tc = stats.t.ppf(0.975, n - 1)
    return stats.nct.sf(tc, n - 1, nc) + stats.nct.cdf(-tc, n - 1, nc)


def pow_two(n, mu):
    """각 n 명인 독립 두 집단의 검정력. 자유도가 2n-2 로 늘어난다."""
    nc = mu * math.sqrt(n) / math.sqrt(2)
    tc = stats.t.ppf(0.975, 2 * n - 2)
    return stats.nct.sf(tc, 2 * n - 2, nc) + stats.nct.cdf(-tc, 2 * n - 2, nc)


print("\n검정력을 같게 만드는 rho (sigma = 1, 참 mu_D 를 바꾸며)")
print("   n    mu_D   독립 검정력    rho*")
for n in (5, 10, 30):
    for mu in (0.5, 1.0, 1.5):
        r = brentq(lambda rr: pow_paired(n, rr, mu) - pow_two(n, mu), -0.5, 0.95)
        print(f"{n:4d}   {mu:.1f}      {pow_two(n, mu):.4f}      {r:.4f}")

  n   t_{n-1}   t_{2n-2}    rho*
   5   2.7764    2.3060    0.3102
  10   2.2622    2.1009    0.1375
  30   2.0452    2.0017    0.0421

검정력을 같게 만드는 rho (sigma = 1, 참 mu_D 를 바꾸며)
   n    mu_D   독립 검정력    rho*
   5   0.5      0.1077      0.2191
   5   1.0      0.2863      0.2368
   5   1.5      0.5494      0.2563
  10   0.5      0.1851      0.1052
  10   1.0      0.5620      0.1112
  10   1.5      0.8870      0.1173
  30   0.5      0.4779      0.0335
  30   1.0      0.9677      0.0346
  30   1.5      0.9999      0.0356

import numpy as np

dbar, sd_d = 0.4, 1.1

t_stat, p_less, rej_less = test_paired_mean(
    n=12, dbar=dbar, sd_d=sd_d, mu_d0=0.0, alt="less"
)
t2, p_gre, rej_gre = test_paired_mean(
    n=12, dbar=dbar, sd_d=sd_d, mu_d0=0.0, alt="greater"
)
_, p_two, rej_two = test_paired_mean(
    n=12, dbar=dbar, sd_d=sd_d, mu_d0=0.0, alt="two-sided"
)
print(f"t = {t_stat:.12f}   (두 호출이 같은가: {t_stat == t2})")
print(f"  alt='less'      p = {p_less:.12f}  기각 {rej_less}")
print(f"  alt='greater'   p = {p_gre:.12f}  기각 {rej_gre}")
print(f"  두 p 의 합       = {p_less + p_gre:.12f}")
print(f"  alt='two-sided' p = {p_two:.12f}  기각 {rej_two}"
      f"   2*min = {2 * min(p_less, p_gre):.12f}")

# (2) dbar 와 sd_d 를 고정하고 쌍의 개수만 늘린다.
d_over_s = dbar / sd_d
print(f"\nd/s = {d_over_s:.6f}")
for q, name in ((0.95, "단측"), (0.975, "양측")):
    crit = stats.norm.ppf(q)
    n_norm = (crit / d_over_s) ** 2
    n_fix = brentq(lambda n: d_over_s * math.sqrt(n) - stats.t.ppf(q, n - 1),
                   2.5, 5000)
    print(f"  {name}: 정규근사 n > {n_norm:.4f} → {math.ceil(n_norm):d},"
          f"   t 고정점 n > {n_fix:.4f} → {math.ceil(n_fix):d}")

# 전수 탐색. "조건을 만족하는 첫 칸" 만 보지 않는다.
ns = np.arange(2, 2001)
tt = d_over_s * np.sqrt(ns)
p_one = stats.t.sf(tt, ns - 1)
p_dbl = 2 * p_one
for p, name in ((p_one, "단측"), (p_dbl, "양측")):
    ok = ns[p < 0.05]
    print(f"\n{name}: p<0.05 인 n 의 개수 = {len(ok)},  가장 작은 n = {ok.min()},"
          f"  연속인가 = {bool(np.all(np.diff(ok) == 1))}")
    print(f"      p 가 n 에 대해 단조감소인가 = {bool(np.all(np.diff(p) < 0))}"
          f"   (최대 p = {p.max():.6f} at n = {ns[p.argmax()]})")

print("\n   n        t      단측 p     양측 p")
for n in (12, 21, 22, 23, 30, 31, 32):
    t = d_over_s * math.sqrt(n)
    print(f"{n:4d}  {t:.6f}  {stats.t.sf(t, n - 1):.6f}"
          f"  {2 * stats.t.sf(t, n - 1):.6f}")

t = 1.259673314596   (두 호출이 같은가: True)
  alt='less'      p = 0.883070909978  기각 False
  alt='greater'   p = 0.116929090022  기각 False
  두 p 의 합       = 1.000000000000
  alt='two-sided' p = 0.233858180045  기각 False   2*min = 0.233858180045

d/s = 0.363636
  단측: 정규근사 n > 20.4607 → 21,   t 고정점 n > 22.3577 → 23
  양측: 정규근사 n > 29.0510 → 30,   t 고정점 n > 31.4990 → 32

단측: p<0.05 인 n 의 개수 = 1978,  가장 작은 n = 23,  연속인가 = True
      p 가 n 에 대해 단조감소인가 = True   (최대 p = 0.348806 at n = 2)

양측: p<0.05 인 n 의 개수 = 1969,  가장 작은 n = 32,  연속인가 = True
      p 가 n 에 대해 단조감소인가 = True   (최대 p = 0.697612 at n = 2)

   n        t      단측 p     양측 p
  12  1.259673  0.116929  0.233858
  21  1.666391  0.055611  0.111221
  22  1.705606  0.051415  0.102831
  23  1.743939  0.047564  0.095128
  30  1.991718  0.027947  0.055894
  31  2.024642  0.025943  0.051885
  32  2.057038  0.024090  0.048180

import math
import numpy as np
from scipy.stats import t as tdist

def test_paired_mean(n, dbar, sd_d, mu_d0=0.0,
                     alt="two-sided", alpha=0.05):
    df = n - 1
    se = sd_d / math.sqrt(n)
    t = (dbar - mu_d0) / se
    if alt == "two-sided":
        p = 2 * min(tdist.cdf(t, df), 1 - tdist.cdf(t, df))
    elif alt == "less":
        p = tdist.cdf(t, df)
    else:
        p = 1 - tdist.cdf(t, df)
    return t, p, bool(p < alpha)

d = np.array([-7, -8, -8, -7, -7, -8, -4, -7, -6, -7], float)
m, s = d.mean(), d.std(ddof=1)
print(f"차이(투여 후 - 투여 전)  평균 {m:.4f}  표준편차 {s:.4f}\n")

for alt in ["two-sided", "less", "greater"]:
    t, p, rej = test_paired_mean(10, m, s, 0, alt)
    print(f"{alt:10s} t={t:8.4f}  p={p:12.6g}  기각={rej}")

print("\n부호를 반대로 (투여 전 - 투여 후) 잡으면:")
for alt in ["less", "greater"]:
    t, p, rej = test_paired_mean(10, -m, s, 0, alt)
    print(f"{alt:10s} t={t:8.4f}  p={p:12.6g}  기각={rej}")

차이(투여 후 - 투여 전)  평균 -6.9000  표준편차 1.1972

two-sided  t=-18.2253  p= 2.05724e-08  기각=True
less       t=-18.2253  p= 1.02862e-08  기각=True
greater    t=-18.2253  p=           1  기각=False

부호를 반대로 (투여 전 - 투여 후) 잡으면:
less       t= 18.2253  p=           1  기각=False
greater    t= 18.2253  p= 1.02862e-08  기각=True

import numpy as np
from scipy import stats

def power_paired(n, delta, sd_d, alpha=0.05):
    """쌍 n개, 참 평균 차이 delta 일 때 양측 대응 t 검정의 검정력."""
    nc = delta / (sd_d / np.sqrt(n))     # 비중심 모수
    df = n - 1
    crit = stats.t.ppf(1 - alpha / 2, df)
    return stats.nct.sf(crit, df, nc) + stats.nct.cdf(-crit, df, nc)

def n_paired(delta, sigma, rho, power=0.80, alpha=0.05):
    sd_d = sigma * np.sqrt(2 * (1 - rho))
    for n in range(5, 100_000):          # n=5 부터: 극소표본은 수치가 불안정
        if power_paired(n, delta, sd_d, alpha) >= power:
            return n, sd_d, power_paired(n, delta, sd_d, alpha)

print(f"{'ρ':>5s} {'SD_D':>8s} {'필요 쌍':>8s} {'측정 횟수':>10s} {'검정력':>8s}")
for rho in [0.0, 0.3, 0.5, 0.7, 0.9]:
    n, sd, pw = n_paired(3.0, 10.0, rho)
    print(f"{rho:5.1f} {sd:8.3f} {n:8d} {2 * n:10d} {pw:8.4f}")

def n_independent(delta, sigma, power=0.80, alpha=0.05):
    """군당 표본크기."""
    for n in range(5, 100_000):
        nc = delta / (sigma * np.sqrt(2 / n))
        df = 2 * n - 2
        crit = stats.t.ppf(1 - alpha / 2, df)
        if stats.nct.sf(crit, df, nc) + stats.nct.cdf(-crit, df, nc) >= power:
            return n

n_ind = n_independent(3.0, 10.0)
print(f"\n독립 이표본: 군당 {n_ind}명, 총 {2 * n_ind}명 (측정 {2 * n_ind}회)")

    ρ     SD_D     필요 쌍      측정 횟수      검정력
  0.0   14.142      177        354   0.8015
  0.3   11.832      125        250   0.8031
  0.5   10.000       90        180   0.8038
  0.7    7.746       55        110   0.8054
  0.9    4.472       20         40   0.8121

독립 이표본: 군당 176명, 총 352명 (측정 352회)

n, _, _ = n_paired(3.0, 10.0, 0.7)
sd_actual = 10.0 * np.sqrt(2 * (1 - 0.5))
print(f"ρ=0.7 가정으로 {n}쌍 모집 → 실제 ρ=0.5 였다면 "
      f"검정력 {power_paired(n, 3.0, sd_actual):.4f}")

ρ=0.7 가정으로 55쌍 모집 → 실제 ρ=0.5 였다면 검정력 0.5891

import numpy as np
from scipy import stats

d = np.array([-1.2, 0.8, -0.3, 1.5, 0.2, -0.9, 0.6, 1.1, -0.5, 0.4,
              0.9, -0.2, 1.3, 0.1, -0.7, 0.5, 0.3, -1.0, 0.7, 0.0])
n, delta_eq = len(d), 1.0          # 동등성 한계 ±1.0
m, s = d.mean(), d.std(ddof=1)
se, df = s / np.sqrt(n), n - 1

t_lo = (m - (-delta_eq)) / se      # H01: μ ≤ -Δ 에 대한 단측
p_lo = stats.t.sf(t_lo, df)
t_hi = (m - delta_eq) / se         # H02: μ ≥ +Δ 에 대한 단측
p_hi = stats.t.cdf(t_hi, df)

print(f"n={n}  평균 차이 {m:.4f}  SD {s:.4f}  SE {se:.4f}\n")
print(f"하한 검정  t={t_lo:7.4f}  p={p_lo:.6f}")
print(f"상한 검정  t={t_hi:7.4f}  p={p_hi:.6f}")
print(f"TOST p = max = {max(p_lo, p_hi):.6f}  → "
      f"동등성 {'인정' if max(p_lo, p_hi) < 0.05 else '미인정'}\n")

tc90 = stats.t.ppf(0.95, df)       # TOST 는 90% 구간과 짝을 이룬다
tc95 = stats.t.ppf(0.975, df)
print(f"90% CI ({m - tc90 * se:7.4f}, {m + tc90 * se:7.4f})  ← ±1.0 안에 들어감")
print(f"95% CI ({m - tc95 * se:7.4f}, {m + tc95 * se:7.4f})")
print(f"보통의 대응 t 검정 p = {stats.ttest_1samp(d, 0).pvalue:.6f}")

n=20  평균 차이 0.1800  SD 0.7770  SE 0.1738

하한 검정  t= 6.7913  p=0.000001
상한 검정  t=-4.7194  p=0.000075
TOST p = max = 0.000075  → 동등성 인정

90% CI (-0.1204,  0.4804)  ← ±1.0 안에 들어감
95% CI (-0.1837,  0.5437)
보통의 대응 t 검정 p = 0.313225

import numpy as np
from scipy import stats

def power_paired(n, delta, sd_d, alpha=0.05):
    nc = delta / (sd_d / np.sqrt(n))
    df = n - 1
    crit = stats.t.ppf(1 - alpha / 2, df)
    return stats.nct.sf(crit, df, nc) + stats.nct.cdf(-crit, df, nc)

def n_for(delta, sd_d, power=0.80, alpha=0.05):
    for n in range(5, 100_000):
        if power_paired(n, delta, sd_d, alpha) >= power:
            return n
    return 100_000

delta, sd_true = 3.0, 8.0
print(f"참 SD_D={sd_true} 를 안다면 필요한 쌍 = {n_for(delta, sd_true)}\n")

rng = np.random.default_rng(7)
print(f"{'예비 m':>7s} {'계획 n 중앙값':>13s} {'5~95%':>12s} "
      f"{'실제 검정력':>12s} {'0.8 미달':>9s}")
for m in [10, 20, 40]:
    ns, pws = [], []
    for _ in range(2_000):
        # 예비연구의 S_D 를 카이제곱 분포에서 생성
        s = sd_true * np.sqrt(stats.chi2.rvs(m - 1, random_state=rng) / (m - 1))
        n = n_for(delta, s)
        ns.append(n)
        pws.append(power_paired(n, delta, sd_true))   # 실제 SD 로 평가
    ns, pws = np.array(ns), np.array(pws)
    print(f"{m:7d} {np.median(ns):13.0f} "
          f"{np.percentile(ns, 5):5.0f}~{np.percentile(ns, 95):<6.0f} "
          f"{pws.mean():12.3f} {np.mean(pws < 0.8):9.3f}")

참 SD_D=8.0 를 안다면 필요한 쌍 = 58

   예비 m      계획 n 중앙값        5~95%       실제 검정력    0.8 미달
     10            52    22~107           0.731     0.590
     20            57    32~91            0.770     0.522
     40            58    39~81            0.790     0.497

import numpy as np
from scipy import stats

d = np.array([-7, -8, -8, -7, -7, -8, -4, -7, -6, -7], float)
n = len(d)
m, s = d.mean(), d.std(ddof=1)
se = s / np.sqrt(n)
t, p = stats.ttest_1samp(d, 0)
tc = stats.t.ppf(0.975, n - 1)
print(f"n={n}  평균 차이 {m:.2f}  SD {s:.2f}  SE {se:.3f}")
print(f"95% CI ({m - tc * se:.2f}, {m + tc * se:.2f})")
print(f"t({n - 1}) = {t:.2f},  p = {p:.3g}")
print(f"d_z = {m / s:.2f}")

n=10  평균 차이 -6.90  SD 1.20  SE 0.379
95% CI (-7.76, -6.04)
t(9) = -18.23,  p = 2.06e-08
d_z = -5.76