이표본 비모수 검정¶

import numpy as np
from scipy.stats import mannwhitneyu, rankdata

# 두 교수법의 시험 점수를 견준다.
group_a = np.array([85, 78, 92, 88, 76, 95, 89, 82, 91, 87])
group_b = np.array([72, 68, 81, 75, 70, 77, 74, 69, 73, 71])

# Mann-Whitney U 는 이표본 t 검정의 비모수 대응이다. 귀무가설이 "두 평균이
# 같다"가 아니라 "무작위로 뽑은 A 가 B 보다 클 확률이 1/2"이라는 점에 주의한다.
stat, p_value = mannwhitneyu(group_a, group_b, alternative='two-sided')

print(f"U statistic: {stat}")        # 97.0
print(f"P-value: {p_value:.6f}")     # 0.000440  (asymptotic)

# 표본이 작고 동점이 없으면 정확분포를 쓸 수 있다. 근사값과 꽤 갈린다.
print(mannwhitneyu(group_a, group_b, method='exact').pvalue)   # 7.58e-05

# 유의하다는 결론이 났을 때 어느 쪽이 큰지는 평균순위로 말한다.
# 동점에 중간순위를 주려면 argsort 가 아니라 rankdata 를 써야 한다.
combined = np.concatenate([group_a, group_b])
ranks = rankdata(combined)
mean_rank_a = ranks[:len(group_a)].mean()
mean_rank_b = ranks[len(group_a):].mean()
print(f"Mean rank Group A: {mean_rank_a:.1f}")   # 15.2
print(f"Mean rank Group B: {mean_rank_b:.1f}")   # 5.8

alpha = 0.05
if p_value < alpha:
    print("Reject H0: Significant difference between groups.")
    if mean_rank_a > mean_rank_b:
        print("Group A tends to have larger values.")
    else:
        print("Group B tends to have larger values.")
else:
    print("Fail to reject H0.")

U statistic: 97.0
P-value: 0.000440
7.577561757128321e-05
Mean rank Group A: 15.2
Mean rank Group B: 5.8
Reject H0: Significant difference between groups.
Group A tends to have larger values.

import itertools
import math
from scipy import stats

n1, n2 = len(group_a), len(group_b)
N = n1 + n2
R1, R2 = ranks[:n1].sum(), ranks[n1:].sum()
U1 = R1 - n1 * (n1 + 1) / 2
U2 = R2 - n2 * (n2 + 1) / 2
print(f"R1 = {R1},  R2 = {R2},  합 = {R1 + R2}  (N(N+1)/2 = {N * (N + 1) / 2})")
print(f"U1 = {U1},  U2 = {U2},  합 = {U1 + U2}  (n1*n2 = {n1 * n2})")
print(f"min(U1, U2) = {min(U1, U2)}   <- scipy 가 돌려준 것은 U1 = {stat}")
print(f"P(A > B) 추정 = U1/(n1 n2) = {U1 / (n1 * n2)}")

# 귀무분포를 전수 열거한다. C(20, 10) = 184756 가지뿐이다.
total = math.comb(N, n1)
Us = np.array([sum(c) for c in itertools.combinations(range(1, N + 1), n1)]) \
     - n1 * (n1 + 1) / 2
E = n1 * n2 / 2
V = n1 * n2 * (N + 1) / 12
hit = int((np.abs(Us - E) >= abs(U1 - E)).sum())
print(f"\n열거한 가지 수 {total},  평균 {Us.mean()},  분산 {Us.var()}"
      f"  (공식: {E}, {V})")
print(f"|U - {E:.0f}| >= {abs(U1 - E):.0f} 인 것 {hit} 개"
      f"  ->  정확 p = {hit}/{total} = {hit / total:.6e}")
print(f"scipy method='exact' : {mannwhitneyu(group_a, group_b, method='exact').pvalue:.6e}")

z = (abs(U1 - E) - 0.5) / np.sqrt(V)
print(f"\nsigma = {np.sqrt(V):.6f},  z = {z:.6f},"
      f"  정규근사 p = {2 * stats.norm.sf(z):.6e}")
