Kruskal-Wallis 검정¶

from scipy import stats
A = [12, 15, 14, 10, 13]; B = [20, 18, 22, 17, 19]; C = [8, 11, 9, 7, 10]
print(stats.kruskal(A, B, C))
# KruskalResult(statistic=11.8161, pvalue=0.0027175)

KruskalResult(statistic=11.816100178890885, pvalue=0.0027174805686494947)

import itertools
from collections import Counter
from fractions import Fraction
from math import comb

import numpy as np

allv = np.array(A + B + C)
N = len(allv)
r = stats.rankdata(allv)
tie = sum(t ** 3 - t for t in Counter(allv.tolist()).values())
corr = 1 - tie / (N ** 3 - N)
R = [r[0:5].sum(), r[5:10].sum(), r[10:15].sum()]
H_raw = 12 / (N * (N + 1)) * sum(x ** 2 / 5 for x in R) - 3 * (N + 1)
print(f"순위 {r.tolist()}")
print(f"R = {R},  보정 전 H = {H_raw:.6f},  보정계수 = {corr:.6f},"
      f"  H/C = {H_raw / corr:.6f}")

# 15 개 순위를 5 개씩 세 집단에 나누는 모든 방법을 열거한다.
sub2 = np.array(list(itertools.combinations(range(10), 5)))
const = 12 / (N * (N + 1)) / 5
total = 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 - s1 - s2
    chunks.append(const * (s1 ** 2 + s2 ** 2 + s3 ** 2) - 3 * (N + 1))
Hs = np.concatenate(chunks)
print(f"\n열거한 가지 수 {len(Hs)}  (= {comb(15, 5)} x {comb(10, 5)})")
print(f"열거한 H 의 평균 {Hs.mean():.6f}  (chi2_2 의 평균 2),  최댓값 {Hs.max()}")

print("\n명목     임계값   실제 수준   실제/명목")
for alpha in (0.10, 0.05, 0.01, 0.001):
    c = stats.chi2.isf(alpha, 2)
    size = float((Hs / corr >= c - 1e-9).mean())
    print(f"{alpha:>6}  {c:>8.4f}  {size:>10.6f}  {size / alpha:>9.3f}")

hit = int((Hs >= H_raw - 1e-9).sum())
print(f"\nH >= {H_raw:.3f} 인 배치 {hit} 개")
print(f"정확     p = {hit}/{len(Hs)} = {Fraction(hit, len(Hs))} = {hit / len(Hs):.7f}")
print(f"chi2_2   p = {stats.chi2.sf(H_raw / corr, 2):.7f}"
      f"   ->  {stats.chi2.sf(H_raw / corr, 2) / (hit / len(Hs)):.1f} 배")

순위 [7.0, 10.0, 9.0, 4.5, 8.0, 14.0, 12.0, 15.0, 11.0, 13.0, 2.0, 6.0, 3.0, 1.0, 4.5]
R = [38.5, 65.0, 16.5],  보정 전 H = 11.795000,  보정계수 = 0.998214,  H/C = 11.816100

열거한 가지 수 756756  (= 3003 x 252)
열거한 H 의 평균 1.996429  (chi2_2 의 평균 2),  최댓값 12.5

명목     임계값   실제 수준   실제/명목
   0.1    4.6052    0.095508      0.955
  0.05    5.9915    0.042783      0.856
  0.01    9.2103    0.003235      0.323
 0.001   13.8155    0.000000      0.000

H >= 11.795 인 배치 24 개
정확     p = 24/756756 = 2/63063 = 0.0000317
chi2_2   p = 0.0027175   ->  85.7 배

from scipy import stats
a = [15, 18, 20, 17]; b = [22, 25, 19, 23]; c = [12, 14, 16, 13]
print(stats.kruskal(a, b, c))          # H=8.6538, p=0.013208
print(stats.f_oneway(a, b, c))         # F=16.130, p=0.001057

KruskalResult(statistic=8.653846153846153, pvalue=0.013208125416383262)
F_onewayResult(statistic=16.129629629629616, pvalue=0.0010574135527856804)

import numpy as np, itertools
from scipy import stats

N, k, n = 12, 3, 4
ranks = np.arange(1, N + 1)
Rbar = (N + 1) / 2
Hs = []
for g1 in itertools.combinations(range(N), n):
    rest = [i for i in range(N) if i not in g1]
    for g2 in itertools.combinations(rest, n):
        g3 = [i for i in rest if i not in g2]
        means = [ranks[list(g)].mean() for g in (g1, g2, g3)]
        Hs.append(12 / (N * (N + 1)) * n * sum((m - Rbar) ** 2 for m in means))
Hs = np.array(Hs)
print(len(Hs))                                  # 34650
print("exact p:", (Hs >= 8.6538 - 1e-9).mean()) # 0.00139
print("chi2  p:", stats.chi2.sf(8.6538, 2))     # 0.01321

34650
exact p: 0.0013852813852813853
chi2  p: 0.013208430222794458

import numpy as np
from scipy import stats
rng = np.random.default_rng(0)
B = 2000

def sim(n):
    rej = 0
    for _ in range(B):
        g1 = rng.exponential(1, n) - np.log(2)          # 오른쪽 치우침
        g2 = rng.normal(0, 1, n)                        # 대칭
        g3 = -(rng.exponential(1, n) - np.log(2))       # 왼쪽 치우침
        rej += stats.kruskal(g1, g2, g3).pvalue < 0.05
    return rej / B

for n in (10, 20, 50, 100):
    print(n, sim(n))

10 0.097
20 0.1355
50 0.3295
100 0.614