이표본 t-검정 (합동과 Welch)¶

import numpy as np
from scipy import stats

# 체류시간(초)
page_a = np.array([185, 188, 142, 160, 161, 157, 182, 181, 159, 167])
page_b = np.array([173, 181, 182, 170, 169, 177, 168, 183, 169, 164])

# scipy의 기본값은 equal_var=True(합동)이다. Welch를 쓰려면 반드시 명시해야 한다.
t_stat, p_value = stats.ttest_ind(page_a, page_b, equal_var=False)

print(f"Page A: mean = {page_a.mean():.2f}, std = {page_a.std(ddof=1):.2f}")
print(f"Page B: mean = {page_b.mean():.2f}, std = {page_b.std(ddof=1):.2f}")
print(f"t-statistic: {t_stat:.4f}")
print(f"p-value (two-sided): {p_value:.4f}")

# 단측검정. 대립가설은 B 쪽 평균이 더 크다는 것이다.
p_one_sided = p_value / 2 if page_b.mean() > page_a.mean() else 1 - p_value / 2
print(f"p-value (one-sided): {p_one_sided:.4f}")

# 효과크기 Cohen 의 d — 평균 차이를 표준편차 단위로 잰 값
pooled_std = np.sqrt(((len(page_a) - 1) * page_a.std(ddof=1)**2 +
                       (len(page_b) - 1) * page_b.std(ddof=1)**2) /
                      (len(page_a) + len(page_b) - 2))
# 효과크기는 표본크기와 무관하다. t는 n이 커지면 함께 커지지만 d는 그렇지 않다.
# 그래서 "유의한가"와 "쓸모 있을 만큼 큰가"를 따로 말할 수 있다.
cohens_d = (page_b.mean() - page_a.mean()) / pooled_std
print(f"Cohen's d: {cohens_d:.3f}")

Page A: mean = 168.20, std = 15.08
Page B: mean = 173.60, std = 6.70
t-statistic: -1.0350
p-value (two-sided): 0.3204
p-value (one-sided): 0.1602
Cohen's d: 0.463

d = cohens_d
print(f"정확한 d = {d:.9f}   (s_p = {pooled_std:.9f})")

za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
n_norm = 2 * (za + zb) ** 2 / d**2
print(f"정규근사가 주는 집단당 n = {n_norm:.4f}  →  {int(np.ceil(n_norm))}")


def power_two_sample(n, d, alpha=0.05):
    """집단당 n 명, 참 효과크기 d 인 균형 이표본 t 검정의 검정력."""
    df = 2 * n - 2
    nc = d * np.sqrt(n / 2)            # 비중심모수
    tc = stats.t.ppf(1 - alpha / 2, df)
    return stats.nct.sf(tc, df, nc) + stats.nct.cdf(-tc, df, nc)


print("\n   n(집단당)   검정력")
for n in (10, 70, 73, 74, 75, 76, 80):
    print(f"{n:9d}   {power_two_sample(n, d):.5f}")

# 격자의 "첫 칸" 만 보지 않는다. 3 부터 400 까지 모두 훑는다.
# (n >= 757 부터는 검정력이 배정밀도로 정확히 1 이 되어 차가 0 이 된다.)
ns = np.arange(3, 401)
pw = np.array([power_two_sample(int(k), d) for k in ns])
ok = ns[pw >= 0.80]
print(f"\n검정력 >= 0.80 인 n 의 개수 = {len(ok)},  가장 작은 n = {ok.min()},"
      f"  연속인가 = {bool(np.all(np.diff(ok) == 1))}")
print(f"검정력이 n 에 대해 단조증가인가 = {bool(np.all(np.diff(pw) > 0))}")
print(f"현재 설계(n = 10)의 검정력 = {power_two_sample(10, d):.4f}")

정확한 d = 0.462857143   (s_p = 11.666666667)
정규근사가 주는 집단당 n = 73.2730  →  74

   n(집단당)   검정력
       10   0.16531
       70   0.77614
       73   0.79324
       74   0.79868
       75   0.80400
       76   0.80920
       80   0.82883

