순열검정: 기초¶

import pandas as pd
import numpy as np
import random

# 예시 자료. "Practical Statistics for Data Scientists" 에서 가져왔다
session_times = pd.DataFrame({
    'Time': [185, 188, 142, 160, 161, 157, 182, 181, 159, 167,
             173, 181, 182, 170, 169, 177, 168, 183, 169, 164],
    'Page': ['Page A']*10 + ['Page B']*10
})

mean_a = session_times[session_times.Page == 'Page A'].Time.mean()
mean_b = session_times[session_times.Page == 'Page B'].Time.mean()
observed_diff = mean_b - mean_a

print(f"Page A mean: {mean_a:.2f} seconds")        # 168.20
print(f"Page B mean: {mean_b:.2f} seconds")        # 173.60
print(f"Observed difference: {observed_diff:.2f}") # 5.40

Page A mean: 168.20 seconds
Page B mean: 173.60 seconds
Observed difference: 5.40

def perm_fun(x, nA, nB):
    """
    Randomly shuffle group labels and compute difference of means.

    Parameters
    ----------
    x : pandas Series with a 0..n-1 integer index
    nA, nB : group sizes

    Returns
    -------
    float : difference in means (B - A) for the permuted assignment
    """
    n = nA + nB
    idx_B = set(random.sample(range(n), nB))
    idx_A = set(range(n)) - idx_B
    return x.loc[list(idx_B)].mean() - x.loc[list(idx_A)].mean()

nA = session_times[session_times.Page == 'Page A'].shape[0]
nB = session_times[session_times.Page == 'Page B'].shape[0]

random.seed(42)
perm_diffs = [perm_fun(session_times.Time, nA, nB) for _ in range(1000)]

p_value = np.mean(np.abs(perm_diffs) >= np.abs(observed_diff))
print(f"Permutation test p-value: {p_value:.4f}")   # 0.3310

Permutation test p-value: 0.3310

from scipy import stats

z = session_times.Time.to_numpy(float)
S = z.std(ddof=1)                          # 합친 20개의 표본표준편차
sd_theory = S * np.sqrt(1 / nA + 1 / nB)   # (1) 에서 유도한 닫힌 꼴

print(f"T = {z.sum():.0f},  xbar = {z.mean():.1f},  S^2 = {z.var(ddof=1):.5f}")
print(f"SD(d*) 닫힌 꼴      = {sd_theory:.6f}")
print(f"B=1000 순열분포의 SD = {np.std(perm_diffs, ddof=1):.6f}"
      f"   (몬테카를로 오차 {sd_theory / np.sqrt(2 * 1000):.6f})")

gaps = np.diff(np.unique(np.round(perm_diffs, 6)))
print(f"순열값의 최소 간격   = {gaps.min():.6f}")

zobs = observed_diff / sd_theory
print(f"\n관측 d / SD        = {zobs:.5f}")
print(f"정규근사 양측 p     = {2 * stats.norm.sf(zobs):.4f}")
print(f"연속성보정 (반 칸 0.1) = {2 * stats.norm.sf((observed_diff - 0.1) / sd_theory):.4f}")
print(f"순열 p (B = 1000)   = {p_value:.4f}")
print(f"t 검정 (합동분산) p  = {stats.ttest_ind(z[10:], z[:10]).pvalue:.4f}")

T = 3418,  xbar = 170.9,  S^2 = 136.62105
SD(d*) 닫힌 꼴      = 5.227257
B=1000 순열분포의 SD = 5.100405   (몬테카를로 오차 0.116885)
순열값의 최소 간격   = 0.200000

관측 d / SD        = 1.03305
정규근사 양측 p     = 0.3016
연속성보정 (반 칸 0.1) = 0.3106
순열 p (B = 1000)   = 0.3310
t 검정 (합동분산) p  = 0.3144

import numpy as np
rng = np.random.default_rng(0)

n_control, n_treat = 23739, 22588
c_control, c_treat = 200, 182

obs_diff = c_treat/n_treat - c_control/n_control
print(f"{obs_diff:.6f}")        # -0.000368

# 1 = converted, 0 = did not convert
total = n_control + n_treat
conversion = np.zeros(total)
conversion[:c_control + c_treat] = 1

