회귀와 상관 그림¶

import numpy as np

np.random.seed(42)
DATA_SIZE = 100

# 앞 넷은 기울기를 바꾸고 잡음을 고정했다. 뒤 넷은 기울기를 3 으로 고정하고
# 잡음만 키웠다. 이렇게 나눠 놓으면 r 이 기울기가 아니라 "직선 둘레에
# 얼마나 몰려 있는가"를 재는 값임이 드러난다.
CONFIGS = [
    # (절편, 기울기, 잡음 크기)
    (2, 0.05, 1),    # 아주 완만한 기울기, 잡음 작음
    (2, -0.6, 1),    # 중간 크기 음의 기울기
    (2, 1.0,  1),    # 기울기 1
    (2, 3.0,  1),    # 가파른 기울기, 잡음 작음
    (2, 3.0,  3),    # 같은 기울기, 잡음 중간
    (2, 3.0, 10),    # 같은 기울기, 잡음 큼
    (2, 3.0, 20),    # 같은 기울기, 잡음 매우 큼
    (2, 3.0, 50),    # 같은 기울기, 잡음 극심
]


def generate(beta1, beta2, error_scale, n=DATA_SIZE):
    """주어진 기울기와 잡음으로 자료를 만들고 상관계수까지 돌려준다."""
    x = np.random.randint(1, n, n).astype(float)
    y = beta1 + beta2 * x + error_scale * np.random.randn(n)
    r = np.corrcoef(x, y)[0, 1]
    return x, y, r

import numpy as np

m = 99                       # randint(1, 100) 은 1..99 를 뽑는다
var_disc = (m**2 - 1) / 12   # 이산균등
var_cont = (100 - 1)**2 / 12 # 연속균등 Uniform(1, 100)
print(f"이산 var {var_disc:.4f}  sd {np.sqrt(var_disc):.4f}")
print(f"연속 var {var_cont:.4f}  sd {np.sqrt(var_cont):.4f}")
print(f"차 {var_cont - var_disc:.6f}  = 1/12 = {1 / 12:.6f}")

# 실제로 뽑아서 확인
rng = np.random.default_rng(0)
big = rng.integers(1, 100, 2_000_000)
print(f"200만 개 표본분산 {big.var():.4f}")

sx = np.sqrt(var_disc)
for b2, s in [(0.05, 1), (3.0, 1), (3.0, 50)]:
    snr = (b2 * sx / s) ** 2
    print(f"\nb2={b2}, sigma={s}:  b2*sx = {b2 * sx:.4f}"
          f"   SNR = {snr:.4f}   rho^2 = {snr / (1 + snr):.4f}"
          f"   rho = {np.sqrt(snr / (1 + snr)):.4f}")

# rho = 0.5 가 되는 sigma
print(f"\nb2=0.05 에서 rho=0.5 가 되는 sigma = {0.05 * sx * np.sqrt(3):.4f}")

이산 var 816.6667  sd 28.5774
연속 var 816.7500  sd 28.5788
차 0.083333  = 1/12 = 0.083333
200만 개 표본분산 817.0021

b2=0.05, sigma=1:  b2*sx = 1.4289   SNR = 2.0417   rho^2 = 0.6712   rho = 0.8193

b2=3.0, sigma=1:  b2*sx = 85.7321   SNR = 7350.0000   rho^2 = 0.9999   rho = 0.9999

b2=3.0, sigma=50:  b2*sx = 85.7321   SNR = 2.9400   rho^2 = 0.7462   rho = 0.8638

b2=0.05 에서 rho=0.5 가 되는 sigma = 2.4749

import matplotlib.pyplot as plt

# 여덟 칸을 한 화면에 놓고 r 값을 각 칸에 적어 둔다. 윗줄에서는 기울기가
# 아무리 달라도 r 이 비슷하고, 아랫줄에서는 기울기가 같은데 r 이 줄어든다.
fig, axes = plt.subplots(2, 4, figsize=(20, 9))

