이표본 평균 검정¶

import math
from scipy.stats import t as tdist

def test_diff_two_means(n1, m1, s1, n2, m2, s2, method="welch",
                        delta0=0.0, alt="two-sided", alpha=0.05):
    """귀무가설 H0: mu1 - mu2 = delta0.

    delta0을 0이 아닌 값으로 둘 수 있게 해 두었다.
    "차이가 있는가"가 아니라 "차이가 5 이상인가"를 묻는 동등성·비열등성
    검정에서 이 자리가 쓰인다.
    method='welch'(기본) 또는 'pooled'. 돌려주는 값은 (t, df, p, 기각 여부).
    """
    diff_hat = m1 - m2
    if method == "welch":
        # 두 분산을 따로 둔 채 더한다. 합동하지 않는다.
        se = math.sqrt(s1**2 / n1 + s2**2 / n2)
        num = (s1**2 / n1 + s2**2 / n2) ** 2
        den = (s1**2 / n1)**2 / (n1 - 1) + (s2**2 / n2)**2 / (n2 - 1)
        df = num / den
    else:
        # 합동: 두 분산이 같다고 보고 자유도로 가중평균한다.
        # 이 가정이 틀리면 표준오차가 편향되고 t가 t분포를 따르지 않는다.
        df = n1 + n2 - 2
        sp2 = ((n1 - 1) * s1**2 + (n2 - 1) * s2**2) / df
        se = math.sqrt(sp2 * (1 / n1 + 1 / n2))

    t = (diff_hat - delta0) / 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, df, p, (p < alpha)

import numpy as np

rng = np.random.default_rng(0)
K = 200_000
n1 = rng.integers(2, 201, K)
n2 = rng.integers(2, 201, K)
s1 = 10.0 ** rng.uniform(-3, 3, K)      # 분산비를 10^12 까지 벌린다
s2 = 10.0 ** rng.uniform(-3, 3, K)

a, b = s1**2 / n1, s2**2 / n2
nu = (a + b) ** 2 / (a**2 / (n1 - 1) + b**2 / (n2 - 1))
lo = np.minimum(n1, n2) - 1
hi = n1 + n2 - 2
print(f"범위 검사 {K:,} 개")
print(f"  nu < min(n1,n2)-1 인 경우 = {int(np.sum(nu < lo - 1e-9))}")
print(f"  nu > n1+n2-2    인 경우 = {int(np.sum(nu > hi + 1e-9))}")
print(f"  하한에 가장 가까운 비 nu/lo 의 최소 = {np.min(nu / lo):.12f}")
print(f"  상한에 가장 가까운 비 nu/hi 의 최대 = {np.max(nu / hi):.12f}")

# 두 끝을 직접 만들어 본다. n1 = 12, n2 = 10 고정.
print("\n  s1      s2        welch nu      합동 df")
for ss1, ss2 in ((1.0, 1.0), (1.0, 1e3), (1e3, 1.0),
                 (math.sqrt(12 * 11), math.sqrt(10 * 9))):
    _, nu_w, _, _ = test_diff_two_means(12, 0, ss1, 10, 0, ss2, method="welch")
    _, nu_p, _, _ = test_diff_two_means(12, 0, ss1, 10, 0, ss2, method="pooled")
    print(f"{ss1:8.4f} {ss2:8.4f}  {nu_w:12.6f}  {nu_p:8.0f}")

# 표준오차 항등식
print("\n표준오차 항등식 검사")
worst = 0.0
for _ in range(200_000):
    m, k = int(rng.integers(2, 201)), int(rng.integers(2, 201))
    u, v = 10.0 ** rng.uniform(-3, 3), 10.0 ** rng.uniform(-3, 3)
    N = m + k
    se_w2 = u**2 / m + v**2 / k
    sp2 = ((m - 1) * u**2 + (k - 1) * v**2) / (N - 2)
    se_p2 = sp2 * (1 / m + 1 / k)
    pred = (N - 1) * (m - k) * (u**2 - v**2) / ((N - 2) * m * k)
    worst = max(worst, abs(se_p2 - se_w2 - pred) / max(se_w2, 1e-300))