B = 5000
perm_diffs = np.empty(B)
for b in range(B):
    perm = rng.permutation(conversion)
    perm_diffs[b] = perm[:n_treat].mean() - perm[n_treat:].mean()

p_value = np.mean(np.abs(perm_diffs) >= abs(obs_diff))
print(f"Conversion A/B test p-value: {p_value:.4f}")   # 0.6784

-0.000368
Conversion A/B test p-value: 0.6784

from scipy import stats

N, K = total, c_control + c_treat
pbar = K / N
step = 1 / n_treat + 1 / n_control             # 순열값의 격자 간격

S2 = N / (N - 1) * pbar * (1 - pbar)           # 합친 0/1 자료의 표본분산
sd_theory = np.sqrt(S2 * step)                 # (1) 에서 유도한 닫힌 꼴
sd_pooled = np.sqrt(pbar * (1 - pbar) * step)  # 합동비율을 쓴 고전 표준오차

print(f"N = {N},  K = {K},  pbar = {pbar:.8f},  격자 간격 = {step:.4e}")
print(f"SD(d*) 닫힌 꼴    = {sd_theory:.10f}")
print(f"고전 합동 표준오차 = {sd_pooled:.10f}   (비 = {sd_theory / sd_pooled:.8f})")
print(f"B=5000 순열분포 SD = {perm_diffs.std(ddof=1):.10f}"
      f"   (몬테카를로 오차 {sd_theory / np.sqrt(2 * B):.10f})")

zobs = obs_diff / sd_theory
print(f"\n관측 d / SD     = {zobs:.5f}")
print(f"정규근사 양측 p  = {2 * stats.norm.sf(abs(zobs)):.4f}")

# 정확 순열분포. 처치군 전환 수 X ~ Hypergeometric(N, K, n_treat) 이고
# d 는 X 의 증가하는 일차식이므로 X 의 분포만으로 정확 p 가 나온다.
x = np.arange(K + 1)
d = x * step - K / n_control
pmf = stats.hypergeom.pmf(x, N, K, n_treat)
sel = np.abs(d) >= abs(obs_diff) - 1e-18
x_lo = x[sel][d[sel] < 0].max()
x_hi = x[sel][d[sel] > 0].min()
print(f"정확 양측 순열 p = {pmf[sel].sum():.4f}"
      f"   (걸리는 X: <= {x_lo} 또는 >= {x_hi})")
print(f"순열 p (B=5000)  = {p_value:.4f}")
print(f"참고: 단측 상측 P(X >= {c_treat}) = "
      f"{stats.hypergeom.sf(c_treat - 1, N, K, n_treat):.4f}")

# 연속성보정. X 의 격자에서 양쪽으로 반 칸을 물려 잰다.
EX = K * n_treat / N
sdX = np.sqrt(n_treat * n_control * K * (N - K) / (N**2 * (N - 1)))
cc = stats.norm.cdf((x_lo + 0.5 - EX) / sdX) + stats.norm.sf((x_hi - 0.5 - EX) / sdX)
print(f"\nE[X] = {EX:.4f},  SD(X) = {sdX:.4f}")
print(f"연속성보정 정규근사 = {cc:.4f}")

N = 46327,  K = 382,  pbar = 0.00824573,  격자 간격 = 8.6396e-05
SD(d*) 닫힌 꼴    = 0.0008405592
고전 합동 표준오차 = 0.0008405501   (비 = 1.00001079)
B=5000 순열분포 SD = 0.0008365779   (몬테카를로 오차 0.0000084056)

관측 d / SD     = -0.43730
정규근사 양측 p  = 0.6619
정확 양측 순열 p = 0.6811   (걸리는 X: <= 182 또는 >= 191)
순열 p (B=5000)  = 0.6784
참고: 단측 상측 P(X >= 182) = 0.6873

E[X] = 186.2546,  SD(X) = 9.7291
연속성보정 정규근사 = 0.6811

import numpy as np
from scipy import stats
rng = np.random.default_rng(0)