for idx, (b1, b2, es) in enumerate(CONFIGS):
    row, col = divmod(idx, 4)
    ax = axes[row, col]
    x, y, r = generate(b1, b2, es)

    ax.scatter(x, y, alpha=0.5, s=15, edgecolors="grey")
    xs = np.sort(x)
    ax.plot(xs, b1 + b2 * xs, color="#FA954D", lw=2, alpha=0.8)
    ax.set_title(f"Y = {b1} + {b2}X + {es}u", fontsize=10)
    ax.annotate(f"r = {r:.3f}", xy=(0.05, 0.9),
                xycoords="axes fraction", fontsize=11,
                bbox=dict(boxstyle="round", fc="wheat", alpha=0.5))

fig.suptitle("How Slope and Noise Affect Pearson Correlation",
             fontsize=14, y=1.01)
plt.tight_layout()
plt.show()

import numpy as np

sx2 = (99**2 - 1) / 12
sx = np.sqrt(sx2)

np.random.seed(42)
print(" b2    s    표본 r     이론 rho      차       SNR      y 폭   표본 s_x")
for b1, b2, es in CONFIGS:
    x = np.random.randint(1, DATA_SIZE, DATA_SIZE).astype(float)
    y = b1 + b2 * x + es * np.random.randn(DATA_SIZE)
    r = np.corrcoef(x, y)[0, 1]
    rho = b2 * sx / np.sqrt(b2**2 * sx2 + es**2)
    print(f"{b2:5.2f} {es:4d}  {r:+.4f}   {rho:+.4f}  {r - rho:+.4f}"
          f"  {(b2 * sx / es)**2:9.2f}  {y.max() - y.min():8.2f}"
          f"   {x.std(ddof=1):7.3f}")

# (3) X 를 좁은 구간에서 뽑으면
sx_small = np.sqrt((10**2 - 1) / 12)
snr = (0.05 * sx_small / 1) ** 2
print(f"\nX ~ {{1..10}} 이면 sx = {sx_small:.4f}  SNR = {snr:.4f}"
      f"  rho = {np.sqrt(snr / (1 + snr)):.4f}")

 b2    s    표본 r     이론 rho      차       SNR      y 폭   표본 s_x
 0.05    1  +0.8224   +0.8193  +0.0031       2.04      7.61    28.956
-0.60    1  -0.9979   -0.9983  +0.0004     294.00     58.56    28.146
 1.00    1  +0.9994   +0.9994  -0.0000     816.67     98.62    29.463
 3.00    1  +0.9999   +0.9999  +0.0000    7350.00    294.24    30.080
 3.00    3  +0.9992   +0.9994  -0.0002     816.67    287.48    27.123
 3.00   10  +0.9926   +0.9933  -0.0006      73.50    297.23    28.002
 3.00   20  +0.9774   +0.9739  +0.0035      18.38    331.80    29.662
 3.00   50  +0.8401   +0.8638  -0.0237       2.94    415.96    29.013

X ~ {1..10} 이면 sx = 2.8723  SNR = 0.0206  rho = 0.1422

import numpy as np
np.random.seed(42)
x = np.random.uniform(1, 100, 1000)
y = 2 + 3 * x + 10 * np.random.randn(1000)
print(f"Simulated r = {np.corrcoef(x, y)[0, 1]:.4f}")   # 0.9935
print("Theoretical rho = 0.9933")

Simulated r = 0.9935
Theoretical rho = 0.9933

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
slopes = [0.5, 2, 5]
noises = [1, 5, 20]

fig, axes = plt.subplots(3, 3, figsize=(12, 12))
for i, b2 in enumerate(slopes):
    for j, sigma in enumerate(noises):
        ax = axes[i, j]
        x = np.random.uniform(1, 50, 100)
        y = 1 + b2 * x + sigma * np.random.randn(100)
        r = np.corrcoef(x, y)[0, 1]

        ax.scatter(x, y, s=10, alpha=0.5)
        xs = np.sort(x)
        ax.plot(xs, 1 + b2 * xs, 'r-', lw=2)
        ax.set_title(f"b2={b2}, sigma={sigma}\nr={r:.3f}", fontsize=9)

