이표본 순열검정¶

import numpy as np
import matplotlib.pyplot as plt

# A 페이지와 B 페이지의 체류시간
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])

obs_diff = np.abs(page_a.mean() - page_b.mean())
print(f"Page A: mean = {page_a.mean():.2f}")    # 168.20
print(f"Page B: mean = {page_b.mean():.2f}")    # 173.60
print(f"Observed |difference|: {obs_diff:.2f}") # 5.40

Page A: mean = 168.20
Page B: mean = 173.60
Observed |difference|: 5.40

def two_sample_permutation_test(x, y, n_perms=1000, seed=0):
    """
    Two-sample permutation test for a difference in means.

    Parameters
    ----------
    x, y : array-like
        The two samples.
    n_perms : int
        Number of random permutations.

    Returns
    -------
    p_value : float
        Two-sided p-value, with the +1 correction.
    perm_diffs : ndarray
        The permutation distribution of |mean difference|.
    """
    rng = np.random.default_rng(seed)
    obs_diff = np.abs(x.mean() - y.mean())
    pooled = np.concatenate([x, y])
    nx = len(x)

    perm_diffs = np.empty(n_perms)
    for i in range(n_perms):
        p = rng.permutation(pooled)
        perm_diffs[i] = np.abs(p[:nx].mean() - p[nx:].mean())

    p_value = ((perm_diffs >= obs_diff).sum() + 1) / (n_perms + 1)
    return p_value, perm_diffs

p_val, perm_stats = two_sample_permutation_test(page_a, page_b, n_perms=10000)

print(f"Permutation test p-value: {p_val:.4f}")
print(f"Conclusion: {'Reject H0' if p_val < 0.05 else 'Fail to reject H0'}")

Permutation test p-value: 0.3307
Conclusion: Fail to reject H0

import itertools
from scipy import stats

z = np.concatenate([page_a, page_b]).astype(float)
n1, N, total = len(page_a), len(z), float(np.sum(page_a) + np.sum(page_b))

# 가능한 라벨 배정 184,756 가지를 모두 열거한다. d* 는 부호까지 담는다.
signed = np.array([(2 * z[list(c)].sum() - total) / n1
                   for c in itertools.combinations(range(N), n1)])
exact = np.abs(signed)
p_exact = (exact >= obs_diff - 1e-9).mean()
print(f"순열 수 = {len(signed)},  정확 p = {p_exact:.6f}")

B = 10_000
mean_hat = (B * p_exact + 1) / (B + 1)
bias = (1 - p_exact) / (B + 1)
sd = np.sqrt(B * p_exact * (1 - p_exact)) / (B + 1)
print(f"E[p-hat] = {mean_hat:.6f}   (치우침 +{bias:.3e},  SD {sd:.6f},"
      f"  치우침/SD = {bias / sd:.4f})")
print(f"실행값   = {p_val:.6f}   -> 예측에서 {(p_val - mean_hat) / sd:+.2f} SD")

sd_null = z.std(ddof=1) * np.sqrt(1 / n1 + 1 / (N - n1))
print(f"\nSD(d*) 닫힌 꼴      = {sd_null:.6f}")
print(f"열거한 d* 의 SD     = {signed.std(ddof=0):.6f}"
      f"   (왜도 {stats.skew(signed):.1e}, 초과첨도 {stats.kurtosis(signed):+.4f})")
print(f"E[T] 반정규 예측    = {sd_null * np.sqrt(2 / np.pi):.6f}")
print(f"E[T] 정확 열거      = {exact.mean():.6f}   (SD 의 {exact.mean() / sd_null:.6f} 배)")
print(f"E[T] 모의 (B=10000) = {perm_stats.mean():.6f}"
      f"   (평균의 몬테카를로 오차 {perm_stats.std(ddof=1) / np.sqrt(B):.6f})")

순열 수 = 184756,  정확 p = 0.325586
E[p-hat] = 0.325654   (치우침 +6.743e-05,  SD 0.004685,  치우침/SD = 0.0144)
실행값   = 0.330667   -> 예측에서 +1.07 SD

SD(d*) 닫힌 꼴      = 5.227257
열거한 d* 의 SD     = 5.227257   (왜도 -8.2e-17, 초과첨도 -0.3076)
E[T] 반정규 예측    = 4.170747
E[T] 정확 열거      = 4.234476   (SD 의 0.810076 배)
E[T] 모의 (B=10000) = 4.262900   (평균의 몬테카를로 오차 0.030767)