검정력 >= 0.80 인 n 의 개수 = 326,  가장 작은 n = 75,  연속인가 = True
검정력이 n 에 대해 단조증가인가 = True
현재 설계(n = 10)의 검정력 = 0.1653

from scipy import stats

group1, group2 = page_a, page_b          # 위 보기의 자료를 그대로 쓴다

# Welch t-검정 (권장). 등분산을 가정하지 않는다.
t_welch, p_welch = stats.ttest_ind(group1, group2, equal_var=False)

# 합동 t-검정. 등분산을 가정한다.
t_pooled, p_pooled = stats.ttest_ind(group1, group2, equal_var=True)

print(f"Welch : t = {t_welch:.4f}, p = {p_welch:.4f}")
print(f"Pooled: t = {t_pooled:.4f}, p = {p_pooled:.4f}")

# 단측 p-값. 통계량의 부호에 따라 처리가 달라진다.
if t_welch > 0:
    p_one_sided = p_welch / 2
else:
    p_one_sided = 1 - p_welch / 2
print(f"one-sided (H1: mu1 > mu2): p = {p_one_sided:.4f}")

Welch : t = -1.0350, p = 0.3204
Pooled: t = -1.0350, p = 0.3144
one-sided (H1: mu1 > mu2): p = 0.8398

n = len(group1)
v1, v2 = group1.var(ddof=1), group2.var(ddof=1)
a, b = v1 / n, v2 / n

se_welch = np.sqrt(a + b)
sp2 = ((n - 1) * v1 + (n - 1) * v2) / (2 * n - 2)
se_pool = np.sqrt(sp2 * (1 / n + 1 / n))
print(f"n1 = n2 = {n}")
print(f"  SE_welch = {se_welch:.12f}")
print(f"  SE_pool  = {se_pool:.12f}")
print(f"  차 = {se_pool - se_welch:.3e}   (0 이어야 한다)")
print(f"  t_welch = {t_welch:.15f}")
print(f"  t_pool  = {t_pooled:.15f}")
print(f"  두 통계량이 같은가: {t_welch == t_pooled}")

nu = (a + b) ** 2 / (a**2 / (n - 1) + b**2 / (n - 1))
rho = v2 / v1
print(f"\n자유도: 합동 {2 * n - 2},  Welch {nu:.9f}")
print(f"  균형 공식 (n-1)(1+rho)^2/(1+rho^2) = "
      f"{(n - 1) * (1 + rho) ** 2 / (1 + rho**2):.9f}   (rho = s2^2/s1^2 = {rho:.6f})")
print(f"  공식이 주는 범위: [{n - 1}, {2 * (n - 1)}]")
print(f"\n같은 |t| = {abs(t_welch):.6f} 를 두 자유도로 읽으면")
print(f"  df = {2 * n - 2}      : p = {2 * stats.t.sf(abs(t_welch), 2 * n - 2):.9f}")
print(f"  df = {nu:.4f} : p = {2 * stats.t.sf(abs(t_welch), nu):.9f}")

n1 = n2 = 10
  SE_welch = 5.217491947500
  SE_pool  = 5.217491947500
  차 = 8.882e-16   (0 이어야 한다)
  t_welch = -1.034980035299904
  t_pool  = -1.034980035299904
  두 통계량이 같은가: True

자유도: 합동 18,  Welch 12.424624520
  균형 공식 (n-1)(1+rho)^2/(1+rho^2) = 12.424624520   (rho = s2^2/s1^2 = 0.197693)
  공식이 주는 범위: [9, 18]

같은 |t| = 1.034980 를 두 자유도로 읽으면
  df = 18      : p = 0.314383791
  df = 12.4246 : p = 0.320403323

import statsmodels.api as sm

# statsmodels는 자유도까지 함께 돌려준다. scipy는 그렇지 않다.
# 결과를 보고할 때 자유도를 함께 적어야 하므로 이 점이 편하다.
t_stat, p_value, df = sm.stats.ttest_ind(group1, group2,
                                         usevar='unequal',
                                         alternative='two-sided')