print(f"근사/정확 = {2 * stats.norm.sf(z) / (hit / total):.2f} 배,"
      f"  절대차 = {2 * stats.norm.sf(z) - hit / total:.6f}")

R1 = 152.0,  R2 = 58.0,  합 = 210.0  (N(N+1)/2 = 210.0)
U1 = 97.0,  U2 = 3.0,  합 = 100.0  (n1*n2 = 100)
min(U1, U2) = 3.0   <- scipy 가 돌려준 것은 U1 = 97.0
P(A > B) 추정 = U1/(n1 n2) = 0.97

열거한 가지 수 184756,  평균 50.0,  분산 175.0  (공식: 50.0, 175.0)
|U - 50| >= 47 인 것 14 개  ->  정확 p = 14/184756 = 7.577562e-05
scipy method='exact' : 7.577562e-05

sigma = 13.228757,  z = 3.515070,  정규근사 p = 4.396388e-04
근사/정확 = 5.80 배,  절대차 = 0.000364

import numpy as np
from scipy import stats

a = np.array([50, 55, 60, 65, 70])
b = np.array([45, 50, 55, 60, 65])
allv = np.concatenate([a, b])
gm = np.median(allv)
print(f"정렬 {np.sort(allv).tolist()}  ->  전체 중앙값 {gm}")
print(f"= {gm} 인 관측값 {int((allv == gm).sum())} 개")
print(f"A: 위 {int((a > gm).sum())}, 아래 {int((a < gm).sum())}  |  "
      f"B: 위 {int((b > gm).sum())}, 아래 {int((b < gm).sum())}")

print("\n가능한 분할표를 모두 열거한다")
for k in range(6):
    T = np.array([[k, 5 - k], [5 - k, k]])
    c0 = stats.chi2_contingency(T, correction=False)
    c1 = stats.chi2_contingency(T, correction=True)
    print(f"  a={k}  {T.tolist()}   chi2 = {c0.statistic:6.2f}, p = {c0.pvalue:.6f}"
          f"   |  Yates chi2 = {c1.statistic:5.2f}, p = {c1.pvalue:.6f}")
print(f"\nscipy median_test: {stats.median_test(a, b)}")

정렬 [45, 50, 50, 55, 55, 60, 60, 65, 65, 70]  ->  전체 중앙값 57.5
= 57.5 인 관측값 0 개
A: 위 3, 아래 2  |  B: 위 2, 아래 3

가능한 분할표를 모두 열거한다
  a=0  [[0, 5], [5, 0]]   chi2 =  10.00, p = 0.001565   |  Yates chi2 =  6.40, p = 0.011412
  a=1  [[1, 4], [4, 1]]   chi2 =   3.60, p = 0.057780   |  Yates chi2 =  1.60, p = 0.205903
  a=2  [[2, 3], [3, 2]]   chi2 =   0.40, p = 0.527089   |  Yates chi2 =  0.00, p = 1.000000
  a=3  [[3, 2], [2, 3]]   chi2 =   0.40, p = 0.527089   |  Yates chi2 =  0.00, p = 1.000000
  a=4  [[4, 1], [1, 4]]   chi2 =   3.60, p = 0.057780   |  Yates chi2 =  1.60, p = 0.205903
  a=5  [[5, 0], [0, 5]]   chi2 =  10.00, p = 0.001565   |  Yates chi2 =  6.40, p = 0.011412

scipy median_test: MedianTestResult(statistic=0.0, pvalue=1.0, median=57.5, table=array([[3, 2],
       [2, 3]]))

import numpy as np
from scipy.stats import chi2_contingency