fig, ax = plt.subplots(figsize=(10, 6))

ax.hist(perm_stats, bins=30, alpha=0.7, color='steelblue', edgecolor='black',
        label='Permuted test statistics')

obs_stat = np.abs(page_a.mean() - page_b.mean())
ax.axvline(obs_stat, color='red', linewidth=2,
           label=f'Observed difference = {obs_stat:.2f}')

rejection_region = perm_stats[perm_stats >= obs_stat]
ax.hist(rejection_region, bins=30, alpha=0.5, color='red',
        label=f'P-value region (p = {p_val:.3f})')

ax.set_xlabel('Absolute Difference in Means')
ax.set_ylabel('Frequency')
ax.set_title('Two-Sample Permutation Test Distribution')
ax.legend()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

plt.tight_layout()
plt.show()

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

# 대조군: 23,739명 중 200명 전환
# 실험군: 22,588명 중 182명 전환
n_control, conv_control = 23739, 200
n_treatment, conv_treat = 22588, 182

rate_control = conv_control / n_control
rate_treatment = conv_treat / n_treatment
obs_diff_rates = rate_treatment - rate_control

print(f"Control conversion rate:   {rate_control:.4f}")     # 0.0084
print(f"Treatment conversion rate: {rate_treatment:.4f}")   # 0.0081
print(f"Observed difference: {obs_diff_rates:.6f}")         # -0.000368

# 이진 반응벡터. 앞의 (대조군 전환 + 실험군 전환) 개가 1 이다
binary_response = np.zeros(n_control + n_treatment)
binary_response[:conv_control + conv_treat] = 1

B = 5000
perm_diffs = np.empty(B)
for b in range(B):
    p = rng.permutation(binary_response)
    perm_diffs[b] = p[:n_treatment].mean() - p[n_treatment:].mean()

p_value_ab = ((np.abs(perm_diffs) >= abs(obs_diff_rates)).sum() + 1) / (B + 1)
print(f"A/B test p-value: {p_value_ab:.4f}")   # 0.68

Control conversion rate:   0.0084
Treatment conversion rate: 0.0081
Observed difference: -0.000368
A/B test p-value: 0.6785

from scipy import stats

N = n_control + n_treatment
K = conv_control + conv_treat
pbar = K / N
step = 1 / n_treatment + 1 / n_control
sd_null = np.sqrt(N / (N - 1) * pbar * (1 - pbar) * step)

mde = stats.norm.ppf(0.975) * sd_null
print(f"SD(d*)        = {sd_null:.8f}")
print(f"탐지 가능한 최소 |차| = {mde:.8f}"
      f"   (기저 전환율의 {100 * mde / pbar:.2f}%)")
print(f"관측 |차| / 문턱      = {abs(obs_diff_rates) / mde:.4f}")

# 정확 순열 p. 실험군 전환 수가 Hypergeometric(N, K, n_treatment) 을 따른다.
x = np.arange(K + 1)
d = x * step - K / n_control
pmf = stats.hypergeom.pmf(x, N, K, n_treatment)
p_exact = pmf[np.abs(d) >= abs(obs_diff_rates) - 1e-18].sum()

mean_hat = (B * p_exact + 1) / (B + 1)
bias = (1 - p_exact) / (B + 1)
sd = np.sqrt(B * p_exact * (1 - p_exact)) / (B + 1)
print(f"\n정확 순열 p   = {p_exact:.6f}")
print(f"E[p-hat]      = {mean_hat:.6f}   (치우침 +{bias:.3e},  SD {sd:.6f},"
      f"  치우침/SD = {bias / sd:.4f})")
print(f"실행값        = {p_value_ab:.6f}   -> 예측에서 {(p_value_ab - mean_hat) / sd:+.2f} SD")

SD(d*)        = 0.00084056
탐지 가능한 최소 |차| = 0.00164747   (기저 전환율의 19.98%)
관측 |차| / 문턱      = 0.2231

정확 순열 p   = 0.681128
E[p-hat]      = 0.681192   (치우침 +6.376e-05,  SD 0.006589,  치우침/SD = 0.0097)
실행값        = 0.678464   -> 예측에서 -0.41 SD

from scipy import stats

