잔차제곱합 곡면 시각화¶

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler

# 광고비와 매출을 흉내 낸 자료. 참 기울기는 0.05 다.
np.random.seed(42)
n_samples = 100
TV = np.random.uniform(0, 300, n_samples)
Sales = 7 + 0.05 * TV + np.random.normal(0, 2, n_samples)

# 중심화만 하고 척도는 건드리지 않는다(with_std=False). 중심화하면 절편과
# 기울기의 추정이 서로 독립이 되어, 아래 등고선이 기울어지지 않고 바로 선다.
X = TV.reshape(-1, 1)
scaler = StandardScaler(with_mean=True, with_std=False)
X_scaled = scaler.fit_transform(X)

model = LinearRegression()
model.fit(X_scaled, Sales)
beta_0 = model.intercept_
beta_1 = model.coef_[0]

# 중심화의 두 효과를 확인한다.
print(f"X_scaled 의 합 = {X_scaled.sum():.3e}   (중심화했으므로 0)")
print(f"beta_0 = {beta_0:.6f},   Sales 평균 = {Sales.mean():.6f},"
      f"   차이 = {abs(beta_0 - Sales.mean()):.2e}")
print(f"beta_1 = {beta_1:.8f}   (참 기울기 0.05)")
print()
# 중심화 좌표에서 참 절편은 7 + 0.05*TV-bar 다
print(f"TV 평균 = {TV.mean():.4f}")
print(f"중심화 좌표에서의 참 절편 = 7 + 0.05*{TV.mean():.4f} = {7 + 0.05 * TV.mean():.6f}")
print()
S_xx = (X_scaled.ravel() ** 2).sum()
residuals = Sales - model.predict(X_scaled)
rss_min = (residuals ** 2).sum()
s_square = rss_min / (n_samples - 2)
print(f"S_xx = {S_xx:.4f}   (이론 (n-1)*300^2/12 = {99 * 7500:.1f})")
print(f"RSS_min = {rss_min:.4f},  s^2 = {s_square:.6f},  s = {np.sqrt(s_square):.6f}  (참 sigma 2)")
print(f"se(beta_0) = s/sqrt(n)    = {np.sqrt(s_square / n_samples):.6f}")
print(f"se(beta_1) = s/sqrt(S_xx) = {np.sqrt(s_square / S_xx):.8f}")
print(f"(beta_1 - 0.05)/se = {(beta_1 - 0.05) / np.sqrt(s_square / S_xx):+.4f}")
print(f"(beta_0 - {7 + 0.05 * TV.mean():.4f})/se = "
      f"{(beta_0 - (7 + 0.05 * TV.mean())) / np.sqrt(s_square / n_samples):+.4f}")

X_scaled 의 합 = -5.400e-13   (중심화했으므로 0)
beta_0 = 14.050550,   Sales 평균 = 14.050550,   차이 = 0.00e+00
beta_1 = 0.04693485   (참 기울기 0.05)

TV 평균 = 141.0542
중심화 좌표에서의 참 절편 = 7 + 0.05*141.0542 = 14.052711

S_xx = 788534.5515   (이론 (n-1)*300^2/12 = 742500.0)
RSS_min = 322.6338,  s^2 = 3.292182,  s = 1.814437  (참 sigma 2)
se(beta_0) = s/sqrt(n)    = 0.181444
se(beta_1) = s/sqrt(S_xx) = 0.00204330
(beta_1 - 0.05)/se = -1.5001
(beta_0 - 14.0527)/se = -0.0119

# 최적해 둘레로 격자를 깔고 칸마다 잔차제곱합을 계산한다.
# 최소제곱이 무엇을 최소화하는지를 눈으로 보려는 것이다.
B0_range = np.linspace(beta_0 - 2, beta_0 + 2, 50)
B1_range = np.linspace(beta_1 - 0.05, beta_1 + 0.05, 50)
B0_mesh, B1_mesh = np.meshgrid(B0_range, B1_range)

RSS = np.zeros_like(B0_mesh)
for i in range(B0_mesh.shape[0]):
    for j in range(B0_mesh.shape[1]):
        # ravel() 이 반드시 있어야 한다. X_scaled 는 (100, 1) 이고 Sales 는
        # (100,) 이므로, 그대로 빼면 (100, 100) 으로 퍼져 10000 개를 더한다.
        # 보기 3 에서 그 덫을 다룬다.
        y_pred = B0_mesh[i, j] + B1_mesh[i, j] * X_scaled.ravel()
        RSS[i, j] = np.sum((Sales - y_pred) ** 2)

# RSS 의 이차 전개가 정확한지, 격자의 네 자리에서 확인한다.
def rss_at(b0, b1):
    return np.sum((Sales - (b0 + b1 * X_scaled.ravel())) ** 2)