def moods_median_test(*groups, correction=True):
    """Mood 의 중앙값 검정. 여러 집단의 중앙값을 견준다.

    자료 전체의 중앙값을 구한 뒤, 집단마다 그보다 큰 값과 작은 값의 개수를
    세어 분할표를 만들고 카이제곱 검정을 돌린다. 값을 "위/아래" 둘로만
    나누므로 순위검정보다도 정보를 더 버리고, 그만큼 검정력이 낮다.
    대신 가정이 거의 없어 아주 거친 자료에도 쓸 수 있다.

    매개변수
    --------
    groups : 집단을 나타내는 배열 여럿
    correction : Yates 연속성 보정 여부 (2x2 표에서만 뜻이 있다)

    돌려주는 값
    ----------
    chi2_stat, p_value, contingency_table
    """
    # 전체 중앙값을 기준선으로 삼는다.
    combined_data = np.concatenate(groups)
    overall_median = np.median(combined_data)

    contingency_table = []
    for group in groups:
        # 중앙값과 정확히 같은 값은 어느 쪽에도 넣지 않는다.
        above_median = np.sum(group > overall_median)
        below_median = np.sum(group < overall_median)
        contingency_table.append([above_median, below_median])

    contingency_table = np.array(contingency_table).T
    chi2_stat, p_value, _, _ = chi2_contingency(contingency_table,
                                                correction=correction)

    return chi2_stat, p_value, contingency_table

# 예시 자료
group_a = np.array([50, 55, 60, 65, 70])
group_b = np.array([45, 50, 55, 60, 65])

# 보정 여부에 따라 p-값이 꽤 갈린다. 2x2 표에서 표본이 작을 때 그렇다.
for corr in (True, False):
    chi2_stat, p_value, table = moods_median_test(group_a, group_b,
                                                  correction=corr)
    print(f"correction={corr}: chi2={chi2_stat:.4f}, p={p_value:.4f}")
print(f"Contingency Table:\n{table}")

correction=True: chi2=0.0000, p=1.0000
correction=False: chi2=0.4000, p=0.5271
Contingency Table:
[[3 2]
 [2 3]]

import math
from scipy.stats import fisher_exact

T = np.array([[3, 2], [2, 3]])
E = np.outer(T.sum(1), T.sum(0)) / T.sum()
print(f"기대도수 {E.tolist()},  |O - E| = {np.abs(T - E).tolist()}")
print(f"예이츠 분자 (|O-E| - 0.5)^2 = {((np.abs(T - E) - 0.5) ** 2).tolist()}")

print(f"\n관측된 표 a=3")
print(f"  보정 없음 p = {chi2_contingency(T, correction=False).pvalue:.6f}")
print(f"  예이츠    p = {chi2_contingency(T, correction=True).pvalue:.6f}")
print(f"  피셔      p = {fisher_exact(T).pvalue:.6f}")

T5 = np.array([[5, 0], [0, 5]])
print(f"가장 치우친 표 a=5")
print(f"  보정 없음 p = {chi2_contingency(T5, correction=False).pvalue:.6f}")
print(f"  예이츠    p = {chi2_contingency(T5, correction=True).pvalue:.6f}")
print(f"  피셔      p = {fisher_exact(T5).pvalue:.6f}")

# H0 아래에서 a 의 분포(초기하)와 alpha=0.05 검정의 실제 수준
probs = np.array([math.comb(5, k) * math.comb(5, 5 - k) / math.comb(10, 5)
                  for k in range(6)])
print(f"\nP(a = k) = {[f'{p:.6f}' for p in probs]}"
      f"   분자 {[math.comb(5, k) * math.comb(5, 5 - k) for k in range(6)]} / {math.comb(10, 5)}")
for name, pv in (("보정 없음", [chi2_contingency(np.array([[k, 5 - k], [5 - k, k]]),
                                              correction=False).pvalue for k in range(6)]),
                 ("예이츠   ", [chi2_contingency(np.array([[k, 5 - k], [5 - k, k]]),
                                              correction=True).pvalue for k in range(6)]),
                 ("피셔     ", [fisher_exact(np.array([[k, 5 - k], [5 - k, k]])).pvalue
                               for k in range(6)])):
    rej = [k for k in range(6) if pv[k] <= 0.05]
    print(f"  {name}: 기각역 a in {rej},  실제 수준 = {probs[rej].sum():.6f}")

