Spearman 순위상관 (재론)¶

import numpy as np
from scipy import stats
X = np.array([2, 5, 3, 8, 10, 8, 10, 6])
Y = np.array([7, 8, 6, 9, 10, 8, 6, 9])
print(stats.rankdata(X))   # [1.  3.  2.  5.5 7.5 5.5 7.5 4. ]
print(stats.rankdata(Y))   # [3.  4.5 1.5 6.5 8.  4.5 1.5 6.5]
print(stats.spearmanr(X, Y))
# SignificanceResult(statistic=0.37424, pvalue=0.36106)

[1.  3.  2.  5.5 7.5 5.5 7.5 4. ]
[3.  4.5 1.5 6.5 8.  4.5 1.5 6.5]
SignificanceResult(statistic=0.37424017169809104, pvalue=0.36106491276336194)

R, S = stats.rankdata(X), stats.rankdata(Y)
n = len(X)
SS_R = ((R - R.mean()) ** 2).sum()
SS_S = ((S - S.mean()) ** 2).sum()
SS_RS = ((R - R.mean()) * (S - S.mean())).sum()
d2 = ((R - S) ** 2).sum()
print(f"R 의 평균 {R.mean()},  S 의 평균 {S.mean()}  ((n+1)/2 = {(n + 1) / 2})")
print(f"SS_R = {SS_R},  SS_S = {SS_S},  SS_RS = {SS_RS},  sum d^2 = {d2}")
print(f"동점이 없을 때의 제곱합 n(n^2-1)/12 = {n * (n * n - 1) / 12}")

rs = SS_RS / np.sqrt(SS_R * SS_S)
print(f"\n완전한 공식   r_s = {SS_RS}/sqrt({SS_R}*{SS_S}) = {rs:.6f}")
print(f"항등식으로도  r_s = ({SS_R}+{SS_S}-{d2})/(2*{np.sqrt(SS_R * SS_S):.6f})"
      f" = {(SS_R + SS_S - d2) / (2 * np.sqrt(SS_R * SS_S)):.6f}")
print(f"축약 공식     r_s = 1 - 6*{d2}/({n}*{n * n - 1}) = "
      f"{1 - 6 * d2 / (n * (n * n - 1)):.6f}   <- 틀린 값")

t = rs * np.sqrt((n - 2) / (1 - rs ** 2))
print(f"\nt = {t:.6f},  양측 p = {2 * stats.t.sf(abs(t), n - 2):.8f}")
print(f"scipy 의 p-값                  = {stats.spearmanr(X, Y).pvalue:.8f}")

R 의 평균 4.5,  S 의 평균 4.5  ((n+1)/2 = 4.5)
SS_R = 41.0,  SS_S = 40.5,  SS_RS = 15.25,  sum d^2 = 51.0
동점이 없을 때의 제곱합 n(n^2-1)/12 = 42.0

완전한 공식   r_s = 15.25/sqrt(41.0*40.5) = 0.374240
항등식으로도  r_s = (41.0+40.5-51.0)/(2*40.749233) = 0.374240
축약 공식     r_s = 1 - 6*51.0/(8*63) = 0.392857   <- 틀린 값

t = 0.988532,  양측 p = 0.36106491
scipy 의 p-값                  = 0.36106491

import itertools
from math import factorial

print("완전 단조인 자료에 spearmanr 을 넣으면")
cases = [("n=8, y = x^3+1", np.arange(8), np.arange(8) ** 3 + 1),
         ("n=8, y = -2x   ", np.arange(8), -2.0 * np.arange(8)),
         ("n=5, y = e^x   ", np.arange(5), np.exp(np.arange(5)))]
for label, x, y in cases:
    r = stats.spearmanr(x, y)
    print(f"  {label}: statistic = {r.statistic!r},  pvalue = {r.pvalue!r}")

