Friedman 검정 (반복측정)¶

import numpy as np
from scipy import stats
d = np.array([[30, 45, 60, 50], [25, 40, 55, 35], [35, 50, 65, 45],
              [20, 30, 50, 40], [40, 55, 70, 60]])
print(stats.friedmanchisquare(*d.T))
# statistic=13.56, pvalue=0.0035696

FriedmanchisquareResult(statistic=13.560000000000002, pvalue=0.0035695719978437604)

import itertools
from collections import Counter
from fractions import Fraction

b, k = 5, 4
perms = list(itertools.permutations(range(1, k + 1)))

# 블록을 하나씩 더해 가며 순위합 벡터의 분포를 쌓는다.
# 7962624 가지를 하나씩 돌지 않아도 같은 분포를 얻는다.
dist = Counter({(0,) * k: 1})
for _ in range(b):
    nxt = Counter()
    for state, c in dist.items():
        for perm in perms:
            nxt[tuple(s + q for s, q in zip(state, perm))] += c
    dist = nxt
total = sum(dist.values())
print(f"총 배정 수 {total}  (= 24^5 = {24 ** 5})")
print(f"R = (5, 10, 15, 20) 인 배정 {dist[(5, 10, 15, 20)]} 가지")

chis = Counter()
for R, c in dist.items():
    chis[round(12 / (b * k * (k + 1)) * sum(x * x for x in R) - 3 * b * (k + 1), 6)] += c
top = sorted(chis)[::-1]
print(f"\nchi2_F 가 가질 수 있는 값 {len(chis)} 개,  최댓값 {top[0]}")
print(" chi2_F   배정 수     누적      정확 p")
for v in top[:5]:
    cum = sum(c for w, c in chis.items() if w >= v - 1e-9)
    print(f"{v:>7.2f} {chis[v]:>9} {cum:>9}  {cum / total:.7f}")

obs = 13.56
cum = sum(c for w, c in chis.items() if w >= obs - 1e-9)
print(f"\n정확 p = {cum}/{total} = {Fraction(cum, total)} = {cum / total:.7f}")
print(f"chi2_3 p = {stats.chi2.sf(obs, 3):.7f}"
      f"  ->  {stats.chi2.sf(obs, 3) / (cum / total):.1f} 배")
print(f"도달 가능한 가장 작은 p = {chis[top[0]]}/{total}"
      f" = {Fraction(chis[top[0]], total)} = {chis[top[0]] / total:.2e}")
print(f"열거한 chi2_F 의 평균 {sum(v * c for v, c in chis.items()) / total}"
      f"  (chi2_3 의 평균 {k - 1})")

총 배정 수 7962624  (= 24^5 = 7962624)
R = (5, 10, 15, 20) 인 배정 1 가지

chi2_F 가 가질 수 있는 값 45 개,  최댓값 15.0
 chi2_F   배정 수     누적      정확 p
  15.00        24        24  0.0000030
  14.04       360       384  0.0000482
  13.56       720      1104  0.0001386
  13.08       600      1704  0.0002140
  12.84       960      2664  0.0003346

정확 p = 1104/7962624 = 23/165888 = 0.0001386
chi2_3 p = 0.0035696  ->  25.7 배
도달 가능한 가장 작은 p = 24/7962624 = 1/331776 = 3.01e-06
열거한 chi2_F 의 평균 3.0  (chi2_3 의 평균 3)

import numpy as np

rng = np.random.default_rng(0)
print(" b   m    R1   R2   chi2_F    (b-2m)^2/b       Z        Z^2")
for _ in range(4):
    b2 = 9
    x = rng.normal(0, 1, b2)
    y = rng.normal(0.3, 1, b2)
    m = int((x > y).sum())
    R1, R2 = b2 + m, 2 * b2 - m
    chi_F = 12 / (b2 * 2 * 3) * (R1 ** 2 + R2 ** 2) - 3 * b2 * 3
    z = (2 * m - b2) / np.sqrt(b2)
    print(f"{b2:>2} {m:>3} {R1:>5} {R2:>4} {chi_F:>9.6f} {(b2 - 2 * m) ** 2 / b2:>13.6f}"
          f" {z:>8.4f} {z ** 2:>10.6f}")

 b   m    R1   R2   chi2_F    (b-2m)^2/b       Z        Z^2
 9   7    16   11  2.777778      2.777778   1.6667   2.777778
 9   5    14   13  0.111111      0.111111   0.3333   0.111111
 9   3    12   15  1.000000      1.000000  -1.0000   1.000000
 9   3    12   15  1.000000      1.000000  -1.0000   1.000000

import numpy as np, itertools
b, k = 4, 3
perms = list(itertools.permutations(range(1, k + 1)))
cnt = tot = 0
for combo in itertools.product(perms, repeat=b):
    R = np.sum(np.array(combo), axis=0)
    chi = 12 / (b * k * (k + 1)) * np.sum(R ** 2) - 3 * b * (k + 1)
    tot += 1
    cnt += chi >= 2.0 - 1e-9
print(cnt / tot, tot)      # 0.4306  1296

0.4305555555555556 1296

import numpy as np
from scipy import stats
d = np.array([[30, 45, 60, 50], [25, 40, 55, 35], [35, 50, 65, 45],
              [20, 30, 50, 40], [40, 55, 70, 60]])
print(stats.friedmanchisquare(*d.T))      # 13.560, p=0.003570
print(stats.kruskal(*d.T))                # 12.237, p=0.006615
print(np.corrcoef(d.T)[0])                # 블록 간 상관

FriedmanchisquareResult(statistic=13.560000000000002, pvalue=0.0035695719978437604)
KruskalResult(statistic=12.236750189825356, pvalue=0.006614630461599845)
[1.         0.98639392 1.         0.82199494]

import numpy as np
from scipy import stats

# 자료 1: 모든 블록에서 효과가 일정하다
d1 = np.array([[10, 11, 12], [20, 21, 22], [30, 31, 32],
               [40, 41, 42], [50, 51, 52], [60, 61, 62]])
# 자료 2: 한 블록의 효과만 압도적으로 크다
d2 = np.array([[10, 11, 12], [20, 21, 22], [30, 31, 32],
               [40, 41, 42], [50, 51, 52], [60, 160, 260]])

for name, d in (("d1", d1), ("d2", d2)):
    print(name, stats.friedmanchisquare(*d.T))

d1 FriedmanchisquareResult(statistic=12.0, pvalue=0.002478752176666357)
d2 FriedmanchisquareResult(statistic=12.0, pvalue=0.002478752176666357)

import numpy as np
from scipy import stats
# 블록마다 동점이 하나씩 있는 자료
d = np.array([[5, 5, 8], [6, 6, 9], [4, 4, 7], [3, 3, 6], [7, 7, 10]])
print(stats.friedmanchisquare(*d.T))

FriedmanchisquareResult(statistic=10.0, pvalue=0.006737946999085468)