순열검정 시연 (코드)¶

import numpy as np

def perm_test_two_sample(x, y, n_perm=9999, rng=None):
    """평균 차이에 대한 이표본 순열검정.

    두 집단을 합친 뒤 이름표를 섞어 차이를 다시 계산하는 일을 되풀이한다.
    분포를 가정하지 않으므로 정규성도 등분산도 필요 없다.
    """
    rng = rng or np.random.default_rng(0)
    obs_diff = x.mean() - y.mean()
    pooled = np.concatenate([x, y])
    m = len(x)
    P = np.array([rng.permutation(pooled) for _ in range(n_perm)])
    perm_diffs = P[:, :m].mean(axis=1) - P[:, m:].mean(axis=1)
    p_value = ((np.abs(perm_diffs) >= abs(obs_diff)).sum() + 1) / (n_perm + 1)
    return obs_diff, p_value, perm_diffs

import itertools
from fractions import Fraction

x = np.array([1.0, 2.0, 3.0])
y = np.array([4.0, 5.0, 6.0])
z = np.concatenate([x, y])

# 20 가지를 모두 열거한다.
exact = np.array([z[list(c)].mean() - np.delete(z, list(c)).mean()
                  for c in itertools.combinations(range(6), 3)])
hit = int((np.abs(exact) >= 3 - 1e-12).sum())
print(f"순열 수 = {len(exact)},  극단 배정 = {hit},  "
      f"정확 p = {Fraction(hit, len(exact))} = {hit / len(exact):.4f}")
print(f"SD 열거 = {exact.std():.7f},  닫힌 꼴 = "
      f"{z.std(ddof=1) * np.sqrt(1 / 3 + 1 / 3):.7f}")

for B in (999, 9_999, 99_999):
    obs, p, perms = perm_test_two_sample(x, y, n_perm=B,
                                         rng=np.random.default_rng(0))
    print(f"B = {B:>6}:  obs = {obs:+.1f},  p = {p:.6f},  "
          f"SD = {perms.std(ddof=1):.7f},  p 의 눈금 = 1/{B + 1}")

순열 수 = 20,  극단 배정 = 2,  정확 p = 1/10 = 0.1000
SD 열거 = 1.5275252,  닫힌 꼴 = 1.5275252
B =    999:  obs = -3.0,  p = 0.103000,  SD = 1.5101673,  p 의 눈금 = 1/1000
B =   9999:  obs = -3.0,  p = 0.102000,  SD = 1.5301857,  p 의 눈금 = 1/10000
B =  99999:  obs = -3.0,  p = 0.100650,  SD = 1.5315233,  p 의 눈금 = 1/100000

def perm_test_multi_group(groups, n_perm=9999, rng=None):
    """집단평균의 분산을 통계량으로 쓰는 다집단 순열검정.

    분산분석의 F 대신 집단평균들의 분산을 쓴다. 순열검정에서는 통계량이
    어떤 분포를 따라야 할 까닭이 없으므로, 뜻만 통하면 무엇이든 쓸 수 있다.
    집단 차이가 클수록 이 값이 커지므로 오른쪽 꼬리만 본다.
    """
    rng = rng or np.random.default_rng(0)
    pooled = np.concatenate(groups)
    sizes = [len(g) for g in groups]
    obs_var = np.var([g.mean() for g in groups])

    perm_vars = np.empty(n_perm)
    for i in range(n_perm):
        p = rng.permutation(pooled)
        idx, means = 0, []
        for s in sizes:
            means.append(p[idx:idx + s].mean())
            idx += s
        perm_vars[i] = np.var(means)

    p_value = ((perm_vars >= obs_var).sum() + 1) / (n_perm + 1)
    return obs_var, p_value, perm_vars

import itertools
from fractions import Fraction

groups = [np.array([1.0, 2.0]), np.array([3.0, 4.0]), np.array([5.0, 6.0])]
pooled = np.concatenate(groups)
k, n_g, N = len(groups), 2, len(pooled)

T_obs = np.var([g.mean() for g in groups])
SSB = n_g * sum((g.mean() - pooled.mean()) ** 2 for g in groups)
print(f"T_obs = {T_obs:.7f},   SSB/N = {SSB / N:.7f}")
print(f"E[T] 닫힌 꼴 (k-1)S^2/N = {(k - 1) * pooled.var(ddof=1) / N:.7f}")

# 6! = 720 가지 자리 배정을 모두 열거한다.
Ts = np.array([np.var([a[0:2].mean(), a[2:4].mean(), a[4:6].mean()])
               for a in (pooled[list(p)] for p in itertools.permutations(range(N)))])