# 앞의 순열검정 결과를 Welch t 검정과 견준다. 자료가 정규에 가깝고
# 표본이 넉넉하면 두 p-값이 거의 같게 나온다. 순열검정이 t 검정을
# 대신하는 것이 아니라, 가정이 미덥지 않을 때 기댈 곳이 된다는 뜻이다.
t_stat, p_ttest = stats.ttest_ind(page_a, page_b, equal_var=False)
print(f"Welch's t-test p-value: {p_ttest:.4f}")
print(f"Permutation test p-value: {p_val:.4f}")

Welch's t-test p-value: 0.3204
Permutation test p-value: 0.3307

import itertools
import numpy as np

z = np.concatenate([page_a, page_b]).astype(float)
n1, n2 = len(page_a), len(page_b)
N, T, Q = n1 + n2, z.sum(), (z ** 2).sum()

t_obs = stats.ttest_ind(page_a, page_b).statistic
sp2 = ((n1 - 1) * page_a.var(ddof=1) + (n2 - 1) * page_b.var(ddof=1)) / (N - 2)
S2 = z.var(ddof=1)
print(f"s_p^2 = {sp2:.6f},   S^2 = {S2:.6f},   t_obs = {t_obs:.6f}")
print(f"항등식  (N-1)S^2      = {(N - 1) * S2:.6f}")
print(f"        (N-2+t^2)s_p^2 = {(N - 2 + t_obs ** 2) * sp2:.6f}")
print(f"SD(d*)/SE_t  예측 = {np.sqrt((N - 2 + t_obs ** 2) / (N - 1)):.8f}")
print(f"SD(d*)/SE_t  실제 = {np.sqrt(S2 / sp2):.8f}")

# 모든 배정에서 |d*| 와 |t*| 를 함께 구한다.
# 둘 다 A 집단의 합 SA 하나로 정해진다.
SA = np.array([z[list(c)].sum() for c in itertools.combinations(range(N), n1)])
d = SA / n1 - (T - SA) / n2
sp2s = (Q - SA ** 2 / n1 - (T - SA) ** 2 / n2) / (N - 2)
t = d / np.sqrt(sp2s * (1 / n1 + 1 / n2))

obs_d = abs(page_a.mean() - page_b.mean())
print(f"\n정확 순열 p (|d*| 기준) = {(np.abs(d) >= obs_d - 1e-9).mean():.6f}")
print(f"정확 순열 p (|t*| 기준) = {(np.abs(t) >= abs(t_obs) - 1e-9).mean():.6f}")
order = np.argsort(np.abs(d))
print(f"|t*| 가 |d*| 의 증가함수인가: "
      f"{bool(np.all(np.diff(np.abs(t)[order]) >= -1e-9))}")
print(f"두 기준이 고르는 배정 집합이 같은가: "
      f"{np.array_equal(np.abs(d) >= obs_d - 1e-9, np.abs(t) >= abs(t_obs) - 1e-9)}")

w = stats.ttest_ind(page_a, page_b, equal_var=False)
print(f"\nWelch  t 검정 p = {w.pvalue:.6f}   (df = {w.df:.4f},  SE = "
      f"{np.sqrt(page_a.var(ddof=1) / n1 + page_b.var(ddof=1) / n2):.6f})")
print(f"합동   t 검정 p = {stats.ttest_ind(page_a, page_b).pvalue:.6f}"
      f"   (df = {N - 2},  SE = {np.sqrt(sp2 * (1 / n1 + 1 / n2)):.6f})")
print(f"순열   (B=10000) = {p_val:.6f}   (귀무분포 SD = {np.sqrt(S2 * (1 / n1 + 1 / n2)):.6f})")

s_p^2 = 136.111111,   S^2 = 136.621053,   t_obs = -1.034980
항등식  (N-1)S^2      = 2595.800000
        (N-2+t^2)s_p^2 = 2595.800000
SD(d*)/SE_t  예측 = 1.00187150
SD(d*)/SE_t  실제 = 1.00187150

정확 순열 p (|d*| 기준) = 0.325586
정확 순열 p (|t*| 기준) = 0.325586
|t*| 가 |d*| 의 증가함수인가: True
두 기준이 고르는 배정 집합이 같은가: True