print(f"RSS_min = {rss_min:.4f},  n = {n_samples},  S_xx = {S_xx:.4f}")
print()
print(f"{'(d0, d1)':>18s}{'실제 RSS':>16s}{'RSS_min + n d0^2 + S_xx d1^2':>32s}{'차이':>12s}")
for d0, d1 in [(0, 0), (2, 0), (0, 0.05), (2, 0.05), (-2, -0.05), (1, -0.025)]:
    actual = rss_at(beta_0 + d0, beta_1 + d1)
    predicted = rss_min + n_samples * d0 ** 2 + S_xx * d1 ** 2
    print(f"{f'({d0}, {d1})':>18s}{actual:16.4f}{predicted:32.4f}{abs(actual - predicted):12.2e}")
print()
print(f"격자 끝에서 두 항의 크기")
print(f"  n * 2^2        = {n_samples * 4:.1f}")
print(f"  S_xx * 0.05^2  = {S_xx * 0.0025:.1f}")
print(f"  비 = {S_xx * 0.0025 / (n_samples * 4):.3f}")
print(f"RSS 가 격자 안에서 커지는 최대 배수 = "
      f"{(rss_min + n_samples * 4 + S_xx * 0.0025) / rss_min:.3f}")

RSS_min = 322.6338,  n = 100,  S_xx = 788534.5515

          (d0, d1)          실제 RSS    RSS_min + n d0^2 + S_xx d1^2          차이
            (0, 0)        322.6338                        322.6338    0.00e+00
            (2, 0)        722.6338                        722.6338    3.41e-13
         (0, 0.05)       2293.9702                       2293.9702    1.82e-12
         (2, 0.05)       2693.9702                       2693.9702    1.82e-12
       (-2, -0.05)       2693.9702                       2693.9702    3.18e-12
       (1, -0.025)        915.4679                        915.4679    1.36e-12

격자 끝에서 두 항의 크기
  n * 2^2        = 400.0
  S_xx * 0.05^2  = 1971.3
  비 = 4.928
RSS 가 격자 안에서 커지는 최대 배수 = 8.350

(100,) - (100, 1)  ->  (100, 100)

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(16, 6))

# 왼쪽: 등고선. 별표가 최소점이고, 그것이 곧 최소제곱추정값이다.
ax1 = fig.add_subplot(121)
contour = ax1.contour(B0_mesh, B1_mesh, RSS / 1000, levels=20, cmap='viridis')
ax1.plot(beta_0, beta_1, 'r*', markersize=20, label='Optimal')
ax1.set_xlabel('beta_0 (Intercept)')
ax1.set_ylabel('beta_1 (Slope)')
ax1.set_title('RSS Contour Plot')

# 오른쪽: 같은 것을 곡면으로. RSS 가 계수의 이차함수이므로 사발 모양이고,
# 그래서 최소점이 하나뿐이며 닫힌 해가 존재한다.
ax2 = fig.add_subplot(122, projection='3d')
ax2.plot_surface(B0_mesh, B1_mesh, RSS / 1000, cmap='viridis', alpha=0.8)
ax2.set_xlabel('beta_0')
ax2.set_ylabel('beta_1')
ax2.set_zlabel('RSS / 1000')
ax2.set_title('RSS 3D Surface')

plt.tight_layout()
plt.show()

# (1) 등고선의 최소점이 정말 최소제곱해인지 확인한다. 곡면 그림을 그릴 때
#     반드시 넣어야 하는 한 줄이다.
i_min, j_min = np.unravel_index(np.argmin(RSS), RSS.shape)
print(f"RSS 배열의 최솟값 자리 = (행 {i_min}, 열 {j_min})   격자 가운데는 (24~25, 24~25)")
print(f"  그 자리의 beta_0 = {B0_mesh[i_min, j_min]:.6f}   (최소제곱해 {beta_0:.6f})")
print(f"  그 자리의 beta_1 = {B1_mesh[i_min, j_min]:.8f}   (최소제곱해 {beta_1:.8f})")
print(f"RSS 배열의 최솟값 = {RSS.min():.4f}   (참 최소제곱 RSS {rss_min:.4f})")

# 격자가 참 최솟값을 얼마나 놓칠 수 있는지는 보기 2 의 이차식이 말해 준다.
g0 = (B0_range[1] - B0_range[0]) / 2
g1 = (B1_range[1] - B1_range[0]) / 2
print(f"  격자 반칸 = ({g0:.4f}, {g1:.5f}),  그래서 최대 "
      f"{n_samples * g0 ** 2 + S_xx * g1 ** 2:.2f} 만큼 뜰 수 있다")
print(f"  실제로 뜬 값 = {RSS.min() - rss_min:.2f}")
print()

