Pearson 상관계수¶

import numpy as np
from scipy import stats

# 예시 자료
x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
y = np.array([2.1, 3.9, 6.2, 7.8, 10.1, 12.3, 13.8, 16.1, 18.0, 20.2])

# 방법 1: NumPy 의 상관행렬
corr_matrix = np.corrcoef(x, y)
r_numpy = corr_matrix[0, 1]
print(f"NumPy r = {r_numpy:.4f}")

# 방법 2: SciPy (p-값도 함께 돌려준다)
r_scipy, p_value = stats.pearsonr(x, y)
print(f"SciPy r = {r_scipy:.4f}, p-value = {p_value:.6f}")

# 이상점 하나를 더하면 어떻게 되는지 본다.
x2 = np.append(x, 10)
y2 = np.append(y, -5)
print(f"이상점 추가 후 r = {np.corrcoef(x2, y2)[0, 1]:.4f}")

# (1) 세 제곱합과 닫힌 꼴
n = len(x)
Sxx = ((x - x.mean())**2).sum()
Syy = ((y - y.mean())**2).sum()
Sxy = ((x - x.mean()) * (y - y.mean())).sum()
print(f"\nSxx = {Sxx}   n(n^2-1)/12 = {n * (n**2 - 1) / 12}")
print(f"Syy = {Syy}   Sxy = {Sxy}")
print(f"r = {Sxy / np.sqrt(Sxx * Syy):.7f}   기울기 = {Sxy / Sxx:.6f}"
      f"   r^2 = {(Sxy / np.sqrt(Sxx * Syy))**2:.7f}")

# (2) 척도불변. 공분산은 그렇지 않다.
xs, ys = x + 100, 2.20462 * y + 37
print(f"\n변환 전 r {np.corrcoef(x, y)[0, 1]:.16f}  공분산 {np.cov(x, y)[0, 1]:.6f}")
print(f"변환 후 r {np.corrcoef(xs, ys)[0, 1]:.16f}  공분산 {np.cov(xs, ys)[0, 1]:.6f}")

# (3) 이상점 한 점의 기여
bx, by = x2.mean(), y2.mean()
Sxy2 = ((x2 - bx) * (y2 - by)).sum()
Syy2 = ((y2 - by)**2).sum()
Sxx2 = ((x2 - bx)**2).sum()
cxy = (x2[-1] - bx) * (y2[-1] - by)
print(f"\n새 평균 ({bx:.6f}, {by:.6f})")
print(f"그 점의 Sxy 기여 {cxy:.6f}  /  Sxy = {Sxy2:.6f}   ({cxy / Sxy2:+.1%})")
print(f"그 점의 Syy 기여 {(y2[-1] - by)**2:.6f}  /  Syy = {Syy2:.6f}"
      f"   ({(y2[-1] - by)**2 / Syy2:.1%})")
print(f"그 점의 Sxx 기여 {(x2[-1] - bx)**2:.6f}  /  Sxx = {Sxx2:.6f}"
      f"   ({(x2[-1] - bx)**2 / Sxx2:.1%})")
print(f"r' = {Sxy2 / np.sqrt(Sxx2 * Syy2):.6f}")

NumPy r = 0.9996
SciPy r = 0.9996, p-value = 0.000000
이상점 추가 후 r = 0.4180

Sxx = 82.5   n(n^2-1)/12 = 82.5
Syy = 332.865   Sxy = 165.64999999999998
r = 0.9996096   기울기 = 2.007879   r^2 = 0.9992193

변환 전 r 0.9996095571302254  공분산 18.405556
변환 후 r 0.9996095571302255  공분산 40.577256

새 평균 (5.909091, 9.590909)
그 점의 Sxy 기여 -59.690083  /  Sxy = 99.990909   (-59.7%)
그 점의 Syy 기여 212.894628  /  Syy = 567.049091   (37.5%)
그 점의 Sxx 기여 16.735537  /  Sxx = 100.909091   (16.6%)
r' = 0.418008

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(20001)
print(f"{'ρ':>6s} {'n':>5s} {'E[r]':>9s} {'편향':>9s} {'SD':>8s} "
      f"{'치우침':>8s} {'이론 SD':>8s}")