print(f"t = {t_stat:.4f}, p = {p_value:.4f}, df = {df:.4f}")

t = -1.0350, p = 0.3204, df = 12.4246

# 손으로 구한 값과 맞추어 본다.
nu_hand = (n - 1) * (1 + rho) ** 2 / (1 + rho**2)
print(f"statsmodels df = {df:.12f}")
print(f"손으로 구한 df = {nu_hand:.12f}   차 = {abs(df - nu_hand):.2e}")
print(f"분산비 s1^2/s2^2 = {v1 / v2:.6f}  (본문의 '다섯 배')")
print(f"잃은 자유도 = {2 * (n - 1) - nu_hand:.6f}  "
      f"(합동 {2 * (n - 1)} 의 {100 * nu_hand / (2 * (n - 1)):.1f}% 만 남는다)")

print("\n분산비를 바꾸면 자유도가 어떻게 되는가 (n = 10 균형)")
print("   s1^2/s2^2     nu      합동과의 비")
for ratio in (1, 2, 5.058358, 10, 100, 10_000):
    rr = 1 / ratio
    nn = (n - 1) * (1 + rr) ** 2 / (1 + rr**2)
    print(f"{ratio:11.4f}  {nn:8.4f}     {nn / (2 * (n - 1)):.4f}")

print(f"\nscipy 와 statsmodels 의 t, p 가 같은가: "
      f"{np.allclose([t_stat, p_value], [t_welch, p_welch])}")

statsmodels df = 12.424624519556
손으로 구한 df = 12.424624519556   차 = 0.00e+00
분산비 s1^2/s2^2 = 5.058358  (본문의 '다섯 배')
잃은 자유도 = 5.575375  (합동 18 의 69.0% 만 남는다)

분산비를 바꾸면 자유도가 어떻게 되는가 (n = 10 균형)
   s1^2/s2^2     nu      합동과의 비
     1.0000   18.0000     1.0000
     2.0000   16.2000     0.9000
     5.0584   12.4246     0.6903
    10.0000   10.7822     0.5990
   100.0000    9.1800     0.5100
 10000.0000    9.0018     0.5001

scipy 와 statsmodels 의 t, p 가 같은가: True

import numpy as np
from scipy import stats

rng = np.random.default_rng(717)
M, n = 5_000, 20

def perm_p(x, y, B, rng):
    """평균 차이를 통계량으로 하는 양측 순열검정."""
    obs = abs(x.mean() - y.mean())
    z = np.r_[x, y]
    cnt = 1
    for _ in range(B):
        rng.shuffle(z)
        cnt += abs(z[:len(x)].mean() - z[len(x):].mean()) >= obs
    return cnt / (B + 1)

print(f"{'분포':>10s} {'상황':>6s} {'Welch':>8s} {'MWU':>8s} {'순열':>8s}")
for name in ["정규", "로그정규", "t(3)", "지수"]:
    for shift, lab in [(0.0, "수준"), (0.8, "검정력")]:
        a = b = c = 0
        for _ in range(M):
            if name == "정규":
                x, y = rng.standard_normal(n), rng.standard_normal(n) + shift
            elif name == "로그정규":
                x, y = rng.lognormal(0, 1, n), rng.lognormal(0, 1, n) + shift
            elif name == "t(3)":
                x, y = rng.standard_t(3, n), rng.standard_t(3, n) + shift
            else:
                x, y = rng.exponential(1, n), rng.exponential(1, n) + shift
            a += stats.ttest_ind(x, y, equal_var=False).pvalue < 0.05
            b += stats.mannwhitneyu(x, y).pvalue < 0.05
            c += perm_p(x, y, 299, rng) < 0.05
        print(f"{name:>10s} {lab:>6s} {a / M:8.4f} {b / M:8.4f} {c / M:8.4f}")

        분포     상황    Welch      MWU       순열
        정규     수준   0.0448   0.0436   0.0404
        정규    검정력   0.6868   0.6618   0.6756
      로그정규     수준   0.0368   0.0480   0.0462
      로그정규    검정력   0.3226   0.7184   0.3424
      t(3)     수준   0.0412   0.0474   0.0450
      t(3)    검정력   0.3630   0.4812   0.3674
        지수     수준   0.0490   0.0524   0.0516
        지수    검정력   0.7036   0.8772   0.7000