hit = int((Ts >= T_obs - 1e-12).sum())
print(f"열거 수 = {len(Ts)},  열거 E[T] = {Ts.mean():.7f},  최댓값 = {Ts.max():.7f}")
print(f"T* >= T_obs 인 배정 = {hit} -> 정확 p = {Fraction(hit, len(Ts))} = {hit / len(Ts):.6f}")

for B in (999, 9999):
    tv, p, pv = perm_test_multi_group(groups, n_perm=B, rng=np.random.default_rng(0))
    print(f"B = {B:>5}: p = {p:.6f},  평균 T* = {pv.mean():.7f}")

T_obs = 2.6666667,   SSB/N = 2.6666667
E[T] 닫힌 꼴 (k-1)S^2/N = 1.1666667
열거 수 = 720,  열거 E[T] = 1.1666667,  최댓값 = 2.6666667
T* >= T_obs 인 배정 = 48 -> 정확 p = 1/15 = 0.066667
B =   999: p = 0.066000,  평균 T* = 1.1281281
B =  9999: p = 0.070300,  평균 T* = 1.1686835

def perm_test_proportion(n_a, conv_a, n_b, conv_b, n_perm=9999, rng=None):
    """두 비율에 대한 순열검정. A/B 검정에 그대로 쓴다.

    전체 전환 수만큼 1 로 채운 배열을 섞으면, 전환이 두 집단에 무작위로
    흩어진 상태가 된다. 표본이 크거나 전환율이 아주 낮아 정규근사가
    미덥지 않을 때 쓸모가 있다.
    """
    rng = rng or np.random.default_rng(0)
    obs_diff = conv_a / n_a - conv_b / n_b
    pooled = np.zeros(n_a + n_b)
    pooled[:conv_a + conv_b] = 1

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

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

from scipy import stats
from fractions import Fraction
from math import comb

n_a, c_a, n_b, c_b = 5, 4, 5, 1
N, K = n_a + n_b, c_a + c_b
step = 1 / n_a + 1 / n_b
obs = c_a / n_a - c_b / n_b

x = np.arange(K + 1)
d = x * step - K / n_b
pmf = stats.hypergeom.pmf(x, N, K, n_a)
sel = np.abs(d) >= abs(obs) - 1e-12
p_frac = sum(Fraction(comb(K, int(v)) * comb(N - K, n_a - int(v)), comb(N, n_a))
             for v in x[sel])
print(f"obs = {obs:.4f},  격자 간격 = {step:.4f},  전체 배정 = C(10,5) = {comb(N, n_a)}")
print(f"걸리는 X = {x[sel]}")
print(f"정확 양측 p = {p_frac} = {float(p_frac):.6f}")
print(f"Fisher 정확검정 = {stats.fisher_exact([[c_a, n_a - c_a], [c_b, n_b - c_b]])[1]:.6f}")
print(f"SD 닫힌 꼴 = {np.sqrt(N / (N - 1) * (K / N) * (1 - K / N) * step):.7f}")
for B in (999, 9999):
    _, p, pv = perm_test_proportion(n_a, c_a, n_b, c_b, n_perm=B,
                                    rng=np.random.default_rng(0))
    print(f"B = {B:>5}: p = {p:.6f},  SD = {pv.std(ddof=1):.7f}")

obs = 0.6000,  격자 간격 = 0.4000,  전체 배정 = C(10,5) = 252
걸리는 X = [0 1 4 5]
정확 양측 p = 13/63 = 0.206349
Fisher 정확검정 = 0.206349
SD 닫힌 꼴 = 0.3333333
B =   999: p = 0.208000,  SD = 0.3358819
B =  9999: p = 0.203000,  SD = 0.3312704

# 자료를 만드는 난수와 순열에 쓰는 난수를 따로 둔다. 이래야 자료를 그대로
# 두고 순열 횟수만 바꿔 보는 식의 실험이 가능하다.
rng_data = np.random.default_rng(11)
rng = np.random.default_rng(3)

def experiment(n_perm, split):
    """자료와 순열에 난수를 나누어 쓸 때와 하나로 쓸 때를 견준다."""
    if split:
        gd, gp = np.random.default_rng(11), np.random.default_rng(3)
    else:
        gd = gp = np.random.default_rng(11)
    a = gd.normal(120, 30, 36)
    b = gd.normal(135, 30, 40)
    perm_test_two_sample(a, b, n_perm=n_perm, rng=gp)     # 순열 난수를 소비한다
    groups = [gd.normal(mu, 25, 30) for mu in [160, 170, 155, 180]]
    return a.mean(), groups[0].mean()