기대도수 [[2.5, 2.5], [2.5, 2.5]],  |O - E| = [[0.5, 0.5], [0.5, 0.5]]
예이츠 분자 (|O-E| - 0.5)^2 = [[0.0, 0.0], [0.0, 0.0]]

관측된 표 a=3
  보정 없음 p = 0.527089
  예이츠    p = 1.000000
  피셔      p = 1.000000
가장 치우친 표 a=5
  보정 없음 p = 0.001565
  예이츠    p = 0.011412
  피셔      p = 0.007937

P(a = k) = ['0.003968', '0.099206', '0.396825', '0.396825', '0.099206', '0.003968']   분자 [1, 25, 100, 100, 25, 1] / 252
  보정 없음: 기각역 a in [0, 5],  실제 수준 = 0.007937
  예이츠   : 기각역 a in [0, 5],  실제 수준 = 0.007937
  피셔     : 기각역 a in [0, 5],  실제 수준 = 0.007937

from scipy.stats import kruskal

group_a = [85, 78, 92, 88, 76]
group_b = [72, 68, 81, 75, 70]
group_c = [90, 95, 88, 92, 87]

# Kruskal-Wallis 는 일원배치 분산분석의 비모수 대응이다. 집단이 셋 이상일
# 때 쓰며, Mann-Whitney U 를 여러 집단으로 넓힌 것이라고 보면 된다.
stat, p_value = kruskal(group_a, group_b, group_c)

print(f"H statistic: {stat:.4f}")   # 9.4136
print(f"P-value: {p_value:.4f}")    # 0.0090

alpha = 0.05
if p_value < alpha:
    print("Reject H0: At least one group differs significantly.")
else:
    print("Fail to reject H0.")

H statistic: 9.4136
P-value: 0.0090
Reject H0: At least one group differs significantly.

import itertools
import math
from collections import Counter

import numpy as np
from scipy import stats

allv = np.array(group_a + group_b + group_c)
N = len(allv)
r = stats.rankdata(allv)
R = [r[0:5].sum(), r[5:10].sum(), r[10:15].sum()]
print(f"순위  {r.tolist()}")
print(f"동점  {[(k, v) for k, v in sorted(Counter(allv.tolist()).items()) if v > 1]}")
print(f"R = {R},  합 = {sum(R)}  (N(N+1)/2 = {N * (N + 1) / 2})")
print(f"평균순위 = {[x / 5 for x in R]}   (전체 {(N + 1) / 2})")

H = 12 / (N * (N + 1)) * sum(x ** 2 / 5 for x in R) - 3 * (N + 1)
tie = sum(t ** 3 - t for t in Counter(allv.tolist()).values())
C = 1 - tie / (N ** 3 - N)
print(f"\n보정 없는 H = {H:.6f}")
print(f"sum(t^3-t) = {tie},  C = 1 - {tie}/{N ** 3 - N} = {C:.6f}")
print(f"H/C = {H / C:.6f}   (scipy {stat:.6f})")
print(f"chi2_2 근사 p = {stats.chi2.sf(H / C, 2):.6f}   (scipy {p_value:.6f})")

# 귀무분포를 전수 열거한다. C(15,5) * C(10,5) = 756756 가지.
sub2 = np.array(list(itertools.combinations(range(10), 5)))
const = 12 / (N * (N + 1)) / 5
total_rank = r.sum()
chunks = []
for first in itertools.combinations(range(N), 5):
    s1 = r[list(first)].sum()
    rest = np.array([r[i] for i in range(N) if i not in first])
    s2 = rest[sub2].sum(axis=1)
    s3 = total_rank - s1 - s2
    chunks.append(const * (s1 ** 2 + s2 ** 2 + s3 ** 2) - 3 * (N + 1))