print("\n전수 열거로 구한 올바른 값")
for m in (5, 8):
    base = np.arange(1, m + 1)
    rs_all = np.array([np.corrcoef(base, p)[0, 1]
                       for p in itertools.permutations(base)])
    hit = int((np.abs(rs_all) >= 1 - 1e-12).sum())
    print(f"  n = {m}: 순열 {len(rs_all)} 개 (= {m}!),"
          f"  |r_s| = 1 인 것 {hit} 개  ->  p = 2/{factorial(m)}"
          f" = {hit / len(rs_all):.8f}")

완전 단조인 자료에 spearmanr 을 넣으면
  n=8, y = x^3+1: statistic = 1.0,  pvalue = 0.0
  n=8, y = -2x   : statistic = -1.0,  pvalue = 0.0
  n=5, y = e^x   : statistic = 0.9999999999999999,  pvalue = 1.4042654220543672e-24

전수 열거로 구한 올바른 값
  n = 5: 순열 120 개 (= 5!),  |r_s| = 1 인 것 2 개  ->  p = 2/120 = 0.01666667
  n = 8: 순열 40320 개 (= 8!),  |r_s| = 1 인 것 2 개  ->  p = 2/40320 = 0.00004960

import numpy as np
from scipy import stats
math = np.arange(1, 9); eng = np.array([3, 1, 2, 5, 4, 8, 6, 7])
print(stats.spearmanr(math, eng))
# SignificanceResult(statistic=0.83333, pvalue=0.010176)

SignificanceResult(statistic=0.8333333333333335, pvalue=0.01017554012345675)

import numpy as np, itertools
from scipy import stats

n = 8
R = np.arange(1, n + 1)
denom = n * (n * n - 1) / 6
rs_all = np.array([1 - ((R - np.array(p)) ** 2).sum() / denom
                   for p in itertools.permutations(R)])
print(len(rs_all), rs_all.mean().round(6), rs_all.std().round(4))

for obs in (0.8333, 0.6190, 0.4762, 0.3810):
    exact = 2 * (rs_all >= obs - 1e-9).mean()
    t = obs * np.sqrt((n - 2) / (1 - obs ** 2))
    approx = 2 * stats.t.sf(t, n - 2)
    print(round(obs, 4), round(exact, 5), round(approx, 5))

40320 0.0 0.378
0.8333 0.01538 0.01018
0.619 0.11498 0.10177
0.4762 0.21617 0.23292
0.381 0.32684 0.35175

import numpy as np
from scipy import stats
rng = np.random.default_rng(0)
n = 30
x = rng.normal(0, 1, n)
y = x + rng.normal(0, 0.5, n)      # 강한 양의 상관

print("원자료:", stats.pearsonr(x, y).statistic.round(3),
      stats.spearmanr(x, y).statistic.round(3))

# (A) 두 변수 모두 극단이지만 추세를 따르는 점
xa, ya = np.append(x, 10), np.append(y, 10)
print("추세 따름:", stats.pearsonr(xa, ya).statistic.round(3),
      stats.spearmanr(xa, ya).statistic.round(3))

# (B) 추세를 거스르는 지렛대점
xb, yb = np.append(x, 10), np.append(y, -10)
print("추세 거스름:", stats.pearsonr(xb, yb).statistic.round(3),
      stats.spearmanr(xb, yb).statistic.round(3))

원자료: 0.829 0.84
추세 따름: 0.972 0.855
추세 거스름: -0.701 0.668

import numpy as np
from scipy import stats
x = np.arange(1, 21)
y = np.exp(x / 2)          # 완전한 단조증가
print(stats.pearsonr(x, y).statistic.round(4))    # 0.6991
print(stats.spearmanr(x, y).statistic.round(4))   # 1.0000

0.6991
1.0

x = np.linspace(-3, 3, 41)
y = x ** 2                 # 완벽한 결정론적 관계
print(stats.pearsonr(x, y).statistic.round(4))    # -0.0000
print(stats.spearmanr(x, y).statistic.round(4))   # -0.0146

-0.0
-0.0146