공분산 밑바닥부터 만들기¶

import numpy as np

def covariance_step_by_step(x, y):
    """표본공분산을 구하고, 중간 계산인 편차까지 함께 돌려준다.

    편차를 돌려주는 까닭은 뒤에서 편차곱을 막대로 그려 보이기 위함이다.
    """
    n = len(x)
    x_mean = x.mean()
    y_mean = y.mean()
    x_dev = x - x_mean
    y_dev = y - y_mean
    # n이 아니라 n-1로 나눈다. 베셀 보정이며, 분산에서와 같은 이유다.
    # 편차를 참 평균이 아니라 표본평균에서 쟀기 때문에 자유도 하나를 잃는다.
    cov = np.sum(x_dev * y_dev) / (n - 1)
    return cov, x_dev, y_dev

def pearson_r_step_by_step(x, y):
    """피어슨 상관계수를 정의대로 구한다.

    공분산을 두 표준편차의 곱으로 나눈다. 이 나눗셈이 단위를 없애므로
    r 은 -1 과 1 사이에 갇힌 값이 된다.
    """
    cov, _, _ = covariance_step_by_step(x, y)
    # ddof=1로 맞춰야 한다. 공분산이 n-1로 나눈 값이므로
    # 표준편차도 같은 규약을 써야 두 n-1이 약분되어 r이 척도와 무관해진다.
    sx = x.std(ddof=1)
    sy = y.std(ddof=1)
    return cov / (sx * sy)

# 한 직선 위에 완벽히 놓인 자료. r 는 정확히 1 이어야 한다.
a = np.arange(10.0)
b = 2 * a + 1
n = len(a)
cov, _, _ = covariance_step_by_step(a, b)

print(f"n = {n}")
print(f"규약을 맞춘 r      = {pearson_r_step_by_step(a, b):.6f}")
print(f"분모만 n 으로 나눈 r = {cov / (a.std(ddof=0) * b.std(ddof=0)):.6f}"
      f"   (= r * n/(n-1) = {n / (n - 1):.6f})")

# 둘을 일관되게 쓰면 규약이 약분된다.
cov_n = np.sum((a - a.mean()) * (b - b.mean())) / n
print(f"\n둘 다 n 으로   r = {cov_n / (a.std(ddof=0) * b.std(ddof=0)):.6f}")
print(f"둘 다 n-1 로  r = {cov / (a.std(ddof=1) * b.std(ddof=1)):.6f}")

n = 10
규약을 맞춘 r      = 1.000000
분모만 n 으로 나눈 r = 1.111111   (= r * n/(n-1) = 1.111111)

둘 다 n 으로   r = 1.000000
둘 다 n-1 로  r = 1.000000

np.random.seed(42)
WEEKS = 48

# 두 주의 가격에 공통으로 실릴 하락 추세. 공분산이 커지는 까닭이 바로 이것이다.
trend = np.linspace(0, -12, WEEKS)

CA = 248.0 + trend + np.random.normal(0, 0.5, WEEKS)
NY = 350.0 + trend * 0.8 + np.random.normal(0, 0.6, WEEKS)

s_T2 = trend.var(ddof=1)
print(f"s_T^2 = {s_T2:.6f}   (12/47)^2 * 196 = {(12 / 47) ** 2 * 196:.6f}")
print(f"E[Cov]    = 0.8 s_T^2        = {0.8 * s_T2:.4f}")
print(f"E[s_CA^2] = s_T^2 + 0.25     = {s_T2 + 0.25:.4f}")
print(f"E[s_NY^2] = 0.64 s_T^2 + 0.36 = {0.64 * s_T2 + 0.36:.4f}")
print(f"대입해 얻은 rho              = {0.8 * s_T2 / np.sqrt((s_T2 + 0.25) * (0.64 * s_T2 + 0.36)):.4f}")

# 씨앗을 바꿔 가며 20만 번 되풀이해 기댓값을 직접 잰다.
rng = np.random.default_rng(99)
A = 248.0 + trend + rng.normal(0, 0.5, (200_000, WEEKS))
B = 350.0 + 0.8 * trend + rng.normal(0, 0.6, (200_000, WEEKS))
Ad = A - A.mean(1, keepdims=True)
Bd = B - B.mean(1, keepdims=True)
cv = (Ad * Bd).sum(1) / (WEEKS - 1)
rr = cv / (Ad.std(1, ddof=1) * Bd.std(1, ddof=1))
print(f"\n20만 회 평균 Cov = {cv.mean():.4f} (표준편차 {cv.std():.4f})")
print(f"20만 회 평균 r   = {rr.mean():.5f} (표준편차 {rr.std():.5f})")