Welch  t 검정 p = 0.320403   (df = 12.4246,  SE = 5.217492)
합동   t 검정 p = 0.314384   (df = 18,  SE = 5.217492)
순열   (B=10000) = 0.330667   (귀무분포 SD = 5.227257)

import numpy as np, itertools
from scipy import stats

A = np.array([14.2, 16.8, 13.5, 15.9, 17.3, 12.8, 16.1, 14.7])
B = np.array([11.3, 13.6, 10.9, 12.4, 14.1, 11.8, 13.2, 12.7])
z = np.concatenate([A, B]); n = 8

obs_mean = A.mean() - B.mean()
obs_med = np.median(A) - np.median(B)
print(round(obs_mean, 4), round(obs_med, 4))   # 2.6625  2.75

dm, dmd = [], []
for c in itertools.combinations(range(16), n):
    a = z[list(c)]; b = np.delete(z, list(c))
    dm.append(a.mean() - b.mean())
    dmd.append(np.median(a) - np.median(b))
dm, dmd = np.array(dm), np.array(dmd)

print("perms:", len(dm))                                    # 12870
print("mean  :", (np.abs(dm)  >= abs(obs_mean) - 1e-12).mean())
print("median:", (np.abs(dmd) >= abs(obs_med)  - 1e-12).mean())
print(stats.ttest_ind(A, B, equal_var=False))

2.6625 2.75
perms: 12870
mean  : 0.002641802641802642
median: 0.003108003108003108
TtestResult(statistic=3.837370691851344, pvalue=0.0022024857203016353, df=12.494270083405828)

import numpy as np
rng = np.random.default_rng(11)
A = np.array([14.2, 16.8, 13.5, 15.9, 17.3, 12.8, 16.1, 14.7])
B = np.array([11.3, 13.6, 10.9, 12.4, 14.1, 11.8, 13.2, 12.7])
z = np.concatenate([A, B]); n = 8
obs = abs(A.mean() - B.mean())

for Bp in (200, 1000, 10000):
    ps = []
    for _ in range(300):
        P = np.array([rng.permutation(z) for _ in range(Bp)])
        d = P[:, :n].mean(1) - P[:, n:].mean(1)
        ps.append(((np.abs(d) >= obs - 1e-12).sum() + 1) / (Bp + 1))
    ps = np.array(ps)
    print(Bp, round(ps.mean(), 5), round(ps.std(), 5),
          round((ps < 0.05).mean(), 3))

200 0.00789 0.00401 1.0
1000 0.00363 0.00155 1.0
10000 0.00274 0.00051 1.0

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

def perm_p(x, y, stat, Bp=499):
    obs = abs(stat(x) - stat(y))
    z = np.concatenate([x, y]); m = len(x)
    c = 0
    for _ in range(Bp):
        p = rng.permutation(z)
        c += abs(stat(p[:m]) - stat(p[m:])) >= obs - 1e-12
    return (c + 1) / (Bp + 1)

trim = lambda a: stats.trim_mean(a, 0.2)

def power(gen, shift, M=800, m=20, n=20):
    r = dict(mean=0, median=0, trim=0, welch=0)
    for _ in range(M):
        x = gen(m) + shift; y = gen(n)
        r['mean']   += perm_p(x, y, np.mean)   < 0.05
        r['median'] += perm_p(x, y, np.median) < 0.05
        r['trim']   += perm_p(x, y, trim)      < 0.05
        r['welch']  += stats.ttest_ind(x, y, equal_var=False).pvalue < 0.05
    return {k: round(v/M, 3) for k, v in r.items()}

norm = lambda k: rng.normal(0, 1, k)
def contam(k):
    a = rng.normal(0, 1, k)
    mask = rng.random(k) < 0.1
    a[mask] = rng.normal(0, 5, mask.sum())
    return a

for name, gen in [("normal", norm), ("contaminated", contam)]:
    for sh in (0.0, 1.0):
        print(name, sh, power(gen, sh))

normal 0.0 {'mean': 0.041, 'median': 0.055, 'trim': 0.045, 'welch': 0.039}
normal 1.0 {'mean': 0.866, 'median': 0.789, 'trim': 0.838, 'welch': 0.87}
contaminated 0.0 {'mean': 0.056, 'median': 0.056, 'trim': 0.048, 'welch': 0.042}
contaminated 1.0 {'mean': 0.47, 'median': 0.664, 'trim': 0.712, 'welch': 0.464}