for split, name in [(False, "난수 하나"), (True, "난수 둘  ")]:
    a1, g1 = experiment(999, split)
    a2, g2 = experiment(9999, split)
    print(f"{name}  page_a.mean()    : {a1:.6f} / {a2:.6f}   같은가 {a1 == a2}")
    print(f"{name}  groups[0].mean() : {g1:.6f} / {g2:.6f}   같은가 {g1 == g2}")

난수 하나  page_a.mean()    : 116.800393 / 116.800393   같은가 True
난수 하나  groups[0].mean() : 160.057513 / 158.822504   같은가 False
난수 둘    page_a.mean()    : 116.800393 / 116.800393   같은가 True
난수 둘    groups[0].mean() : 159.293090 / 159.293090   같은가 True

# 두 페이지의 체류시간. 참 평균이 15 만큼 다르다.
page_a = rng_data.normal(120, 30, size=36)
page_b = rng_data.normal(135, 30, size=40)
diff, p, perms = perm_test_two_sample(page_a, page_b, rng=rng)
print(page_a.mean(), page_b.mean(), diff, p)
# 116.80  139.14  -22.34  0.0006

116.80039289004895 139.1365870469421 -22.336194156893157 0.0006

from scipy import stats

z = np.concatenate([page_a, page_b])
m, n = len(page_a), len(page_b)
se_true = 30 * np.sqrt(1 / m + 1 / n)
sd_perm = z.std(ddof=1) * np.sqrt(1 / m + 1 / n)
print(f"참 SE = {se_true:.6f}   (관측 차이 {diff:.6f} 는 참값 -15 에서 "
      f"{(diff + 15) / se_true:+.3f} SE)")
print(f"순열 귀무 SD 닫힌 꼴 = {sd_perm:.6f}   (합친 S = {z.std(ddof=1):.4f})")
print(f"순열 귀무 SD 모의    = {perms.std(ddof=1):.6f}"
      f"   (몬테카를로 오차 {sd_perm / np.sqrt(2 * 9999):.6f})")
print(f"z = {diff / sd_perm:.6f},  정규근사 양측 p = "
      f"{2 * stats.norm.sf(abs(diff / sd_perm)):.6f}")
print(f"순열 p = {p:.6f}  ->  초과 횟수 c = {round(p * 10000) - 1}")
print(f"Welch t 검정 p = {stats.ttest_ind(page_a, page_b, equal_var=False).pvalue:.6f}")
print(f"B = 9999 에서 가능한 최소 p = {1 / 10000}")

참 SE = 6.892024   (관측 차이 -22.336194 는 참값 -15 에서 -1.064 SE)
순열 귀무 SD 닫힌 꼴 = 6.441618   (합친 S = 28.0394)
순열 귀무 SD 모의    = 6.419883   (몬테카를로 오차 0.045551)
z = -3.467482,  정규근사 양측 p = 0.000525
순열 p = 0.000600  ->  초과 횟수 c = 5
Welch t 검정 p = 0.000345
B = 9999 에서 가능한 최소 p = 0.0001

# 네 처치군. 참 평균이 모두 다르다.
groups = [rng_data.normal(mu, 25, 30) for mu in [160, 170, 155, 180]]
var_obs, p_multi, perm_vars = perm_test_multi_group(groups, rng=rng)
print([round(g.mean(), 2) for g in groups], var_obs, p_multi)
# [159.29, 172.22, 160.02, 175.48]  51.759  0.0140

[159.29, 172.22, 160.02, 175.48] 51.75863218978212 0.014

pooled = np.concatenate(groups)
k, n_g, N = len(groups), 30, len(pooled)
print(f"T_obs = {var_obs:.6f},  SSB/N = "
      f"{n_g * sum((g.mean() - pooled.mean()) ** 2 for g in groups) / N:.6f}")
print(f"E[T] 닫힌 꼴 = {(k - 1) * pooled.var(ddof=1) / N:.6f}")
print(f"모의 평균 T* = {perm_vars.mean():.6f}"
      f"   (평균의 몬테카를로 오차 {perm_vars.std(ddof=1) / np.sqrt(9999):.6f})")
print(f"순열 p = {p_multi:.6f}  ->  초과 횟수 c = {round(p_multi * 10000) - 1}")
f, pa = stats.f_oneway(*groups)
print(f"일원분산분석 F = {f:.6f},  p = {pa:.6f}")

