A/B 검정 순열검정 (코드)¶

import numpy as np

# 앞 10개가 A 페이지, 뒤 10개가 B 페이지의 체류시간이다.
times = np.array([185, 188, 142, 160, 161, 157, 182, 181, 159, 167,   # A 페이지
                  173, 181, 182, 170, 169, 177, 168, 183, 169, 164])  # B 페이지

def perm_test_two_sample_means(data, nA, n_perms=9999, rng=None):
    """두 집단의 평균 차이에 대한 순열검정.

    귀무가설이 참이라면 어느 값이 A 에서 나왔고 어느 값이 B 에서 나왔는지가
    아무 뜻이 없다. 그래서 이름표를 마구 섞어 가며 차이를 다시 계산하면,
    "차이가 없을 때 이 정도 차이가 얼마나 흔한가"를 직접 셀 수 있다.
    """
    rng = rng or np.random.default_rng(0)
    obs_diff = data[:nA].mean() - data[nA:].mean()
    perm_diffs = np.empty(n_perms)
    for i in range(n_perms):
        p = rng.permutation(data)
        perm_diffs[i] = p[:nA].mean() - p[nA:].mean()
    p_value = ((np.abs(perm_diffs) >= abs(obs_diff)).sum() + 1) / (n_perms + 1)
    return p_value, perm_diffs, obs_diff

import itertools
from scipy import stats

p_value, perm_diffs, obs_diff = perm_test_two_sample_means(times, 10)
print(f"관측 차이 = {obs_diff:.2f},  순열 p (B=9999) = {p_value:.4f}")
print(f"Welch t 검정 p = "
      f"{stats.ttest_ind(times[:10], times[10:], equal_var=False).pvalue:.4f}")

# C(20,10) = 184,756 가지를 모두 열거한다.
z = times.astype(float)
T = z.sum()
d = np.array([(2 * z[list(c)].sum() - T) / 10
              for c in itertools.combinations(range(20), 10)])
ad = np.abs(d)
print(f"\n정확 p = {(ad >= abs(obs_diff) - 1e-9).mean():.6f}"
      f"   (SD(d*) = {d.std():.6f},  1.96 x SD = {1.96 * d.std():.4f})")

# alpha = 0.05 기각 문턱을 격자 위에서 찾는다.
for v in np.unique(np.round(ad, 6)):
    if (ad >= v - 1e-9).mean() <= 0.05:
        break
print(f"문턱 |d| = {v:.1f} 의 정확 p   = {(ad >= v - 1e-9).mean():.6f}")
print(f"한 칸 아래 {v - 0.2:.1f} 의 정확 p = {(ad >= v - 0.2 - 1e-9).mean():.6f}")
print(f"관측 |d| / 문턱 = {abs(obs_diff) / v:.4f}")

관측 차이 = -5.40,  순열 p (B=9999) = 0.3306
Welch t 검정 p = 0.3204

정확 p = 0.325586   (SD(d*) = 5.227257,  1.96 x SD = 10.2454)
문턱 |d| = 10.2 의 정확 p   = 0.049352
한 칸 아래 10.0 의 정확 p = 0.054580
관측 |d| / 문턱 = 0.5294

def perm_test_proportion(n_control, conv_control, n_treatment, conv_treatment,
                         n_perms=9999, rng=None):
    """전환율 차이에 대한 순열검정.

    이진 자료도 다를 것이 없다. 전체 전환 수만큼 1 을 채운 배열을 만들고
    그것을 섞으면, "전환이 두 집단에 무작위로 흩어진" 상태가 된다.
    """
    rng = rng or np.random.default_rng(0)
    binary = np.zeros(n_control + n_treatment)
    binary[:conv_control + conv_treatment] = 1
    obs_diff = conv_treatment / n_treatment - conv_control / n_control

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

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

from scipy import stats

p_value, perm_diffs, obs_diff = perm_test_proportion(23739, 200, 22588, 182,
                                                     n_perms=2000)