s_T^2 = 12.776822   (12/47)^2 * 196 = 12.776822
E[Cov]    = 0.8 s_T^2        = 10.2215
E[s_CA^2] = s_T^2 + 0.25     = 13.0268
E[s_NY^2] = 0.64 s_T^2 + 0.36 = 8.5372
대입해 얻은 rho              = 0.9692

20만 회 평균 Cov = 10.2217 (표준편차 0.3788)
20만 회 평균 r   = 0.96986 (표준편차 0.00631)

import pandas as pd

cov, x_dev, y_dev = covariance_step_by_step(CA, NY)
r = pearson_r_step_by_step(CA, NY)

print(f"CA mean     = {CA.mean():.4f}")
print(f"NY mean     = {NY.mean():.4f}")
print(f"Covariance  = {cov:.4f}")
print(f"Pearson r   = {r:.4f}")

# 직접 구한 값이 라이브러리와 맞는지 확인한다. pandas 의 cov 는 ddof=1 이므로
# 위 구현도 n-1 로 나눠야 값이 맞는다.
df = pd.DataFrame({"CA": CA, "NY": NY})
print(f"pandas cov  = {df['CA'].cov(df['NY']):.4f}")
print(f"pandas corr = {df['CA'].corr(df['NY']):.4f}")
print(f"numpy corr  = {np.corrcoef(CA, NY)[0, 1]:.4f}")

# (1) 나누는 수의 규약
print(f"\nnp.cov 기본(n-1)      = {np.cov(CA, NY)[0, 1]:.4f}")
print(f"np.cov bias=True(n)   = {np.cov(CA, NY, bias=True)[0, 1]:.4f}"
      f"   = 위 값 x 47/48 = {np.cov(CA, NY)[0, 1] * 47 / 48:.4f}")

# (2) 보기 2 의 이론값·표집 표준편차에 대어 보기
E_cov, sd_cov = 10.2215, 0.3788
E_r, sd_r = 0.96986, 0.00631
print(f"\nCov  이 표본 {cov:.4f}  이론 {E_cov:.4f}"
      f"  차 {cov - E_cov:+.4f}  SD {sd_cov}  -> {(cov - E_cov) / sd_cov:+.2f}")
print(f"r    이 표본 {r:.5f}  이론 {E_r:.5f}"
      f"  차 {r - E_r:+.5f}  SD {sd_r}  -> {(r - E_r) / sd_r:+.2f}")
print(f"상대적 흔들림  Cov {sd_cov / E_cov:.4f}   r {sd_r / E_r:.4f}")

CA mean     = 241.8974
NY mean     = 345.1893
Covariance  = 10.5691
Pearson r   = 0.9753
pandas cov  = 10.5691
pandas corr = 0.9753
numpy corr  = 0.9753

np.cov 기본(n-1)      = 10.5691
np.cov bias=True(n)   = 10.3489   = 위 값 x 47/48 = 10.3489

Cov  이 표본 10.5691  이론 10.2215  차 +0.3476  SD 0.3788  -> +0.92
r    이 표본 0.97525  이론 0.96986  차 +0.00539  SD 0.00631  -> +0.85
상대적 흔들림  Cov 0.0371   r 0.0065

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))

# 왼쪽: 산점도와 회귀직선
axes[0].scatter(CA, NY, alpha=0.6, edgecolors="grey")
z = np.polyfit(CA, NY, 1)
axes[0].plot(np.sort(CA), np.polyval(z, np.sort(CA)),
             color="red", linewidth=2)
axes[0].set_xlabel("CA Price (\\$)")
axes[0].set_ylabel("NY Price (\\$)")
axes[0].set_title(f"Scatter (r = {r:.3f})")