print(f"  상대오차의 최대값 = {worst:.3e}")

범위 검사 200,000 개
  nu < min(n1,n2)-1 인 경우 = 0
  nu > n1+n2-2    인 경우 = 0
  하한에 가장 가까운 비 nu/lo 의 최소 = 1.000000000000
  상한에 가장 가까운 비 nu/hi 의 최대 = 0.999999999798

  s1      s2        welch nu      합동 df
  1.0000   1.0000     19.289855        20
  1.0000 1000.0000      9.000015        20
1000.0000   1.0000     11.000026        20
 11.4891   9.4868     20.000000        20

표준오차 항등식 검사
  상대오차의 최대값 = 3.773e-14

t, df, p, reject = test_diff_two_means(
    n1=12, m1=0.0, s1=1.0, n2=10, m2=0.5, s2=1.5,
    method="welch", alt="greater"
)
print("t:", t, "df:", df, "p:", p, "reject:", reject)

# 같은 자료를 합동 t로도 해 본다. 자유도가 어떻게 달라지는지 보라.
t_p, df_p, p_p, reject_p = test_diff_two_means(
    n1=12, m1=0.0, s1=1.0, n2=10, m2=0.5, s2=1.5,
    method="pooled", alt="greater"
)
print("t:", t_p, "df:", df_p, "p:", p_p, "reject:", reject_p)

t: -0.9004503377814964 df: 15.195761856710394 p: 0.8090351110315042 reject: False
t: -0.9341987329938274 df: 20 p: 0.8193278973550725 reject: False

import numpy as np
from scipy import stats

n1, n2, sg1, sg2 = 12, 10, 1.0, 1.5
d1, d2, N = n1 - 1, n2 - 1, n1 + n2
c = 1 / n1 + 1 / n2

# 두 표준오차와 항등식
se_w = math.sqrt(sg1**2 / n1 + sg2**2 / n2)
sp2 = (d1 * sg1**2 + d2 * sg2**2) / (N - 2)
se_p = math.sqrt(sp2 * c)
ident = (N - 1) * (n1 - n2) * (sg1**2 - sg2**2) / ((N - 2) * n1 * n2)
print(f"SE_welch = {se_w:.9f}   SE_pool = {se_p:.9f}")
print(f"SE_p^2 - SE_w^2 = {se_p**2 - se_w**2:.12f}   항등식 = {ident:.12f}")
print(f"t_welch = {-0.5 / se_w:.9f}   t_pool = {-0.5 / se_p:.9f}")

# 5.3 절의 두 손잡이 R 과 m
tau2 = sg1**2 / n1 + sg2**2 / n2
EVp = c * (d1 * sg1**2 + d2 * sg2**2) / (d1 + d2)
R = EVp / tau2
m = (d1 * sg1**2 + d2 * sg2**2) ** 2 / (d1 * sg1**4 + d2 * sg2**4)
crit = stats.t.ppf(0.975, N - 2)
print(f"\nR = {R:.6f},  m = {m:.6f},  t_(0.975, 20) = {crit:.6f}")
print(f"다듬은 어림 P(|t_m| > crit*sqrt(R)) = "
      f"{2 * stats.t.sf(crit * math.sqrt(R), m):.6f}")
print(f"거친 어림   P(|Z|   > crit*sqrt(R)) = "
      f"{2 * stats.norm.sf(crit * math.sqrt(R)):.6f}")

# 모의실험
rng = np.random.default_rng(7)
B = 200_000
x1 = rng.normal(0, sg1, (B, n1))
x2 = rng.normal(0, sg2, (B, n2))
dd = x1.mean(1) - x2.mean(1)
v1, v2 = x1.var(1, ddof=1), x2.var(1, ddof=1)