print(f"관측 차이 = {obs_diff:.6f},  순열 p (B=2000) = {p_value:.4f}")
print(f"Fisher 정확검정 = "
      f"{stats.fisher_exact([[200, 23739 - 200], [182, 22588 - 182]])[1]:.6f}")

n_c, c_c, n_t, c_t = 23739, 200, 22588, 182
N, K = n_c + n_t, c_c + c_t
step = 1 / n_t + 1 / n_c
x = np.arange(K + 1)
d = x * step - K / n_c
pmf = stats.hypergeom.pmf(x, N, K, n_t)

def p_exact(x_obs):
    d_obs = x_obs * step - K / n_c
    return pmf[np.abs(d) >= abs(d_obs) - 1e-18].sum()

EX = K * n_t / N
sd_null = np.sqrt(N / (N - 1) * (K / N) * (1 - K / N) * step)
print(f"\n전환 한 건이 움직이는 비율차 = {step:.3e},  E[X] = {EX:.2f}")
lo = max(v for v in range(0, 187) if p_exact(v) <= 0.05)
hi = min(v for v in range(187, K + 1) if p_exact(v) <= 0.05)
print(f"기각역: X <= {lo} (p = {p_exact(lo):.4f})  또는  X >= {hi} (p = {p_exact(hi):.4f})")
print(f"  경계 바로 안쪽: X = {lo + 1} 이면 p = {p_exact(lo + 1):.4f},"
      f"  X = {hi - 1} 이면 p = {p_exact(hi - 1):.4f}")
print(f"정규근사 문턱 = {1.96 * sd_null / step:.2f} 건"
      f"  ->  X <= {EX - 1.96 * sd_null / step:.2f} 또는 X >= {EX + 1.96 * sd_null / step:.2f}")
print(f"관측 X = {c_t}: 아래 경계까지 {c_t - lo} 건,  위 경계까지 {hi - c_t} 건")
print(f"정확 p (X = {c_t}) = {p_exact(c_t):.6f}")

관측 차이 = -0.000368,  순열 p (B=2000) = 0.6667
Fisher 정확검정 = 0.681128

전환 한 건이 움직이는 비율차 = 8.640e-05,  E[X] = 186.25
기각역: X <= 166 (p = 0.0397)  또는  X >= 206 (p = 0.0450)
  경계 바로 안쪽: X = 167 이면 p = 0.0508,  X = 205 이면 p = 0.0572
정규근사 문턱 = 19.07 건  ->  X <= 167.19 또는 X >= 205.32
관측 X = 182: 아래 경계까지 16 건,  위 경계까지 24 건
정확 p (X = 182) = 0.681128

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

def power_estimate(n, true_diff=15, sigma=20, M=1000, B=999):
    reject = 0
    for _ in range(M):
        a = rng.normal(170, sigma, n)
        b = rng.normal(170 + true_diff, sigma, n)
        obs = a.mean() - b.mean()
        z = np.concatenate([a, b])
        P = np.array([rng.permutation(z) for _ in range(B)])
        d = P[:, :n].mean(1) - P[:, n:].mean(1)
        reject += ((np.abs(d) >= abs(obs)).sum() + 1)/(B+1) < 0.05
    return reject / M

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

times = np.array([185, 188, 142, 160, 161, 157, 182, 181, 159, 167,
                  173, 181, 182, 170, 169, 177, 168, 183, 169, 164])

obs_mean = times[:10].mean() - times[10:].mean()
obs_med = np.median(times[:10]) - np.median(times[10:])
print(obs_mean, obs_med)          # -5.4  -7.5

B = 99999
P = np.array([rng.permutation(times) for _ in range(B)])
pm  = P[:, :10].mean(1) - P[:, 10:].mean(1)
pmd = np.median(P[:, :10], axis=1) - np.median(P[:, 10:], axis=1)

p_mean = ((np.abs(pm)  >= abs(obs_mean) - 1e-12).sum() + 1) / (B + 1)
p_med  = ((np.abs(pmd) >= abs(obs_med)  - 1e-12).sum() + 1) / (B + 1)

-5.400000000000006 -7.5