for rho in [0.0, 0.5, 0.8, 0.95]:
    for n in [10, 30]:
        rs = []
        for _ in range(10_000):
            z = rng.standard_normal((n, 2))
            x = z[:, 0]
            y = rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]
            rs.append(np.corrcoef(x, y)[0, 1])
        rs = np.array(rs)
        print(f"{rho:6.2f} {n:5d} {rs.mean():9.4f} {rs.mean() - rho:9.4f} "
              f"{rs.std(ddof=1):8.4f} {stats.skew(rs):8.4f} "
              f"{(1 - rho**2) / np.sqrt(n - 1):8.4f}")

     ρ     n      E[r]        편향       SD      치우침    이론 SD
  0.00    10    0.0047    0.0047   0.3365  -0.0398   0.3333
  0.00    30    0.0002    0.0002   0.1849  -0.0028   0.1857
  0.50    10    0.4752   -0.0248   0.2664  -0.8488   0.2500
  0.50    30    0.4954   -0.0046   0.1429  -0.5994   0.1393
  0.80    10    0.7835   -0.0165   0.1474  -1.9197   0.1200
  0.80    30    0.7951   -0.0049   0.0715  -0.9843   0.0669
  0.95    10    0.9436   -0.0064   0.0456  -3.7411   0.0325
  0.95    30    0.9482   -0.0018   0.0199  -1.3341   0.0181

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

def fisher_ci(r, n, a=0.05):
    z = np.arctanh(r)
    se = 1 / np.sqrt(n - 3)
    q = stats.norm.ppf(1 - a / 2)
    return np.tanh(z - q * se), np.tanh(z + q * se)

def naive_ci(r, n, a=0.05):
    se = (1 - r**2) / np.sqrt(n - 1)
    q = stats.norm.ppf(1 - a / 2)
    return r - q * se, r + q * se

rng = np.random.default_rng(20002)
B = 5_000
print("참 ρ 를 덮는 비율 (명목 95%)")
print(f"{'ρ':>6s} {'n':>5s} {'피셔 z':>9s} {'단순 근사':>10s}")
for rho in [0.0, 0.5, 0.8, 0.95]:
    for n in [10, 30, 100]:
        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]
            r = np.corrcoef(x, y)[0, 1]
            lo, hi = fisher_ci(r, n)
            a += lo <= rho <= hi
            lo, hi = naive_ci(r, n)
            b += lo <= rho <= hi
        print(f"{rho:6.2f} {n:5d} {a / B:9.4f} {b / B:10.4f}")

참 ρ 를 덮는 비율 (명목 95%)
     ρ     n      피셔 z      단순 근사
  0.00    10    0.9468     0.8530
  0.00    30    0.9486     0.9216
  0.00   100    0.9488     0.9388
  0.50    10    0.9550     0.8610
  0.50    30    0.9492     0.9188
  0.50   100    0.9464     0.9352
  0.80    10    0.9500     0.8518
  0.80    30    0.9458     0.9124
  0.80   100    0.9520     0.9416
  0.95    10    0.9588     0.8462
  0.95    30    0.9526     0.9114
  0.95   100    0.9510     0.9356

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]          # 참 상관 0
print(f"원 자료 (n={n}):  피어슨 r = {np.corrcoef(x, y)[0, 1]:+.4f},  "
      f"스피어만 = {stats.spearmanr(x, y).statistic:+.4f}")

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

print("\n반대 방향의 점 하나")
print(f"{'추가한 점 (a, -a)':>18s} {'피어슨 r':>10s} {'스피어만':>10s} {'켄들':>9s}")
for a in [10, 50, 100]:
    X, Y = np.append(x, a), np.append(y, -a)
    print(f"{f'({a}, {-a})':>18s} {np.corrcoef(X, Y)[0, 1]:+10.4f} "
          f"{stats.spearmanr(X, Y).statistic:+10.4f} "
          f"{stats.kendalltau(X, Y).statistic:+9.4f}")

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

      추가한 점 (a, a)      피어슨 r       스피어만        켄들
            (3, 3)    +0.5816    +0.4117   +0.2762
            (5, 5)    +0.7415    +0.4117   +0.2762
          (10, 10)    +0.9049    +0.4117   +0.2762
          (20, 20)    +0.9726    +0.4117   +0.2762
          (50, 50)    +0.9954    +0.4117   +0.2762
        (100, 100)    +0.9988    +0.4117   +0.2762