# (2) ravel() 을 빼면 무엇이 되는가. 모양이 (100,1) 과 (100,) 이라 퍼진다.
y_pred_demo = beta_0 + beta_1 * X_scaled          # (100, 1)
print(f"Sales.shape = {Sales.shape},  (beta_0 + beta_1*X_scaled).shape = {y_pred_demo.shape}")
print(f"  (Sales - y_pred).shape = {(Sales - y_pred_demo).shape}   <- 10000 개를 더한다")
print(f"  ravel() 을 쓰면 {(Sales - y_pred_demo.ravel()).shape}   <- 100 개")
print()
S_yy = ((Sales - Sales.mean()) ** 2).sum()
print("퍼진 합은 n*[S_yy + n(b0-ybar)^2 + b1^2 S_xx] 와 같다")
for d0, d1 in [(0, 0), (2, 0), (0, 0.05)]:
    actual = np.sum((Sales - ((beta_0 + d0) + (beta_1 + d1) * X_scaled)) ** 2)
    predicted = n_samples * (S_yy + n_samples * d0 ** 2 + (beta_1 + d1) ** 2 * S_xx)
    print(f"  ({d0}, {d1})  실제 {actual:14.4f}  예측 {predicted:14.4f}  차이 {abs(actual - predicted):.2e}")
print()

# ravel() 을 뺀 격자를 만들어 최소점이 어디로 가는지 본다.
RSS_broken = np.zeros_like(B0_mesh)
for i in range(B0_mesh.shape[0]):
    for j in range(B0_mesh.shape[1]):
        RSS_broken[i, j] = np.sum((Sales - (B0_mesh[i, j] + B1_mesh[i, j] * X_scaled)) ** 2)
i2, j2 = np.unravel_index(np.argmin(RSS_broken), RSS_broken.shape)
print(f"ravel() 을 빼면 최솟값 자리 = (행 {i2}, 열 {j2}),  값 {RSS_broken.min():.4f}")
print(f"  그 자리의 beta_1 = {B1_mesh[i2, j2]:.8f}  (최소제곱해 {beta_1:.8f})")
print(f"  B1 격자에서 0 에 가장 가까운 칸 = {int(np.argmin(np.abs(B1_range)))}")
print(f"  참 최솟값의 {RSS_broken.min() / rss_min:.1f} 배")
print(f"  RSS/1000 의 범위:  옳은 판 {RSS.min() / 1000:.3f} ~ {RSS.max() / 1000:.3f},"
      f"   ravel() 뺀 판 {RSS_broken.min() / 1000:.3f} ~ {RSS_broken.max() / 1000:.3f}")

RSS 배열의 최솟값 자리 = (행 24, 열 24)   격자 가운데는 (24~25, 24~25)
  그 자리의 beta_0 = 14.009734   (최소제곱해 14.050550)
  그 자리의 beta_1 = 0.04591444   (최소제곱해 0.04693485)
RSS 배열의 최솟값 = 323.6215   (참 최소제곱 RSS 322.6338)
  격자 반칸 = (0.0408, 0.00102),  그래서 최대 0.99 만큼 뜰 수 있다
  실제로 뜬 값 = 0.99

Sales.shape = (100,),  (beta_0 + beta_1*X_scaled).shape = (100, 1)
  (Sales - y_pred).shape = (100, 100)   <- 10000 개를 더한다
  ravel() 을 쓰면 (100,)   <- 100 개

퍼진 합은 n*[S_yy + n(b0-ybar)^2 + b1^2 S_xx] 와 같다
  (0, 0)  실제    379672.7322  예측    379672.7322  차이 0.00e+00
  (2, 0)  실제    419672.7322  예측    419672.7322  차이 0.00e+00
  (0, 0.05)  실제    946903.8408  예측    946903.8408  차이 1.16e-10

ravel() 을 빼면 최솟값 자리 = (행 2, 열 25),  값 206066.1906
  그 자리의 beta_1 = 0.00101648  (최소제곱해 0.04693485)
  B1 격자에서 0 에 가장 가까운 칸 = 2
  참 최솟값의 638.7 배
  RSS/1000 의 범위:  옳은 판 0.324 ~ 2.694,   ravel() 뺀 판 206.066 ~ 986.904

from sklearn.preprocessing import StandardScaler

Xz = StandardScaler().fit_transform(TV.reshape(-1, 1)).flatten()

beta = np.array([0.0, 0.0])
lr = 0.1
for step in range(200):
    residuals = Sales - beta[0] - beta[1] * Xz
    grad = np.array([-2 * residuals.mean(),
                     -2 * (residuals * Xz).mean()])
    beta -= lr * grad
print(f"GD solution: beta_0={beta[0]:.6f}, beta_1={beta[1]:.6f}")

GD solution: beta_0=14.050550, beta_1=4.167789