plt.suptitle("Slope vs Noise: Effect on r", fontsize=14)
plt.tight_layout()
plt.show()

import numpy as np

rng = np.random.default_rng(27001)
print("같은 기울기 β=2, x 의 범위만 바꾼다 (σ=5 고정)")
print(f"{'x 범위':>10s} {'SD(x)':>8s} {'기울기 추정':>10s} {'r':>8s} {'r²':>8s}")
for lo, hi in [(-0.5, 0.5), (-2, 2), (-5, 5), (-20, 20)]:
    x = rng.uniform(lo, hi, 20_000)
    y = 2 * x + rng.normal(0, 5, 20_000)
    print(f"{f'[{lo},{hi}]':>10s} {x.std():8.3f} {np.polyfit(x, y, 1)[0]:10.4f} "
          f"{np.corrcoef(x, y)[0, 1]:8.4f} {np.corrcoef(x, y)[0, 1]**2:8.4f}")

같은 기울기 β=2, x 의 범위만 바꾼다 (σ=5 고정)
      x 범위    SD(x)     기울기 추정        r       r²
[-0.5,0.5]    0.289     1.9000   0.1090   0.0119
    [-2,2]    1.151     2.0032   0.4190   0.1755
    [-5,5]    2.888     1.9790   0.7503   0.5629
  [-20,20]   11.493     1.9970   0.9772   0.9549

print("\n같은 r≈0.6, 기울기만 바꾼다")
print(f"{'목표 β':>8s} {'σ_x':>7s} {'σ_ε':>8s} {'기울기 추정':>10s} {'r':>8s}")
for beta in [0.5, 2, 10, 50]:
    sx = 2.0
    se = np.sqrt(beta**2 * sx**2 * (1 / 0.36 - 1))   # r=0.6 이 되도록
    x = rng.normal(0, sx, 50_000)
    y = beta * x + rng.normal(0, se, 50_000)
    print(f"{beta:8.1f} {sx:7.1f} {se:8.2f} {np.polyfit(x, y, 1)[0]:10.4f} "
          f"{np.corrcoef(x, y)[0, 1]:8.4f}")

같은 r≈0.6, 기울기만 바꾼다
    목표 β     σ_x      σ_ε     기울기 추정        r
     0.5     2.0     1.33     0.5039   0.6043
     2.0     2.0     5.33     1.9985   0.6030
    10.0     2.0    26.67    10.0993   0.6050
    50.0     2.0   133.33    49.9082   0.5971

같은 n 으로 기울기를 더 정밀하게 추정하려면
  → x 를 양 극단에 몰아 배치한다 (D-최적 설계)
  → 단, 선형성을 확인할 수 없게 된다는 대가가 있다

import numpy as np

rng = np.random.default_rng(27001)
for rho in [0.3, 0.6, 0.9]:
    z = rng.standard_normal((100_000, 2))
    x = z[:, 0]
    y = rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]
    b_yx = np.polyfit(x, y, 1)[0]
    b_xy = np.polyfit(y, x, 1)[0]
    print(f"  ρ={rho}: b(y~x)={b_yx:.4f}, b(x~y)={b_xy:.4f}, "
          f"곱={b_yx * b_xy:.4f}, r²={np.corrcoef(x, y)[0, 1]**2:.4f}")

  ρ=0.3: b(y~x)=0.3001, b(x~y)=0.3025, 곱=0.0908, r²=0.0908
  ρ=0.6: b(y~x)=0.5993, b(x~y)=0.5994, 곱=0.3592, r²=0.3592
  ρ=0.9: b(y~x)=0.9006, b(x~y)=0.9018, 곱=0.8122, r²=0.8122

표준 회귀:    y 에만 오차가 있다고 가정
데밍 회귀:    x, y 의 오차 분산비를 알고 있을 때
직교회귀:     오차 분산이 같다고 가정 (총최소제곱)

