Kendall의 타우¶

import numpy as np
from scipy import stats

# x 는 오름차순인데 y 는 올랐다 내려간다. 일치쌍과 불일치쌍이 섞인 자료다.
x = np.array([1, 2, 3, 4, 5])
y = np.array([3, 5, 4, 2, 1])

# 타우는 (일치쌍 - 불일치쌍) / 전체 쌍이다. "무작위로 두 점을 골랐을 때
# 같은 방향으로 움직일 확률에서 반대 방향일 확률을 뺀 값"으로 읽을 수 있다.
tau, p_value = stats.kendalltau(x, y)
print(f"Kendall tau-b = {tau:.4f}, p-value = {p_value:.4f}")

# 스피어만과 견준다. 둘 다 순위만 쓰지만 눈금이 달라 타우 쪽이 늘 더 작다.
# 표본이 작거나 이상치가 있을 때는 타우가 더 안정적이다.
r_s, p_s = stats.spearmanr(x, y)
print(f"Spearman r_s  = {r_s:.4f}, p-value = {p_s:.4f}")

# (1) 간편 공식으로 r_s 를 직접
d = stats.rankdata(x) - stats.rankdata(y)
n = len(x)
print(f"\nd = {d.astype(int).tolist()}   sum d^2 = {int((d**2).sum())}"
      f"   r_s = {1 - 6 * (d**2).sum() / (n * (n**2 - 1)):.4f}")
print(f"3*tau - 2*r_s = {3 * tau - 2 * r_s:.4f}")

# (2) n! 가지 순열을 전부 돌려 |3tau - 2r_s| 의 최댓값을 본다
from itertools import permutations

for m in (5, 6, 7):
    xs = np.arange(1, m + 1)
    worst = max(abs(3 * stats.kendalltau(xs, q).statistic
                    - 2 * stats.spearmanr(xs, q).statistic)
                for q in permutations(xs))
    print(f"n = {m}:  max |3tau - 2r_s| = {worst:.6f}")

# (3) 정확 p-값. 120 가지를 모두 센다.
perms = list(permutations(range(1, 6)))
taus = np.array([stats.kendalltau(x, q).statistic for q in perms])
rss = np.array([stats.spearmanr(x, q).statistic for q in perms])
n_tau = int((np.abs(taus) >= abs(tau) - 1e-9).sum())
n_rs = int((np.abs(rss) >= abs(r_s) - 1e-9).sum())
print(f"\n|tau| >= 0.6 인 배열 {n_tau}/120   정확 p = {n_tau / 120:.6f}"
      f"   scipy = {p_value:.6f}")
print(f"|r_s| >= 0.7 인 배열 {n_rs}/120   정확 p = {n_rs / 120:.6f}"
      f"   scipy = {p_s:.6f}")
print(f"가장 작은 양측 p-값 2/120 = {2 / 120:.6f}")

Kendall tau-b = -0.6000, p-value = 0.2333
Spearman r_s  = -0.7000, p-value = 0.1881

d = [-2, -3, -1, 2, 4]   sum d^2 = 34   r_s = -0.7000
3*tau - 2*r_s = -0.4000
n = 5:  max |3tau - 2r_s| = 1.000000
n = 6:  max |3tau - 2r_s| = 1.000000
n = 7:  max |3tau - 2r_s| = 1.000000

|tau| >= 0.6 인 배열 28/120   정확 p = 0.233333   scipy = 0.233333
|r_s| >= 0.7 인 배열 28/120   정확 p = 0.233333   scipy = 0.188120
가장 작은 양측 p-값 2/120 = 0.016667

import warnings
warnings.filterwarnings("ignore")

import numpy as np

rng = np.random.default_rng(21001)
for rho in [0.0, 0.5, 0.8]:
    z = rng.standard_normal((400, 2))
    x = z[:, 0]
    y = rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]
    conc = disc = 0
    for i in range(400):
        for j in range(i + 1, 400):
            s = np.sign(x[i] - x[j]) * np.sign(y[i] - y[j])
            if s > 0:
                conc += 1
            elif s < 0:
                disc += 1
    tot = conc + disc
    print(f"  ρ={rho:.1f}: P(일치)={conc / tot:.4f}, "
          f"P(불일치)={disc / tot:.4f}, 차이={(conc - disc) / tot:.4f}, "
          f"τ 이론={2 / np.pi * np.arcsin(rho):.4f}")

  ρ=0.0: P(일치)=0.4642, P(불일치)=0.5358, 차이=-0.0716, τ 이론=0.0000
  ρ=0.5: P(일치)=0.6773, P(불일치)=0.3227, 차이=0.3546, τ 이론=0.3333
  ρ=0.8: P(일치)=0.8013, P(불일치)=0.1987, 차이=0.6025, τ 이론=0.5903

"τ = 0.4 입니다"
  → "두 사람을 무작위로 골랐을 때, 한쪽이 X 가 크면 Y 도 클 확률이
     70% 입니다. 반대인 경우가 30% 이고요."

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