반대 방향의 점 하나
     추가한 점 (a, -a)      피어슨 r       스피어만        켄들
         (10, -10)    -0.8123    +0.1390   +0.0857
         (50, -50)    -0.9915    +0.1390   +0.0857
       (100, -100)    -0.9979    +0.1390   +0.0857

r 을 계산한다
각 i 에 대해 i 를 뺀 r_(i) 를 계산한다
max|r - r_(i)| 를 본다
  → 0.1 을 넘으면 그 점을 조사한다

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(20006)
B = 4_000
rho = 0.6
print("모집단 ρ=0.60. x 의 상위 p 분위만 뽑으면 r 이 얼마가 되는가")
print(f"{'선발 비율':>10s} {'평균 r':>9s} {'이론 예측':>10s}")
for frac in [1.0, 0.5, 0.3, 0.2, 0.1]:
    rs = []
    for _ in range(B):
        z = rng.standard_normal((2_000, 2))
        x = z[:, 0]
        y = rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]
        m = x >= np.quantile(x, 1 - frac)
        if m.sum() > 5:
            rs.append(np.corrcoef(x[m], y[m])[0, 1])
    if frac < 1.0:
        z0 = stats.norm.ppf(1 - frac)
        lam = stats.norm.pdf(z0) / frac
        k = np.sqrt(1 - lam * (lam - z0))       # 절단정규의 SD 비
        pred = rho * k / np.sqrt(1 - rho**2 + rho**2 * k**2)
    else:
        pred = rho
    print(f"{frac:10.0%} {np.mean(rs):9.4f} {pred:10.4f}")

모집단 ρ=0.60. x 의 상위 p 분위만 뽑으면 r 이 얼마가 되는가
     선발 비율      평균 r      이론 예측
      100%    0.5999     0.6000
       50%    0.4124     0.4120
       30%    0.3589     0.3599
       20%    0.3306     0.3309
       10%    0.2931     0.2947

관측: 합격생만 보면 수능 점수와 학점의 상관이 0.30
      → "수능은 대학 성적을 잘 예측하지 못한다"?

합격률이 20% 이고 모집단 SD 대비 합격생 SD 가 0.62 이면
  → 보정된 상관은 약 0.55
  → "중간 정도로 예측한다"

두 결론의 차이가 정책을 가른다.

import warnings
warnings.filterwarnings("ignore")

import numpy as np

rng = np.random.default_rng(20007)
rho = 0.6
print("참 상관 ρ=0.60, 신뢰도 r_xx, r_yy 일 때 관측 상관")
print(f"{'r_xx':>7s} {'r_yy':>7s} {'관측 r (이론)':>13s} {'모의':>8s}")
for rxx, ryy in [(1.0, 1.0), (0.9, 0.9), (0.8, 0.8), (0.7, 0.7), (0.6, 0.5)]:
    pred = rho * np.sqrt(rxx * ryy)
    rs = []
    for _ in range(2_000):
        z = rng.standard_normal((500, 2))
        tx = z[:, 0]
        ty = rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]
        ox = tx + rng.normal(0, np.sqrt(1 / rxx - 1), 500)
        oy = ty + rng.normal(0, np.sqrt(1 / ryy - 1), 500)
        rs.append(np.corrcoef(ox, oy)[0, 1])
    print(f"{rxx:7.2f} {ryy:7.2f} {pred:13.4f} {np.mean(rs):8.4f}")

참 상관 ρ=0.60, 신뢰도 r_xx, r_yy 일 때 관측 상관
   r_xx    r_yy     관측 r (이론)       모의
   1.00    1.00        0.6000   0.6000
   0.90    0.90        0.5400   0.5405
   0.80    0.80        0.4800   0.4782
   0.70    0.70        0.4200   0.4203
   0.60    0.50        0.3286   0.3281

참 ρ=0.60 을 α=0.05, 검정력 0.80 으로 잡으려면
  신뢰도 1.0:  n ≈ 19
  신뢰도 0.7:  관측 ρ=0.42 → n ≈ 42
  신뢰도 0.5:  관측 ρ=0.30 → n ≈ 85