rng = np.random.default_rng(818)
M, n = 5_000, 30
mu2 = 0.5 - 1.4**2 / 2                   # 두 평균이 같아지도록 맞춘다
print(f"X: 로그정규(0, 1)     평균 {np.exp(0.5):.4f}  중앙값 {1.0:.4f}")
print(f"Y: 로그정규({mu2:.2f}, 1.4)  평균 {np.exp(mu2 + 1.4**2 / 2):.4f}  "
      f"중앙값 {np.exp(mu2):.4f}")

a = b = 0
for _ in range(M):
    x, y = rng.lognormal(0, 1, n), rng.lognormal(mu2, 1.4, n)
    a += stats.ttest_ind(x, y, equal_var=False).pvalue < 0.05
    b += stats.mannwhitneyu(x, y).pvalue < 0.05
print(f"\n평균은 같고 중앙값이 다른 상황")
print(f"  Welch 기각률 {a / M:.4f}   MWU 기각률 {b / M:.4f}")

xs, ys = rng.lognormal(0, 1, 200_000), rng.lognormal(mu2, 1.4, 200_000)
print(f"  P(Y > X) ≈ {np.mean(ys > xs):.4f}")

X: 로그정규(0, 1)     평균 1.6487  중앙값 1.0000
Y: 로그정규(-0.48, 1.4)  평균 1.6487  중앙값 0.6188

평균은 같고 중앙값이 다른 상황
  Welch 기각률 0.0720   MWU 기각률 0.3160
  P(Y > X) ≈ 0.3886

import numpy as np
from scipy import stats

x1, s1, n1 = 55_000.0, 7_500.0, 14      # 프로그램 A
x2, s2, n2 = 60_000.0, 8_000.0, 16      # 프로그램 B
diff = x1 - x2

# ── 합동 t
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 = stats.t.ppf(0.975, df_p)
print(f"합동  sp={sp:.2f}  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 * se_p:.1f}, {diff + tc * se_p:.1f})")

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

# ── 효과크기
d = diff / sp
g = d * (1 - 3 / (4 * (n1 + n2) - 9))            # 헤지스의 소표본 보정
se_d = np.sqrt((n1 + n2) / (n1 * n2) + d**2 / (2 * (n1 + n2 - 2)))
print(f"\nCohen d = {d:.4f},  Hedges g = {g:.4f}")
print(f"  d의 95% CI ({d - 1.96 * se_d:.4f}, {d + 1.96 * se_d:.4f})")

# ── 검정력
def power_ind(n, d, alpha=0.05):
    nc = abs(d) * np.sqrt(n / 2)
    df = 2 * n - 2
    c = stats.t.ppf(1 - alpha / 2, df)
    return stats.nct.sf(c, df, nc) + stats.nct.cdf(-c, df, nc)

print(f"\n군당 15명일 때 참 d={d:.3f}에 대한 검정력 {power_ind(15, d):.4f}")
for n in range(5, 2_000):
    if power_ind(n, d) >= 0.80:
        print(f"80% 검정력에 필요한 군당 표본 {n}명")
        break

합동  sp=7771.86  t=-1.7580  df=28  p=0.0897
      95% CI (-10826.1, 826.1)
Welch t=-1.7658  df=27.85  p=0.0884
      95% CI (-10801.7, 801.7)

Cohen d = -0.6433,  Hedges g = -0.6260
  d의 95% CI (-1.3802, 0.0935)

군당 15명일 때 참 d=-0.643에 대한 검정력 0.3981
80% 검정력에 필요한 군당 표본 39명