x = np.array([1, 2, 2, 3, 4, 4, 4, 5, 6, 7], float)
y = np.array([2, 1, 3, 3, 5, 4, 6, 5, 8, 7], float)
print(f"x = {x.astype(int).tolist()}")
print(f"y = {y.astype(int).tolist()}")
for v in ["b", "c"]:
    r = stats.kendalltau(x, y, variant=v)
    print(f"  tau-{v}: {r.statistic:+.6f}  (p = {r.pvalue:.4f})")

n = len(x)
C = D = Tx = Ty = 0
for i in range(n):
    for j in range(i + 1, n):
        dx, dy = np.sign(x[i] - x[j]), np.sign(y[i] - y[j])
        if dx == 0 and dy == 0:
            pass                      # 양쪽 동점은 어디에도 세지 않는다
        elif dx == 0:
            Tx += 1
        elif dy == 0:
            Ty += 1
        elif dx * dy > 0:
            C += 1
        else:
            D += 1
print(f"\n수동 계산:  일치 C={C}, 불일치 D={D}, x 동점={Tx}, y 동점={Ty}")
print(f"  tau-a = (C-D)/[n(n-1)/2] = {(C - D) / (n * (n - 1) / 2):+.6f}")
print(f"  tau-b = (C-D)/√[(C+D+Tx)(C+D+Ty)] = "
      f"{(C - D) / np.sqrt((C + D + Tx) * (C + D + Ty)):+.6f}")

x = [1, 2, 2, 3, 4, 4, 4, 5, 6, 7]
y = [2, 1, 3, 3, 5, 4, 6, 5, 8, 7]
  tau-b: +0.785937  (p = 0.0024)
  tau-c: +0.770000  (p = 0.0024)

수동 계산:  일치 C=36, 불일치 D=3, x 동점=4, y 동점=2
  tau-a = (C-D)/[n(n-1)/2] = +0.733333
  tau-b = (C-D)/√[(C+D+Tx)(C+D+Ty)] = +0.785937

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(21004)
print("arctanh(τ̂) 의 실제 SD 와 피엘러 공식")
print(f"{'ρ':>6s} {'n':>5s} {'실제 SD':>9s} {'√(0.437/(n-4))':>15s}")
for rho in [0.0, 0.5, 0.8]:
    for n in [20, 50]:
        ts = []
        for _ in range(4_000):
            z = rng.standard_normal((n, 2))
            x = z[:, 0]
            y = rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]
            ts.append(np.arctanh(stats.kendalltau(x, y).statistic))
        print(f"{rho:6.2f} {n:5d} {np.std(ts, ddof=1):9.4f} "
              f"{np.sqrt(0.437 / (n - 4)):15.4f}")

arctanh(τ̂) 의 실제 SD 와 피엘러 공식
     ρ     n     실제 SD  √(0.437/(n-4))
  0.00    20    0.1666          0.1653
  0.00    50    0.0994          0.0975
  0.50    20    0.1633          0.1653
  0.50    50    0.0961          0.0975
  0.80    20    0.1659          0.1653
  0.80    50    0.0962          0.0975

def tau_ci_jack(x, y, a=0.05):
    """잭나이프 표준오차를 쓴 켄들 구간."""
    n = len(x)
    t = stats.kendalltau(x, y).statistic
    ts = np.empty(n)
    for i in range(n):
        m = np.ones(n, bool)
        m[i] = False
        ts[i] = stats.kendalltau(x[m], y[m]).statistic
    se = np.sqrt((n - 1) / n * ((ts - ts.mean())**2).sum())
    q = stats.norm.ppf(1 - a / 2)
    return t - q * se, t + q * se

rng = np.random.default_rng(21004)
B = 2_000
print("피복확률 (명목 95%)")
print(f"{'ρ':>6s} {'n':>5s} {'잭나이프':>9s} {'부트스트랩':>11s}")
for rho in [0.0, 0.5, 0.8]:
    for n in [20, 50]:
        kt_th = 2 / np.pi * np.arcsin(rho)
        a = b = 0
        for _ in range(B):
            z = rng.standard_normal((n, 2))
            x = z[:, 0]
            y = rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]
            lo, hi = tau_ci_jack(x, y)
            a += lo <= kt_th <= hi
            bs = [stats.kendalltau(x[i], y[i]).statistic
                  for i in (rng.integers(0, n, n) for _ in range(200))]
            lo, hi = np.quantile(bs, [0.025, 0.975])
            b += lo <= kt_th <= hi
        print(f"{rho:6.2f} {n:5d} {a / B:9.4f} {b / B:11.4f}")

피복확률 (명목 95%)
     ρ     n      잭나이프       부트스트랩
  0.00    20    0.9380      0.9450
  0.00    50    0.9505      0.9405
  0.50    20    0.9350      0.9510
  0.50    50    0.9575      0.9540
  0.80    20    0.9445      0.9595
  0.80    50    0.9490      0.9535

import numpy as np
from scipy import stats

rng = np.random.default_rng(20005)
n = 20
z = rng.standard_normal((n, 2))
x, y = z[:, 0], z[:, 1]
print(f"원 자료 (n={n}): 피어슨 {np.corrcoef(x, y)[0, 1]:+.4f}, "
      f"스피어만 {stats.spearmanr(x, y).statistic:+.4f}, "
      f"켄들 {stats.kendalltau(x, y).statistic:+.4f}")

