Wilcoxon 검정¶

UserWarning: Exact p-value calculation does not work if there are zeros.
Switching to normal approximation.

import numpy as np
from scipy import stats

# 같은 사람을 두 번 잰 대응자료다. 짝지어져 있다는 것이 핵심이라,
# 두 열을 따로 떼어 독립표본 검정을 쓰면 안 된다.
paired_data = np.array([
    [93, 76], [70, 72], [81, 75], [65, 68], [79, 65],
    [54, 54], [94, 88], [91, 81], [77, 65], [65, 57],
    [95, 86], [89, 87], [78, 78], [80, 77], [76, 76]
])

# zero_method="pratt" 는 차이가 0 인 쌍을 순위 매기기에는 넣되 통계량에서는
# 뺀다. "wilcox"(기본값)는 아예 버린다. 동점이 여럿일 때 결과가 갈린다.
statistic, p_value = stats.wilcoxon(
    paired_data[:, 0], paired_data[:, 1],
    alternative="two-sided",
    method="approx",       # 옛 SciPy의 mode= 는 제거되었다
    zero_method="pratt"
)
print(f"W = {statistic}, p = {p_value:.4f}")
# W = 11.0, p = 0.0086

W = 11.0, p = 0.0086

import itertools
import warnings
from collections import Counter

d = paired_data[:, 0] - paired_data[:, 1]
with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    r_exact = stats.wilcoxon(d, method="exact", zero_method="wilcox")
print(f"method='exact'  -> W = {r_exact.statistic}, p = {r_exact.pvalue:.12f}")
print(f"  경고: {str(caught[0].message)}")
print(f"method='approx' -> p = "
      f"{stats.wilcoxon(d, method='approx', zero_method='wilcox').pvalue:.12f}"
      f"   <- 같은 값이다")

# 동점만 있고 0 은 없는 자료: 경고도 없이 틀린 정확값을 준다
tie_only = np.array([1, 1, 2, 3, -1, 5])
rk = stats.rankdata(np.abs(tie_only))
S = np.array(list(itertools.product([0, 1], repeat=len(tie_only))))
Wp = S @ rk
obs = min(rk[tie_only > 0].sum(), rk[tie_only < 0].sum())
enum_p = (np.minimum(Wp, rk.sum() - Wp) <= obs + 1e-9).mean()
print(f"\n동점만 있는 자료 {tie_only.tolist()}  중간순위 {rk.tolist()}")
print(f"  scipy method='exact' = "
      f"{stats.wilcoxon(tie_only, method='exact').pvalue:.6f}  (경고 없음)")
print(f"  중간순위로 전수 열거 = {enum_p:.6f}   <- 참값")

method='exact'  -> W = 5.0, p = 0.007579201614
  경고: Exact p-value calculation does not work if there are zeros. Switching to normal approximation.
method='approx' -> p = 0.007579201614   <- 같은 값이다

동점만 있는 자료 [1, 1, 2, 3, -1, 5]  중간순위 [2.0, 2.0, 4.0, 5.0, 2.0, 6.0]
  scipy method='exact' = 0.093750  (경고 없음)
  중간순위로 전수 열거 = 0.125000   <- 참값

nonzero = d[d != 0]
r = stats.rankdata(np.abs(nonzero))
n_eff = len(nonzero)
total = r.sum()
print(f"n' = {n_eff},  중간순위 {r.tolist()}")
print(f"순위합 {total}  (n'(n'+1)/2 = {n_eff * (n_eff + 1) / 2})")

# 2^12 = 4096 가지 부호 배정을 전수 열거한다.
signs = np.array(list(itertools.product([0, 1], repeat=n_eff)))
W_plus = signs @ r
W_min = np.minimum(W_plus, total - W_plus)
W_obs = r[nonzero < 0].sum()
hit = int((W_min <= W_obs + 1e-9).sum())
print(f"\n열거한 가지 수 {len(W_plus)}  (2^{n_eff} = {2 ** n_eff})")
print(f"관측 W- = {W_obs},  min(W+, W-) <= {W_obs} 인 배정 {hit} 개")
print(f"정확   양측 p = {hit}/{len(W_plus)} = {hit / len(W_plus):.6f}")
approx = stats.wilcoxon(d, method="approx", zero_method="wilcox").pvalue
print(f"정규근사 양측 p = {approx:.6f}   ->  근사/정확 = "
      f"{approx / (hit / len(W_plus)):.2f} 배")

tie = sum(c ** 3 - c for c in Counter(r).values())
V = (n_eff * (n_eff + 1) * (2 * n_eff + 1) - tie / 2) / 24
print(f"\n열거한 E[W+] = {W_plus.mean()},  공식 n'(n'+1)/4 = {n_eff * (n_eff + 1) / 4}")
print(f"열거한 Var    = {W_plus.var()},  동점 보정 공식 = {V}"
      f"   (sum(t^3-t) = {tie})")