t_pool = dd / np.sqrt(((d1 * v1 + d2 * v2) / (N - 2)) * c)
a, b = v1 / n1, v2 / n2
t_welch = dd / np.sqrt(a + b)
nu = (a + b) ** 2 / (a**2 / d1 + b**2 / d2)
print(f"\n반복 {B:,} 회   (오류율 0.05 둘레의 몬테카를로 오차 ±0.0005)")
print(f"  합동 t 의 실제 오류율  = {np.mean(np.abs(t_pool) > crit):.5f}")
print(f"  Welch 의 실제 오류율   = "
      f"{np.mean(np.abs(t_welch) > stats.t(nu).ppf(0.975)):.5f}")
print(f"  실현된 nu: 평균 {nu.mean():.4f}, 최소 {nu.min():.4f}, 최대 {nu.max():.4f}"
      f"   (한계 {min(n1, n2) - 1} 과 {N - 2})")

SE_welch = 0.555277708   SE_pool = 0.535218024
SE_p^2 - SE_w^2 = -0.021875000000   항등식 = -0.021875000000
t_welch = -0.900450338   t_pool = -0.934198733

R = 0.929054,  m = 17.265193,  t_(0.975, 20) = 2.085963
다듬은 어림 P(|t_m| > crit*sqrt(R)) = 0.060252
거친 어림   P(|Z|   > crit*sqrt(R)) = 0.044367

반복 200,000 회   (오류율 0.05 둘레의 몬테카를로 오차 ±0.0005)
  합동 t 의 실제 오류율  = 0.05967
  Welch 의 실제 오류율   = 0.04979
  실현된 nu: 평균 15.4196, 최소 9.2603, 최대 20.0000   (한계 9 과 20)

import numpy as np
from scipy import stats

def welch_df(s1, n1, s2, n2):
    a, b = s1**2 / n1, s2**2 / n2
    return (a + b)**2 / (a**2 / (n1 - 1) + b**2 / (n2 - 1))

x1, s1, n1 = 82.0, 6.0, 15
x2, s2, n2 = 76.0, 8.0, 12
diff = x1 - x2

se_w = np.sqrt(s1**2 / n1 + s2**2 / n2)
df_w = welch_df(s1, n1, s2, n2)
t_w = diff / se_w
tc_w = stats.t.ppf(0.975, df_w)
print(f"Welch  t={t_w:.4f}  df={df_w:.3f}  "
      f"p={2 * stats.t.sf(abs(t_w), df_w):.4f}")
print(f"       95% CI ({diff - tc_w * se_w:.3f}, {diff + tc_w * se_w:.3f})")

sp = np.sqrt(((n1 - 1) * s1**2 + (n2 - 1) * s2**2) / (n1 + n2 - 2))
se_p, df_p = sp * np.sqrt(1 / n1 + 1 / n2), n1 + n2 - 2
t_p = diff / se_p
tc_p = stats.t.ppf(0.975, df_p)
print(f"합동   t={t_p:.4f}  df={df_p}  "
      f"p={2 * stats.t.sf(abs(t_p), df_p):.4f}")
print(f"       95% CI ({diff - tc_p * se_p:.3f}, {diff + tc_p * se_p:.3f})")

print(f"\n자유도가 움직이는 범위: [{min(n1, n2) - 1}, {n1 + n2 - 2}]")
print(f"{'s1':>5s} {'s2':>5s} {'Welch df':>10s}")
for a, b in [(1, 1), (1, 2), (1, 5), (1, 20)]:
    print(f"{a:5d} {b:5d} {welch_df(a, n1, b, n2):10.3f}")

Welch  t=2.1576  df=19.953  p=0.0433
       95% CI (0.198, 11.802)
합동   t=2.2287  df=25  p=0.0351
       95% CI (0.455, 11.545)

자유도가 움직이는 범위: [11, 25]
   s1    s2   Welch df
    1     1     23.715
    1     2     15.357
    1     5     11.706
    1    20     11.044

import numpy as np
from scipy import stats