# 가운데: 주마다의 편차곱. 공분산은 이 막대들의 평균이다.
# 파란 막대(양)가 빨간 막대(음)를 압도하면 공분산이 양이 된다.
products = x_dev * y_dev
colours = ["steelblue" if p > 0 else "salmon" for p in products]
axes[1].bar(range(WEEKS), products, color=colours, edgecolor="white")
axes[1].axhline(0, color="black", linewidth=0.5)
axes[1].set_xlabel("Week")
axes[1].set_ylabel("$(x - \\bar{x})(y - \\bar{y})$")
axes[1].set_title("Deviation Products")

# 오른쪽: 두 시계열. 함께 내려가는 모습이 위 편차곱의 부호를 설명한다.
weeks = np.arange(WEEKS)
axes[2].plot(weeks, CA, label="CA", marker="o", markersize=3)
axes[2].plot(weeks, NY, label="NY", marker="s", markersize=3)
axes[2].set_xlabel("Week")
axes[2].set_ylabel("Price (\\$)")
axes[2].set_title("Common Trend")
axes[2].legend()

plt.tight_layout()
plt.show()

import numpy as np
from scipy import stats

# (1) 왼쪽 판
print(f"r = {r:.4f}   r^2 = {r**2:.4f}   기울기 = {np.polyfit(CA, NY, 1)[0]:.4f}")
print(f"CA {CA.min():.1f}~{CA.max():.1f}   NY {NY.min():.1f}~{NY.max():.1f}")

# (1) 가운데 판
prod = x_dev * y_dev
print(f"\n편차곱: 양수 {(prod > 0).sum()}개  음수 {(prod < 0).sum()}개")
print(f"  합/n = {prod.mean():.4f}   합/(n-1) = {prod.sum() / (WEEKS - 1):.4f} = Cov")
print(f"  최대 {prod.max():.2f}   최소 {prod.min():.2f}"
      f"   상위 5개가 차지하는 몫 {np.sort(prod)[-5:].sum() / prod.sum():.4f}")

# (1) 오른쪽 판 — 주당 기울기
t = np.arange(WEEKS)
bCA, bNY = np.polyfit(t, CA, 1)[0], np.polyfit(t, NY, 1)[0]
print(f"\n주당 기울기  CA {bCA:.4f}   NY {bNY:.4f}   비 {bNY / bCA:.4f}")

# (3) 추세를 빼고 다시 재기
res_CA = CA - np.polyval(np.polyfit(t, CA, 1), t)
res_NY = NY - np.polyval(np.polyfit(t, NY, 1), t)
rr, pp = stats.pearsonr(res_CA, res_NY)
print(f"\n추세 제거 전 r = {r:.4f}")
print(f"추세 제거 후 r = {rr:.4f}   p = {pp:.4f}"
      f"   SE ~ 1/sqrt(n-3) = {1 / np.sqrt(WEEKS - 3):.4f}")
print(f"잔차 표준편차 {res_CA.std(ddof=1):.4f} / {res_NY.std(ddof=1):.4f}"
      f"   (만들 때 쓴 값 0.5 / 0.6)")

r = 0.9753   r^2 = 0.9511   기울기 = 0.7816
CA 236.0~248.2   NY 339.5~350.2

편차곱: 양수 46개  음수 2개
  합/n = 10.3489   합/(n-1) = 10.5691 = Cov
  최대 31.86   최소 -0.50   상위 5개가 차지하는 몫 0.2915

주당 기울기  CA -0.2606   NY -0.2067   비 0.7932

추세 제거 전 r = 0.9753
추세 제거 후 r = 0.0320   p = 0.8292   SE ~ 1/sqrt(n-3) = 0.1491
잔차 표준편차 0.4553 / 0.5557   (만들 때 쓴 값 0.5 / 0.6)

import numpy as np

np.random.seed(42)
WEEKS = 48
trend_CA = np.linspace(0, -12, WEEKS)
trend_NY = np.linspace(0, -8, WEEKS)   # 기울기만 다른 추세

CA = 248.0 + trend_CA + np.random.normal(0, 3, WEEKS)
NY = 350.0 + trend_NY + np.random.normal(0, 3, WEEKS)

print(f"Covariance = {np.cov(CA, NY)[0, 1]:.4f}")
print(f"Pearson r  = {np.corrcoef(CA, NY)[0, 1]:.4f}")

Covariance = 10.5231
Pearson r  = 0.5733

np.random.seed(42)
CA = 248.0 + np.linspace(0, -12, WEEKS) + np.random.normal(0, 0.5, WEEKS)
NY = 350.0 + np.random.normal(0, 0.6, WEEKS)   # 추세 없음