print(f"\n한 점 (a, a) 를 추가하면")
print(f"{'a':>6s} {'피어슨':>9s} {'스피어만':>10s} {'켄들':>9s}")
for a in [3, 10, 100]:
    X, Y = np.append(x, a), np.append(y, a)
    print(f"{a:6d} {np.corrcoef(X, Y)[0, 1]:+9.4f} "
          f"{stats.spearmanr(X, Y).statistic:+10.4f} "
          f"{stats.kendalltau(X, Y).statistic:+9.4f}")

print(f"\n단조 비선형: y = exp(3x)")
xs = np.linspace(-1, 1, 30)
ys = np.exp(3 * xs)
print(f"  피어슨 {np.corrcoef(xs, ys)[0, 1]:+.4f}, "
      f"스피어만 {stats.spearmanr(xs, ys).statistic:+.4f}, "
      f"켄들 {stats.kendalltau(xs, ys).statistic:+.4f}")

print(f"\n단조 비선형: y = x^9")
ys = xs**9
print(f"  피어슨 {np.corrcoef(xs, ys)[0, 1]:+.4f}, "
      f"스피어만 {stats.spearmanr(xs, ys).statistic:+.4f}, "
      f"켄들 {stats.kendalltau(xs, ys).statistic:+.4f}")

print(f"\n비단조: y = x²")
ys = xs**2
print(f"  피어슨 {np.corrcoef(xs, ys)[0, 1]:+.4f}, "
      f"스피어만 {stats.spearmanr(xs, ys).statistic:+.4f}, "
      f"켄들 {stats.kendalltau(xs, ys).statistic:+.4f}")

원 자료 (n=20): 피어슨 +0.3119, 스피어만 +0.3188, 켄들 +0.2000

한 점 (a, a) 를 추가하면
     a       피어슨       스피어만        켄들
     3   +0.5816    +0.4117   +0.2762
    10   +0.9049    +0.4117   +0.2762
   100   +0.9988    +0.4117   +0.2762

단조 비선형: y = exp(3x)
  피어슨 +0.8131, 스피어만 +1.0000, 켄들 +1.0000

단조 비선형: y = x^9
  피어슨 +0.6944, 스피어만 +1.0000, 켄들 +1.0000

비단조: y = x²
  피어슨 -0.0000, 스피어만 -0.0127, 켄들 -0.0070

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(21006)
print("이변량 정규에서 두 측도의 비")
print(f"{'ρ':>6s} {'스피어만':>9s} {'켄들':>8s} {'r_s/τ':>8s} {'3/2':>6s}")
for rho in [0.2, 0.4, 0.6, 0.8, 0.95]:
    z = rng.standard_normal((5_000, 2))
    x = z[:, 0]
    y = rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]
    rs = stats.spearmanr(x, y).statistic
    tk = stats.kendalltau(x, y).statistic
    print(f"{rho:6.2f} {rs:9.4f} {tk:8.4f} {rs / tk:8.4f} {1.5:6.2f}")

print("\n검정력 비교 (ρ=0.3, 정규)")
B = 6_000
print(f"{'n':>5s} {'스피어만':>9s} {'켄들':>8s}")
for n in [20, 50, 100]:
    b = c = 0
    for _ in range(B):
        z = rng.standard_normal((n, 2))
        x = z[:, 0]
        y = 0.3 * z[:, 0] + np.sqrt(0.91) * z[:, 1]
        b += stats.spearmanr(x, y).pvalue < 0.05
        c += stats.kendalltau(x, y).pvalue < 0.05
    print(f"{n:5d} {b / B:9.4f} {c / B:8.4f}")

이변량 정규에서 두 측도의 비
     ρ      스피어만       켄들    r_s/τ    3/2
  0.20    0.1854   0.1242   1.4929   1.50
  0.40    0.3884   0.2640   1.4710   1.50
  0.60    0.5659   0.3978   1.4225   1.50
  0.80    0.7922   0.5963   1.3286   1.50
  0.95    0.9451   0.7970   1.1858   1.50

검정력 비교 (ρ=0.3, 정규)
    n      스피어만       켄들
   20    0.2262   0.2147
   50    0.5280   0.5235
  100    0.8310   0.8272

stats.kendalltau(x, y)                  → τ_b (기본값), p
stats.kendalltau(x, y, variant="c")     → τ_c
stats.kendalltau(x, y, method="exact")  → 정확 p (동점이 없을 때)

주의: 동점이 있으면 method="exact" 를 쓸 수 없다

교육 수준(6단계)과 소득 구간(8단계)의 연관

  켄들 τ_b = 0.34,  95% CI [0.27, 0.41],  n = 420,  p < 0.001
  → 두 사람을 무작위로 뽑으면 67% 에서 순서가 일치한다
  (대응하는 피어슨 척도로는 약 0.51)

  두 변수 모두 서열 척도이고 동점이 많아 τ_b 를 사용했다.