T_obs = 51.758632,  SSB/N = 51.758632
E[T] 닫힌 꼴 = 14.884560
모의 평균 T* = 14.808953   (평균의 몬테카를로 오차 0.118600)
순열 p = 0.014000  ->  초과 횟수 c = 139
일원분산분석 F = 3.715378,  p = 0.013538

# 전환율 A/B 검정. 표본이 2만 이상인데 전환은 200 안팎이라 전환율이
# 1% 아래다. 이런 자료에서 순열검정이 쓸모 있다.
diff_ab, p_ab, perms_ab = perm_test_proportion(23739, 200, 22588, 182, rng=rng)
print(diff_ab)      # 0.000368

0.0003675791182059275

N, K = 23739 + 22588, 200 + 182
step = 1 / 23739 + 1 / 22588
x = np.arange(K + 1)
d = x * step - K / 22588
pmf = stats.hypergeom.pmf(x, N, K, 23739)
p_exact = pmf[np.abs(d) >= abs(diff_ab) - 1e-18].sum()
sd_closed = np.sqrt(N / (N - 1) * (K / N) * (1 - K / N) * step)

print(f"순열 p (B=9999) = {p_ab:.6f}")
print(f"정확 순열 p     = {p_exact:.6f}   (E[p-hat] = {(9999 * p_exact + 1) / 10000:.6f},"
      f"  SD = {np.sqrt(9999 * p_exact * (1 - p_exact)) / 10000:.6f})")
print(f"Fisher 정확검정 = "
      f"{stats.fisher_exact([[200, 23739 - 200], [182, 22588 - 182]])[1]:.6f}")
print(f"귀무 SD 닫힌 꼴 = {sd_closed:.8f},  모의 = {perms_ab.std(ddof=1):.8f}")
print(f"섞은 원소의 수 = {9999 * N:,}")

순열 p (B=9999) = 0.681900
정확 순열 p     = 0.681128   (E[p-hat] = 0.681160,  SD = 0.004660)
Fisher 정확검정 = 0.681128
귀무 SD 닫힌 꼴 = 0.00084056,  모의 = 0.00084891
섞은 원소의 수 = 463,223,673

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

pvals = []
for _ in range(1000):
    x = rng.normal(0, 1, 50)
    y = rng.normal(0, 1, 50)
    _, p, _ = perm_test_two_sample(x, y, n_perm=999, rng=rng)
    pvals.append(p)

plt.hist(pvals, bins=20, edgecolor='k')
plt.xlabel('p-value'); plt.ylabel('Frequency')
plt.title('Distribution of p-values under H0')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
rng = np.random.default_rng(1)

diff_ab, p_ab, perms_ab = perm_test_proportion(23739, 200, 22588, 182,
                                               n_perm=5000, rng=rng)
mu, sigma = perms_ab.mean(), perms_ab.std()

fig, ax = plt.subplots(figsize=(8, 4))
ax.hist(perms_ab * 100, bins=40, density=True, alpha=0.7, edgecolor='k')
x = np.linspace(perms_ab.min() * 100, perms_ab.max() * 100, 200)
ax.plot(x, stats.norm.pdf(x, mu * 100, sigma * 100), 'r-', lw=2,
        label='Normal approximation')
ax.set_xlabel('Difference in rate (%)')
ax.set_ylabel('Density')
ax.legend()
plt.tight_layout()
plt.show()

import numpy as np, itertools
from scipy import stats

def paired_perm_exact(before, after):
    d = np.asarray(after) - np.asarray(before)
    n = len(d)
    S = np.array(list(itertools.product([1, -1], repeat=n)))
    m = (S * d).mean(axis=1)
    count = (np.abs(m) >= abs(d.mean()) - 1e-12).sum()
    return d.mean(), count, len(m), count / len(m)

before = [82, 78, 91, 85, 73]
after  = [88, 82, 95, 89, 78]
print(paired_perm_exact(before, after))
# (4.6, 2, 32, 0.0625)
print(stats.ttest_rel(after, before))       # t = 11.5,  p = 0.00033
print(stats.wilcoxon(np.array(after) - np.array(before), method='exact'))
# W = 0.0,  p = 0.0625

(4.6, 2, 32, 0.0625)
TtestResult(statistic=11.5, pvalue=0.00032642125636699325, df=4)
WilcoxonResult(statistic=0.0, pvalue=0.0625)