print(f"Covariance = {np.cov(CA, NY)[0, 1]:.4f}")
print(f"Pearson r  = {np.corrcoef(CA, NY)[0, 1]:.4f}")

Covariance = 0.1348
Pearson r  = 0.0659

import numpy as np

rng = np.random.default_rng(25001)
print(f"{'평균 이동량':>12s} {'단순 공식':>16s} {'2단계 공식':>16s} {'상대오차':>12s}")
for shift in [0, 1e3, 1e6, 1e8, 1e9]:
    x = rng.standard_normal(1000) + shift
    y = rng.standard_normal(1000) * 2 + 0.7 * (x - shift) + shift
    n = len(x)
    naive = ((x * y).mean() - x.mean() * y.mean()) * n / (n - 1)
    two = ((x - x.mean()) * (y - y.mean())).sum() / (n - 1)
    print(f"{shift:12.0e} {naive:16.8f} {two:16.8f} "
          f"{abs(naive - two) / abs(two):12.2e}")

      평균 이동량            단순 공식           2단계 공식         상대오차
       0e+00       0.71369665       0.71369665     1.56e-16
       1e+03       0.76347205       0.76347205     6.90e-11
       1e+06       0.64957536       0.64952755     7.36e-05
       1e+08       0.00000000       0.63626123     1.00e+00
       1e+09       0.00000000       0.76178758     1.00e+00

print("\nfloat32 에서")
for shift in [0, 1e3, 1e5, 1e6]:
    x = (rng.standard_normal(1000) + shift).astype(np.float32)
    y = (rng.standard_normal(1000) * 2 + shift).astype(np.float32)
    naive = float((x * y).mean()) - float(x.mean()) * float(y.mean())
    two = float(((x - x.mean()) * (y - y.mean())).sum() / (len(x) - 1))
    true = float(np.cov(x.astype(np.float64), y.astype(np.float64))[0, 1])
    print(f"  shift={shift:8.0e}: 단순 {naive:14.6f}  2단계 {two:12.6f}  "
          f"참값 {true:10.6f}")

float32 에서
  shift=   0e+00: 단순       0.032305  2단계     0.032337  참값   0.032337
  shift=   1e+03: 단순       0.010370  2단계    -0.033220  참값  -0.033220
  shift=   1e+05: 단순       0.000549  2단계    -0.024096  참값  -0.024100
  shift=   1e+06: 단순   61440.000000  2단계    -0.031164  참값  -0.031155

def online_cov(xs, ys):
    """한 번의 통과로 안정적으로 공분산을 구한다 (웰퍼드)."""
    n = 0
    mx = my = C = 0.0
    for x, y in zip(xs, ys):
        n += 1
        dx = x - mx
        mx += dx / n
        my += (y - my) / n
        C += dx * (y - my)          # 갱신된 my 를 쓰는 것이 요점
    return C / (n - 1)

rng = np.random.default_rng(25001)
x = rng.standard_normal(10_000) + 1e8
y = rng.standard_normal(10_000) + 0.5 * (x - 1e8) + 1e8
print(f"\n  단순 공식: {(x * y).mean() - x.mean() * y.mean():.6f}")
print(f"  2단계:     {((x - x.mean()) * (y - y.mean())).sum() / (len(x) - 1):.6f}")
print(f"  웰퍼드:    {online_cov(x, y):.6f}")
print(f"  np.cov:    {np.cov(x, y)[0, 1]:.6f}")

  단순 공식: 0.000000
  2단계:     0.509465
  웰퍼드:    0.509465
  np.cov:    0.509465

import numpy as np

rng = np.random.default_rng(25002)
print(f"{'p (변수 수)':>11s} {'n (표본)':>9s} {'계수(rank)':>11s} "
      f"{'최소 고윳값':>13s} {'역행렬':>8s}")
for p, n in [(5, 100), (5, 10), (10, 10), (20, 15), (50, 30), (50, 60)]:
    X = rng.standard_normal((n, p))
    S = np.cov(X, rowvar=False)
    ev = np.linalg.eigvalsh(S)
    print(f"{p:11d} {n:9d} {np.linalg.matrix_rank(S):11d} {ev.min():13.2e} "
          f"{'가능' if ev.min() > 1e-10 else '불가':>8s}")

   p (변수 수)    n (표본)    계수(rank)        최소 고윳값      역행렬
          5       100           5      7.03e-01       가능
          5        10           5      2.28e-01       가능
         10        10           9     -2.11e-16       불가
         20        15          14     -3.68e-16       불가
         50        30          29     -1.91e-15       불가
         50        60          50      5.15e-03       가능