Hs = np.concatenate(chunks)
hit = int((Hs >= H - 1e-9).sum())
print(f"\n열거한 가지 수 {len(Hs)}  (= {math.comb(15, 5)} x {math.comb(10, 5)})")
print(f"열거한 H 의 평균 {Hs.mean():.6f}   (chi2_2 의 평균 = 2)")
print(f"H >= {H:.2f} 인 것 {hit} 개  ->  정확 p = {hit}/{len(Hs)} = {hit / len(Hs):.7f}")
print(f"chi2_2 근사 / 정확 = {stats.chi2.sf(H / C, 2) / (hit / len(Hs)):.2f} 배")

순위  [8.0, 6.0, 13.5, 10.5, 5.0, 3.0, 1.0, 7.0, 4.0, 2.0, 12.0, 15.0, 10.5, 13.5, 9.0]
동점  [(88, 2), (92, 2)]
R = [43.0, 17.0, 60.0],  합 = 120.0  (N(N+1)/2 = 120.0)
평균순위 = [8.6, 3.4, 12.0]   (전체 8.0)

보정 없는 H = 9.380000
sum(t^3-t) = 12,  C = 1 - 12/3360 = 0.996429
H/C = 9.413620   (scipy 9.413620)
chi2_2 근사 p = 0.009034   (scipy 0.009034)

열거한 가지 수 756756  (= 3003 x 252)
열거한 H 의 평균 1.992857   (chi2_2 의 평균 = 2)
H >= 9.38 인 것 1890 개  ->  정확 p = 1890/756756 = 0.0024975
chi2_2 근사 / 정확 = 3.62 배

import numpy as np
from scipy.stats import mannwhitneyu, median_test
a = np.array([50, 55, 60, 65, 70]); b = np.array([45, 50, 55, 60, 65])
print(mannwhitneyu(a, b, method='exact'))    # U=17.0, p=0.4206
print(median_test(a, b))                     # stat=0.0, p=1.0

MannwhitneyuResult(statistic=17.0, pvalue=0.42063492063492064)
MedianTestResult(statistic=0.0, pvalue=1.0, median=57.5, table=array([[3, 2],
       [2, 3]]))

import numpy as np
from scipy.stats import median_test
a = np.array([3, 3, 3, 4, 5, 5])
b = np.array([1, 2, 3, 3, 3, 3])
print(np.median(np.concatenate([a, b])))    # 3.0

for ties in ("below", "above", "ignore"):
    r = median_test(a, b, ties=ties)
    print(ties, r.statistic.round(4), r.pvalue.round(4))
    print(r.table)

3.0
below 1.7778 0.1824
[[3 0]
 [3 6]]
above 0.6 0.4386
[[6 4]
 [0 2]]
ignore 1.7014 0.1921
[[3 0]
 [0 2]]

import numpy as np
from scipy.stats import kruskal, mannwhitneyu, rankdata
rng = np.random.default_rng(0)
x = rng.normal(0, 1, 15); y = rng.normal(0.7, 1, 18)

H, pH = kruskal(x, y)
r = mannwhitneyu(x, y, method='asymptotic', use_continuity=False)
n1, n2 = len(x), len(y)
mu = n1 * n2 / 2
sd = np.sqrt(n1 * n2 * (n1 + n2 + 1) / 12)
Z = (r.statistic - mu) / sd
print(H, Z**2, H - Z**2)        # 6.4052  6.4052  ~0
print(pH, r.pvalue)             # 0.01138  0.01138

6.405228758169926 6.405228758169934 -7.993605777301127e-15
0.011378476531193307 0.01137847653119324

import numpy as np, itertools
from scipy.stats import mannwhitneyu
ga = [85, 78, 92, 88, 76, 95, 89, 82, 91, 87]
gb = [72, 68, 81, 75, 70, 77, 74, 69, 73, 71]
print(mannwhitneyu(ga, gb, method='exact').pvalue)        # 7.5776e-05
print(mannwhitneyu(ga, gb, method='asymptotic',
                   use_continuity=True).pvalue)           # 4.3964e-04
print(mannwhitneyu(ga, gb, method='asymptotic',
                   use_continuity=False).pvalue)          # 3.8106e-04

7.577561757128321e-05
0.00043963875262656465
0.00038105845205068555