N = 100
print(f"{'σ1:σ2':>8s} {'균등 SE':>10s} {'최적 배분':>14s} "
      f"{'최적 SE':>10s} {'분산 효율':>10s}")
for s1, s2 in [(1, 1), (1, 2), (1, 3), (1, 5)]:
    se_eq = np.sqrt(s1**2 / (N / 2) + s2**2 / (N / 2))
    n1 = N * s1 / (s1 + s2)
    se_op = np.sqrt(s1**2 / n1 + s2**2 / (N - n1))
    print(f"{s1}:{s2:<6d} {se_eq:10.4f} {n1:6.1f}:{N - n1:<7.1f} "
          f"{se_op:10.4f} {(se_eq / se_op)**2:10.4f}")

def power(n1, n2, delta, s1, s2, alpha=0.05):
    a, b = s1**2 / n1, s2**2 / n2
    se, df = np.sqrt(a + b), (a + b)**2 / (a**2 / (n1 - 1) + b**2 / (n2 - 1))
    c = stats.t.ppf(1 - alpha / 2, df)
    return stats.nct.sf(c, df, delta / se) + stats.nct.cdf(-c, df, delta / se)

print(f"\nσ1=1, σ2=3, 참 차이 1, 총 {N}명일 때")
for n1 in [25, 33, 50, 60, 75]:
    print(f"  n1={n1:3d}, n2={N - n1:3d}  →  검정력 "
          f"{power(n1, N - n1, 1.0, 1.0, 3.0):.4f}")

   σ1:σ2      균등 SE          최적 배분      최적 SE      분산 효율
1:1          0.2000   50.0:50.0        0.2000     1.0000
1:2          0.3162   33.3:66.7        0.3000     1.1111
1:3          0.4472   25.0:75.0        0.4000     1.2500
1:5          0.7211   16.7:83.3        0.6000     1.4444

σ1=1, σ2=3, 참 차이 1, 총 100명일 때
  n1= 25, n2= 75  →  검정력 0.6969
  n1= 33, n2= 67  →  검정력 0.6837
  n1= 50, n2= 50  →  검정력 0.5949
  n1= 60, n2= 40  →  검정력 0.5122
  n1= 75, n2= 25  →  검정력 0.3506

import numpy as np
from scipy import stats

def yuen(x, y, tr=0.2):
    """절사비율 tr 의 유엔 검정. (t, p, df) 를 돌려준다."""
    def wins_var(a, g):
        s = np.sort(a)
        n = len(a)
        s[:g] = s[g]                 # 아래 g 개를 g번째 값으로 대치
        s[n - g:] = s[n - g - 1]     # 위 g 개도 마찬가지
        return s.var(ddof=1)

    g1, g2 = int(np.floor(tr * len(x))), int(np.floor(tr * len(y)))
    h1, h2 = len(x) - 2 * g1, len(y) - 2 * g2
    d1 = (len(x) - 1) * wins_var(x, g1) / (h1 * (h1 - 1))
    d2 = (len(y) - 1) * wins_var(y, g2) / (h2 * (h2 - 1))
    t = (stats.trim_mean(x, tr) - stats.trim_mean(y, tr)) / np.sqrt(d1 + d2)
    df = (d1 + d2)**2 / (d1**2 / (h1 - 1) + d2**2 / (h2 - 1))
    return t, 2 * stats.t.sf(abs(t), df), df

rng = np.random.default_rng(626)
M, n = 5_000, 25
print(f"{'상황':>10s} {'Welch 수준':>11s} {'Yuen 수준':>11s} "
      f"{'Welch 검정력':>13s} {'Yuen 검정력':>13s}")