from sklearn.covariance import LedoitWolf
lw = LedoitWolf().fit(X)
lw.covariance_       → 축소된 공분산행렬 (항상 양정부호)
lw.precision_        → 그 역행렬
lw.shrinkage_        → 자동으로 고른 λ

import numpy as np

print("개별 분산 σ²=1, 모든 쌍 상관 ρ, 동일 비중 n 자산")
print(f"{'ρ':>6s} {'n=10':>9s} {'n=50':>9s} {'n=200':>9s} {'n→∞':>9s}")
for rho in [0.0, 0.1, 0.3, 0.5]:
    row = [(1 + (n - 1) * rho) / n for n in [10, 50, 200]]
    print(f"{rho:6.2f} {row[0]:9.4f} {row[1]:9.4f} {row[2]:9.4f} {rho:9.4f}")

rng = np.random.default_rng(25002)
print("\n모의실험으로 확인 (ρ=0.3)")
for n in [10, 50, 200]:
    S = np.full((n, n), 0.3)
    np.fill_diagonal(S, 1.0)
    L = np.linalg.cholesky(S)
    r = (rng.standard_normal((200_000, n)) @ L.T).mean(1)
    print(f"  n={n:3d}: 모의 분산 {r.var(ddof=1):.4f}, "
          f"이론 {(1 + (n - 1) * 0.3) / n:.4f}")

개별 분산 σ²=1, 모든 쌍 상관 ρ, 동일 비중 n 자산
     ρ      n=10      n=50     n=200       n→∞
  0.00    0.1000    0.0200    0.0050    0.0000
  0.10    0.1900    0.1180    0.1045    0.1000
  0.30    0.3700    0.3140    0.3035    0.3000
  0.50    0.5500    0.5100    0.5025    0.5000

모의실험으로 확인 (ρ=0.3)
  n= 10: 모의 분산 0.3696, 이론 0.3700
  n= 50: 모의 분산 0.3136, 이론 0.3140
  n=200: 모의 분산 0.3032, 이론 0.3035

import numpy as np

rng = np.random.default_rng(25002)
S = np.array([[4.0, 2.0], [2.0, 3.0]])
ev, V = np.linalg.eigh(S)
print(f"  공분산행렬 = {S.tolist()}")
print(f"  고윳값 = {ev[::-1].round(4).tolist()}   (합 = {ev.sum():.1f} = 대각합)")
print(f"  제1주성분 방향 = {V[:, -1].round(4).tolist()}")
print(f"  제1주성분이 설명하는 비율 = {ev[-1] / ev.sum():.4f}")

d = rng.multivariate_normal([0, 0], S, size=300_000)
print(f"  그 방향으로 투영한 분산 = {(d @ V[:, -1]).var(ddof=1):.4f}  "
      f"(고윳값 {ev[-1]:.4f})")

  공분산행렬 = [[4.0, 2.0], [2.0, 3.0]]
  고윳값 = [5.5616, 1.4384]   (합 = 7.0 = 대각합)
  제1주성분 방향 = [-0.7882, -0.6154]
  제1주성분이 설명하는 비율 = 0.7945
  그 방향으로 투영한 분산 = 5.5667  (고윳값 5.5616)

라이브러리마다 부호가 다를 수 있다
  → 주성분 점수의 부호가 뒤집힌다
  → 해석할 때는 "어느 변수와 같은 방향인지"를 기준으로 삼는다

1. E[XY] - E[X]E[Y] 를 그대로 쓰지 않는다
2. 자료를 두 번 읽을 수 있으면  → 2단계 공식
3. 한 번만 읽어야 하면          → 웰퍼드 온라인 갱신
4. 그냥 np.cov / pandas.cov 를 쓴다  ← 대부분의 경우 정답

np.var(x)             → n 으로 나눔   (ddof=0)
np.cov(x, y)          → n-1 로 나눔  (ddof=1)
np.cov(x, y)[0,0]  ≠  np.var(x)      ← 주의