상관에 대한 순열검정¶

import numpy as np
from scipy import stats

temp = np.array([14.5, 15.6, 16.7, 17.1, 17.5, 17.5, 18.7, 18.7, 20.0, 21.3,
                 24.8, 25.9, 26.2, 27.5, 30.1, 30.1, 30.3, 31.4, 32.2, 32.5])
sales = np.array([276, 233, 288, 309, 328, 246, 322, 310, 309, 323,
                  332, 344, 415, 390, 389, 319, 377, 401, 381, 382])

r = np.corrcoef(temp, sales)[0, 1]
print(f"r = {r:.4f}")                      # 0.8268

rng = np.random.default_rng(0)
B = 9999
perm = np.array([np.corrcoef(temp, rng.permutation(sales))[0, 1]
                 for _ in range(B)])
p = ((np.abs(perm) >= abs(r)).sum() + 1) / (B + 1)
print(f"permutation p = {p:.5f}")          # 0.0001
print(stats.pearsonr(temp, sales))         # p = 7.4e-06

r = 0.8268
permutation p = 0.00010
PearsonRResult(statistic=0.8267919627184452, pvalue=6.976166744153747e-06)

n = len(temp)
print(f"초과 횟수 c = {int((np.abs(perm) >= abs(r)).sum())},  "
      f"최대 |r*| = {np.abs(perm).max():.4f}")
print(f"SD(r*) 닫힌 꼴 = {1 / np.sqrt(n - 1):.7f},  모의 = {perm.std(ddof=1):.7f}"
      f"   (몬테카를로 오차 {1 / np.sqrt(n - 1) / np.sqrt(2 * B):.7f})")
print(f"평균 r*  = {perm.mean():+.7f}   (평균의 오차 {perm.std(ddof=1) / np.sqrt(B):.7f})")

print(f"\n점추정  (0+1)/(B+1)       = {1 / (B + 1):.3e}")
print(f"95% 신뢰상한 1-0.05^(1/B) = {1 - 0.05 ** (1 / B):.3e}")
print(f"t 검정의 p                = {stats.pearsonr(temp, sales).pvalue:.3e}")

tstar = perm * np.sqrt((n - 2) / (1 - perm ** 2))
print(f"\nt_obs = {r * np.sqrt((n - 2) / (1 - r ** 2)):.4f},  "
      f"순열 |t*| 의 최대 = {np.abs(tstar).max():.4f}")
for q in (0.95, 0.99, 0.999):
    print(f"  분위 {q}: 순열 |t*| = {np.quantile(np.abs(tstar), q):.4f},"
          f"   t_18 = {stats.t(n - 2).ppf((1 + q) / 2):.4f}")

초과 횟수 c = 0,  최대 |r*| = 0.8118
SD(r*) 닫힌 꼴 = 0.2294157,  모의 = 0.2286138   (몬테카를로 오차 0.0016223)
평균 r*  = +0.0012602   (평균의 오차 0.0022863)

점추정  (0+1)/(B+1)       = 1.000e-04
95% 신뢰상한 1-0.05^(1/B) = 2.996e-04
t 검정의 p                = 6.976e-06

t_obs = 6.2360,  순열 |t*| 의 최대 = 5.8990
  분위 0.95: 순열 |t*| = 2.0931,   t_18 = 2.1009
  분위 0.99: 순열 |t*| = 2.9373,   t_18 = 2.8784
  분위 0.999: 순열 |t*| = 4.0474,   t_18 = 3.9216

import numpy as np, itertools
from scipy import stats

x = np.array([1.2, 2.5, 3.1, 4.8, 5.3, 6.9])
y = np.array([2.1, 2.9, 4.4, 4.0, 6.8, 7.5])

r = stats.pearsonr(x, y)
print(round(r.statistic, 4), round(r.pvalue, 4))      # 0.9243  0.0084

rs = np.array([np.corrcoef(x, y[list(p)])[0, 1]
               for p in itertools.permutations(range(6))])
print(len(rs), (np.abs(rs) >= abs(r.statistic) - 1e-12).mean())   # 720  0.00694

sp = stats.spearmanr(x, y)
print(round(sp.statistic, 4), round(sp.pvalue, 4))    # 0.9429  0.0048
rss = np.array([stats.spearmanr(x, y[list(p)]).statistic
                for p in itertools.permutations(range(6))])
print((np.abs(rss) >= abs(sp.statistic) - 1e-12).mean())          # 0.01667

0.9243 0.0084
720 0.006944444444444444
0.9429 0.0048
0.016666666666666666

import numpy as np
from scipy import stats
rng = np.random.default_rng(5)

def rows_corr(A, B):
    A = A - A.mean(1, keepdims=True); B = B - B.mean(1, keepdims=True)
    return (A*B).sum(1) / np.sqrt((A**2).sum(1) * (B**2).sum(1))

def size(gen, n=15, M=4000, B=999):
    ct = cp = 0
    for _ in range(M):
        x = gen(n); y = gen(n)
        ct += stats.pearsonr(x, y).pvalue < 0.05
        Y = np.array([rng.permutation(y) for _ in range(B)])
        allr = rows_corr(np.vstack([x, np.tile(x, (B, 1))]),
                         np.vstack([y, Y]))
        r, rp = abs(allr[0]), np.abs(allr[1:])
        cp += ((rp >= r - 1e-9).sum() + 1) / (B + 1) < 0.05
    return round(ct/M, 4), round(cp/M, 4)

import numpy as np
from scipy import stats
rng = np.random.default_rng(3)

def power(n=30, M=1500, B=499):
    cp = cs = cq = 0
    for _ in range(M):
        x = rng.normal(0, 1, n)
        y = x**2 + rng.normal(0, 0.5, n)
        # (1) 피어슨 상관
        r = np.corrcoef(x, y)[0, 1]
        Y = np.array([rng.permutation(y) for _ in range(B)])
        cp += ((np.abs(rows_corr(np.tile(x, (B,1)), Y)) >= abs(r)).sum()+1)/(B+1) < 0.05
        # (2) 스피어만 상관
        rs = stats.spearmanr(x, y).statistic
        rr = np.array([stats.spearmanr(x, rng.permutation(y)).statistic
                       for _ in range(B)])
        cs += ((np.abs(rr) >= abs(rs)).sum() + 1)/(B+1) < 0.05
        # (3) |x - med(x)| 과 y 의 상관
        a = np.abs(x - np.median(x))
        q = np.corrcoef(a, y)[0, 1]
        cq += ((np.abs(rows_corr(np.tile(a, (B,1)), Y)) >= abs(q)).sum()+1)/(B+1) < 0.05
    return round(cp/M,3), round(cs/M,3), round(cq/M,3)