def size(m, n, s1, s2, M=2000, B=999, studentized=False):
    rej = 0
    for _ in range(M):
        x = rng.normal(0, s1, m); y = rng.normal(0, s2, n)
        z = np.concatenate([x, y])
        if studentized:
            obs = (x.mean() - y.mean()) / np.sqrt(
                  x.var(ddof=1)/m + y.var(ddof=1)/n)
        else:
            obs = x.mean() - y.mean()
        cnt = 0
        for _ in range(B):
            p = rng.permutation(z)
            a, b = p[:m], p[m:]
            d = (a.mean() - b.mean())
            if studentized:
                d /= np.sqrt(a.var(ddof=1)/m + b.var(ddof=1)/n)
            cnt += abs(d) >= abs(obs)
        rej += (cnt + 1) / (B + 1) < 0.05
    return round(rej/M, 3)

for (m, n) in [(20, 20), (10, 30), (30, 10)]:
    print(m, n, size(m, n, 1, 3), size(m, n, 1, 3, studentized=True))

20 20 0.058 0.059
10 30 0.004 0.038
30 10 0.196 0.052

import numpy as np, itertools
rng = np.random.default_rng(1)
# 각 집단 6명, 처치군 5명 성공 / 대조군 2명 성공
x = np.array([1, 1, 1, 1, 1, 0])
y = np.array([1, 1, 0, 0, 0, 0])
obs = x.mean() - y.mean()

z = np.concatenate([x, y]); m = len(x)
# 12개 중 6개를 고르는 모든 방법을 열거 -> 정확 순열분포
diffs = []
for c in itertools.combinations(range(12), m):
    a = z[list(c)]; b = np.delete(z, list(c))
    diffs.append(a.mean() - b.mean())
diffs = np.array(diffs)

print("총 순열 수:", len(diffs))
print("p (>)   =", round((np.abs(diffs) >  abs(obs)).mean(), 5))
print("p (>=)  =", round((np.abs(diffs) >= abs(obs)).mean(), 5))

총 순열 수: 924
p (>)   = 0.01515
p (>=)  = 0.24242

import numpy as np
from scipy import stats
rng = np.random.default_rng(7)

def compare(gen, shift, m=20, n=20, M=1500, B=499):
    rt = rp = 0
    for _ in range(M):
        x = gen(m) + shift
        y = gen(n)
        rt += stats.ttest_ind(x, y).pvalue < 0.05
        obs = x.mean() - y.mean()
        z = np.concatenate([x, y])
        cnt = 0
        for _ in range(B):
            p = rng.permutation(z)
            cnt += abs(p[:m].mean() - p[m:].mean()) >= abs(obs)
        rp += (cnt + 1)/(B + 1) < 0.05
    return round(rt/M, 3), round(rp/M, 3)

norm = lambda k: rng.normal(0, 1, k)
t3   = lambda k: rng.standard_t(3, k) / np.sqrt(3)   # 분산 1 로 맞춤

for name, gen in [("normal", norm), ("t(3)", t3)]:
    for sh in (0.0, 0.6, 1.0):
        print(name, sh, compare(gen, sh))

normal 0.0 (0.045, 0.041)
normal 0.6 (0.463, 0.447)
normal 1.0 (0.869, 0.864)
t(3) 0.0 (0.047, 0.043)
t(3) 0.6 (0.543, 0.545)
t(3) 1.0 (0.884, 0.889)

import numpy as np
from scipy import stats
rng = np.random.default_rng(0)
n_control, n_treat = 23739, 22588
c_control, c_treat = 200, 182
total = n_control + n_treat
K = c_control + c_treat

# 몬테카를로 순열
conversion = np.zeros(total); conversion[:K] = 1
B = 20000
cnt = np.array([rng.permutation(conversion)[:n_treat].sum()
                for _ in range(B)])
obs = c_treat
p_mc = ((np.abs(cnt - K*n_treat/total) >=
         abs(obs - K*n_treat/total)).sum() + 1) / (B + 1)

# 초기하분포 정확값
p_exact = 2 * min(stats.hypergeom.cdf(obs, total, K, n_treat),
                  stats.hypergeom.sf(obs - 1, total, K, n_treat))

# Fisher 정확검정
p_fisher = stats.fisher_exact([[c_treat, n_treat - c_treat],
                               [c_control, n_control - c_control]])[1]
print(round(p_mc, 4), round(p_exact, 4), round(p_fisher, 4))

0.6813 0.6999 0.6811