for name in ["정규", "t(3)", "오염 10%", "로그정규"]:
    res = []
    for shift in [0.0, 0.8]:
        a = b = 0
        for _ in range(M):
            if name == "정규":
                x, y = rng.standard_normal(n), rng.standard_normal(n) + shift
            elif name == "t(3)":
                x, y = rng.standard_t(3, n), rng.standard_t(3, n) + shift
            elif name == "오염 10%":
                # 10% 의 관측값이 표준편차 5배로 흩어진다
                x = rng.standard_normal(n) * np.where(rng.random(n) < 0.1, 5, 1)
                y = rng.standard_normal(n) * np.where(rng.random(n) < 0.1, 5, 1) + shift
            else:
                x, y = rng.lognormal(0, 1, n), rng.lognormal(0, 1, n) + shift
            a += stats.ttest_ind(x, y, equal_var=False).pvalue < 0.05
            b += yuen(x, y)[1] < 0.05
        res += [a / M, b / M]
    print(f"{name:>10s} {res[0]:11.4f} {res[1]:11.4f} "
          f"{res[2]:13.4f} {res[3]:13.4f}")

        상황    Welch 수준     Yuen 수준     Welch 검정력      Yuen 검정력
        정규      0.0528      0.0466        0.7882        0.7148
      t(3)      0.0432      0.0474        0.4344        0.5790
    오염 10%      0.0432      0.0482        0.3790        0.6166
      로그정규      0.0378      0.0388        0.3668        0.6068

import numpy as np
from scipy import stats

x1, s1, n1 = 82.0, 6.0, 15
x2, s2, n2 = 76.0, 8.0, 12
diff = x1 - x2
a, b = s1**2 / n1, s2**2 / n2
se = np.sqrt(a + b)
df = (a + b)**2 / (a**2 / (n1 - 1) + b**2 / (n2 - 1))

print(f"평균 차이 {diff:.1f},  SE {se:.4f},  df {df:.2f}\n")
for D in [5.0, 10.0, 15.0]:
    p_lo = stats.t.sf((diff + D) / se, df)      # H01: μ1-μ2 ≤ -Δ
    p_hi = stats.t.cdf((diff - D) / se, df)     # H02: μ1-μ2 ≥ +Δ
    tost = max(p_lo, p_hi)
    print(f"Δ={D:5.1f}  TOST p={tost:.4f}  → "
          f"{'동등 인정' if tost < 0.05 else '동등 미인정'}")

tc90 = stats.t.ppf(0.95, df)
print(f"\n90% CI ({diff - tc90 * se:.3f}, {diff + tc90 * se:.3f})")

D = 5.0
t_ni = (diff + D) / se
print(f"비열등성(Δ=5, '집단 1이 5점 이상 나쁘지는 않다')  "
      f"t={t_ni:.4f}, 단측 p={stats.t.sf(t_ni, df):.6f}")

평균 차이 6.0,  SE 2.7809,  df 19.95

Δ=  5.0  TOST p=0.6385  → 동등 미인정
Δ= 10.0  TOST p=0.0829  → 동등 미인정
Δ= 15.0  TOST p=0.0021  → 동등 인정

90% CI (1.203, 10.797)
비열등성(Δ=5, '집단 1이 5점 이상 나쁘지는 않다')  t=3.9556, 단측 p=0.000392

두 집단의 중심을 비교하고 싶다
    │
    ├─ 자료가 짝지어져 있는가 ──── 예 ──→ 대응 검정 (앞 절)
    │
    └─ 아니오 (독립 두 집단)
         │
         ├─ 1단계: 그림을 그린다
         │     상자그림 · 정규분위수그림 · 산점도
         │     → 이상점, 치우침, 분산 차이, 다봉성을 눈으로 확인
         │
         ├─ 2단계: 무엇을 비교할 것인지 정한다
         │     ├─ 평균 (총액·합계가 관심) ──→ Welch t
         │     ├─ 절사평균 (이상점 있음) ──→ 유엔의 검정
         │     ├─ 순위 (전형적 값) ────────→ 만·휘트니
         │     └─ 분포 전체 ──────────────→ 콜모고로프·스미르노프
         │
         ├─ 3단계: 검정 + 신뢰구간 + 효과크기
         │
         └─ 4단계: 동등성이 관심이면 TOST 추가