import numpy as np

print(f"{'r':>6s} {'r²':>7s} {'잔차 SD / y 의 SD':>17s} {'폭의 감소':>10s}")
for r in [0.3, 0.5, 0.7, 0.9, 0.95, 0.99]:
    ratio = np.sqrt(1 - r**2)
    print(f"{r:6.2f} {r**2:7.3f} {ratio:17.4f} {1 - ratio:10.4f}")

     r      r²    잔차 SD / y 의 SD      폭의 감소
  0.30   0.090            0.9539     0.0461
  0.50   0.250            0.8660     0.1340
  0.70   0.490            0.7141     0.2859
  0.90   0.810            0.4359     0.5641
  0.95   0.902            0.3122     0.6878
  0.99   0.980            0.1411     0.8589

r = 0.7, 몸무게의 SD = 12 kg 이라 하자

  키를 모를 때:  예측구간 폭 ≈ 2 × 1.96 × 12   = 47.0 kg
  키를 알 때:    예측구간 폭 ≈ 2 × 1.96 × 8.57 = 33.6 kg

"상관이 0.7 로 강하다"고 해도
개별 예측은 ±17 kg 범위다

import warnings
warnings.filterwarnings("ignore")

import numpy as np

rng = np.random.default_rng(27002)
n = 200_000
sx, se = 1.0, 1.0
xt = rng.normal(0, sx, n)                # 참 x
yt = 2.0 * xt + rng.normal(0, se, n)     # 참 y
b0 = np.polyfit(xt, yt, 1)[0]
print(f"  오차 없음: 기울기 {b0:.4f}, r {np.corrcoef(xt, yt)[0, 1]:.4f}\n")
print(f"{'오차 위치':>12s} {'오차 SD':>8s} {'기울기':>9s} {'r':>8s} "
      f"{'신뢰도 λ':>9s} {'이론 기울기':>11s}")
for su in [0.5, 1.0, 2.0]:
    xo = xt + rng.normal(0, su, n)
    lam = sx**2 / (sx**2 + su**2)
    print(f"{'x 에만':>12s} {su:8.1f} {np.polyfit(xo, yt, 1)[0]:9.4f} "
          f"{np.corrcoef(xo, yt)[0, 1]:8.4f} {lam:9.4f} {b0 * lam:11.4f}")
for sv in [0.5, 1.0, 2.0]:
    yo = yt + rng.normal(0, sv, n)
    print(f"{'y 에만':>12s} {sv:8.1f} {np.polyfit(xt, yo, 1)[0]:9.4f} "
          f"{np.corrcoef(xt, yo)[0, 1]:8.4f} {'—':>9s} {b0:11.4f}")

  오차 없음: 기울기 2.0007, r 0.8944

       오차 위치    오차 SD       기울기        r     신뢰도 λ      이론 기울기
        x 에만      0.5    1.5971   0.7997    0.8000      1.6006
        x 에만      1.0    0.9980   0.6326    0.5000      1.0003
        x 에만      2.0    0.4005   0.4013    0.2000      0.4001
        y 에만      0.5    2.0003   0.8729         —      2.0007
        y 에만      1.0    2.0009   0.8154         —      2.0007
        y 에만      2.0    1.9963   0.6648         —      2.0007

y 의 오차:  y = βx + ε + v   →  오차항이 ε+v 로 커질 뿐, β 는 그대로
x 의 오차:  관측 x* = x + u  →  x* 와 오차항이 상관을 갖게 된다
                                → 내생성 → 편향

광고 노출 횟수(x)와 구매 확률 점수(y), n = 400

  기울기 b = 0.42 점/회,  95% CI [0.33, 0.51]
  r = 0.44,  r² = 0.19
  잔차 SD = 8.6 점  (y 의 SD 는 9.6 점)

  x 의 범위는 0~50 회였다. 범위가 다른 연구와 r 을 직접
  비교하지 않았고, 기울기로 비교했다.

  노출 횟수는 로그에 기록되어 측정오차가 작다고 보았다.