n' = 12,  중간순위 [12.0, 1.5, 5.5, 3.5, 11.0, 5.5, 9.0, 10.0, 7.0, 8.0, 1.5, 3.5]
순위합 78.0  (n'(n'+1)/2 = 78.0)

열거한 가지 수 4096  (2^12 = 4096)
관측 W- = 5.0,  min(W+, W-) <= 5.0 인 배정 20 개
정확   양측 p = 20/4096 = 0.004883
정규근사 양측 p = 0.007579   ->  근사/정확 = 1.55 배

열거한 E[W+] = 39.0,  공식 n'(n'+1)/4 = 39.0
열거한 Var    = 162.125,  동점 보정 공식 = 162.125   (sum(t^3-t) = 18)

from scipy import stats

# 서로 다른 두 집단에서 독립적으로 얻은 관측값이어야 한다.
# 크기가 달라도 되는 것이 대응표본 검정과의 큰 차이다.
a = [12, 15, 18, 22, 25]
b = [8, 10, 14, 19, 21, 24]

# ranksums 는 순위합 검정이다. Mann-Whitney U 와 수학적으로 같은 검정이지만
# 정규근사만 쓰고 동점 보정을 하지 않는다. 표본이 작으면 mannwhitneyu 쪽이 낫다.
statistic, p_value = stats.ranksums(a, b, alternative="two-sided")
print(f"Z = {statistic:.4f}, p = {p_value:.4f}")
# Z = 0.7303, p = 0.4652

Z = 0.7303, p = 0.4652

import numpy as np

rng = np.random.default_rng(0)
B, n = 2000, 15

print("sigma_s  delta   부호순위   순위합")
for sigma_s, delta in ((10.0, 1.0), (1.0, 1.0), (0.0, 1.0), (10.0, 0.0)):
    p_paired, p_indep = [], []
    for _ in range(B):
        S = rng.normal(0, sigma_s, n)           # 개체 간 변동
        X = S + delta + rng.normal(0, 1, n)
        Y = S + rng.normal(0, 1, n)
        p_paired.append(stats.wilcoxon(X, Y).pvalue)
        p_indep.append(stats.ranksums(X, Y).pvalue)
    print(f"{sigma_s:>7} {delta:>6}   "
          f"{np.mean(np.array(p_paired) < 0.05):>8.3f} "
          f"{np.mean(np.array(p_indep) < 0.05):>8.3f}")
print(f"\n반복 {B} 회의 몬테카를로 오차: 0.7 근처에서 "
      f"{np.sqrt(0.7 * 0.3 / B):.4f}, 0.05 근처에서 {np.sqrt(0.05 * 0.95 / B):.4f}")

sigma_s  delta   부호순위   순위합
   10.0    1.0      0.690    0.000
    1.0    1.0      0.701    0.439
    0.0    1.0      0.700    0.742
   10.0    0.0      0.049    0.000

반복 2000 회의 몬테카를로 오차: 0.7 근처에서 0.0102, 0.05 근처에서 0.0049

import numpy as np
from scipy import stats
D = np.array([4, -1, 7, 3, -2, 5])
print(stats.wilcoxon(D, method='exact'))
# WilcoxonResult(statistic=3.0, pvalue=0.15625)
print(stats.wilcoxon(D, alternative='greater', method='exact').pvalue)
# 0.078125

WilcoxonResult(statistic=3.0, pvalue=0.15625)
0.078125

from scipy import stats

a = [12, 15, 18, 22, 25]
b = [8, 10, 14, 19, 21, 24]

print(stats.ranksums(a, b))
# RanksumsResult(statistic=0.7303, pvalue=0.4652)
print(stats.mannwhitneyu(a, b, method='exact'))
# MannwhitneyuResult(statistic=19.0, pvalue=0.5368)

RanksumsResult(statistic=0.7302967433402214, pvalue=0.4652088184521418)
MannwhitneyuResult(statistic=19.0, pvalue=0.5367965367965368)

import numpy as np
from scipy import stats
D = np.array([4, -1, 7, 3, -2, 5])
r = stats.rankdata(np.abs(D))
print("W+ =", r[D > 0].sum(), " W- =", r[D < 0].sum())   # 18.0  3.0

for alt in ("two-sided", "greater", "less"):
    res = stats.wilcoxon(D, alternative=alt, method='exact')
    print(alt, res.statistic, round(res.pvalue, 5))

W+ = 18.0  W- = 3.0
two-sided 3.0 0.15625
greater 18.0 0.07812
less 18.0 0.95312