관측 상관 r = 0.42
신뢰도 r_xx = 0.72 (크론바흐 α), r_yy = 0.68 (재검사)
감쇠 보정 상관 = 0.42 / √(0.72×0.68) = 0.60

주의: 보정값은 신뢰도 추정의 불확실성을 반영하지 않는다.

import numpy as np
from scipy import stats

def n_corr(rho, power=0.8, alpha=0.05):
    z = np.arctanh(rho)
    return int(np.ceil(((stats.norm.ppf(1 - alpha / 2)
                         + stats.norm.ppf(power)) / z)**2 + 3))

print("α=0.05 양측, 검정력 0.80")
print(f"{'ρ':>6s} {'필요한 n':>9s} {'r²':>7s} {'해석':>10s}")
for rho, lab in [(0.1, "매우 작음"), (0.2, "작음"), (0.3, "중간"),
                 (0.4, ""), (0.5, "큼"), (0.7, ""), (0.9, "")]:
    print(f"{rho:6.1f} {n_corr(rho):9,d} {rho**2:7.3f} {lab:>10s}")

print("\nr = 0.3 을 관측했을 때 95% 피셔 구간")
for n in [20, 50, 100, 300, 1000]:
    z = np.arctanh(0.3)
    se = 1 / np.sqrt(n - 3)
    q = stats.norm.ppf(0.975)
    lo, hi = np.tanh(z - q * se), np.tanh(z + q * se)
    print(f"  n={n:5d}:  [{lo:+.3f}, {hi:+.3f}]   폭 {hi - lo:.3f}")

def n_diff(r1, r2, power=0.8, alpha=0.05):
    d = abs(np.arctanh(r1) - np.arctanh(r2))
    return int(np.ceil(2 * ((stats.norm.ppf(1 - alpha / 2)
                             + stats.norm.ppf(power)) / d)**2 + 3))

print("\n두 상관의 차이를 탐지하는 데 필요한 n (독립 표본, 각 군)")
for r1, r2 in [(0.3, 0.5), (0.5, 0.7), (0.2, 0.4), (0.6, 0.7), (0.0, 0.3)]:
    print(f"  r1={r1}, r2={r2}:  군당 n = {n_diff(r1, r2):,d}")

α=0.05 양측, 검정력 0.80
     ρ     필요한 n      r²         해석
   0.1       783   0.010      매우 작음
   0.2       194   0.040         작음
   0.3        85   0.090         중간
   0.4        47   0.160           
   0.5        30   0.250          큼
   0.7        14   0.490           
   0.9         7   0.810           

r = 0.3 을 관측했을 때 95% 피셔 구간
  n=   20:  [-0.164, +0.655]   폭 0.820
  n=   50:  [+0.024, +0.534]   폭 0.510
  n=  100:  [+0.110, +0.469]   폭 0.359
  n=  300:  [+0.193, +0.400]   폭 0.206
  n= 1000:  [+0.243, +0.355]   폭 0.113

두 상관의 차이를 탐지하는 데 필요한 n (독립 표본, 각 군)
  r1=0.3, r2=0.5:  군당 n = 277
  r1=0.5, r2=0.7:  군당 n = 159
  r1=0.2, r2=0.4:  군당 n = 325
  r1=0.6, r2=0.7:  군당 n = 521
  r1=0.0, r2=0.3:  군당 n = 167

□ 산점도를 그렸는가        ← 가장 중요
□ 관계가 선형인가          ← 곡선이면 r 이 오도한다
□ 이상점이 있는가          ← 점 하나가 r 을 지배할 수 있다
□ 범위가 제한되었는가      ← 선발·절단이 있는가
□ 측정 신뢰도는 얼마인가   ← 감쇠

키와 몸무게의 상관: r = 0.68, 95% CI [0.58, 0.76], n = 120
  (피어슨, 피셔 z 구간)
산점도에서 선형 관계가 확인되었고 뚜렷한 이상점은 없었다.
스피어만 순위상관도 0.66 으로 유사하다.
결정계수 r² = 0.46 — 몸무게 변동의 46% 를 키가 설명한다.