주성분회귀와 부분최소제곱 보기¶

import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score

url = ("https://raw.githubusercontent.com/gedeck/"
       "practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/house_sales.csv")
house = pd.read_csv(url, sep='\t')

# 수치형 변수만 쓴다. 주성분은 분산을 기준으로 방향을 찾으므로
# 범주형 가변수를 섞으면 뜻이 흐려진다.
numeric_features = [
    'SqFtTotLiving', 'SqFtLot', 'Bathrooms', 'Bedrooms',
    'BldgGrade', 'NbrLivingUnits', 'SqFtFinBasement', 'YrBuilt', 'YrRenovated'
]
X = house[numeric_features].values
y = house['AdjSalePrice'].values

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

ols_model = LinearRegression().fit(X_scaled, y)
ols_r2 = r2_score(y, ols_model.predict(X_scaled))
ols_rmse = np.sqrt(mean_squared_error(y, ols_model.predict(X_scaled)))

print(f"n = {len(y)}, p = {X_scaled.shape[1]}")
print(f"OLS: R^2 = {ols_r2:.6f},  RMSE = {ols_rmse:,.1f}")
d = np.linalg.svd(X_scaled, compute_uv=False)
lam = d ** 2 / (len(y) - 1)
print(f"상관행렬 고윳값의 합 = {lam.sum():.6f}  (p = {X_scaled.shape[1]})")
print(f"고윳값: {np.round(lam, 4)}")

n = 22687, p = 9
OLS: R^2 = 0.593225,  RMSE = 245,800.5
상관행렬 고윳값의 합 = 9.000397  (p = 9)
고윳값: [3.2464 1.4904 1.0379 0.9616 0.87   0.5566 0.4485 0.25   0.1388]

from sklearn.decomposition import PCA
from sklearn.model_selection import cross_val_score, KFold

# 주성분회귀(PCR): 먼저 주성분을 뽑고 그중 앞의 M 개로 회귀한다.
# 주성분은 반응 y 를 전혀 보지 않고 X 의 분산만 보고 정해진다는 점이
# 아래 PLS 와 갈리는 지점이다.
pca = PCA()
X_pca = pca.fit_transform(X_scaled)

explained_var = pca.explained_variance_ratio_
cumsum_var = np.cumsum(explained_var)

kfold = KFold(n_splits=10, shuffle=True, random_state=42)
pcr_mse_scores = []

# 몇 개의 성분을 쓸지는 교차검증으로 고른다.
for M in range(1, X_scaled.shape[1] + 1):
    reg = LinearRegression()
    cv_scores = cross_val_score(
        reg, X_pca[:, :M], y,
        cv=kfold, scoring='neg_mean_squared_error'
    )
    pcr_mse_scores.append(-cv_scores.mean())

M_opt_pcr = np.argmin(pcr_mse_scores) + 1
pcr_cv_rmse = np.sqrt(pcr_mse_scores[M_opt_pcr - 1])

print("성분별 설명분산:", np.round(explained_var, 4))
print("누적 설명분산:  ", np.round(cumsum_var, 4))
print("PCR CV RMSE:", [f"{np.sqrt(m):,.0f}" for m in pcr_mse_scores])
print(f"M_opt(PCR) = {M_opt_pcr},  CV RMSE = {pcr_cv_rmse:,.1f}")

# 모든 성분을 쓰면 OLS 와 같아야 한다 (보기 1 의 (2))
r2_full = r2_score(y, LinearRegression().fit(X_pca, y).predict(X_pca))
print(f"모든 성분을 쓴 PCR 의 R^2 = {r2_full:.10f},  OLS = {ols_r2:.10f},  "
      f"차이 {abs(r2_full - ols_r2):.2e}")

성분별 설명분산: [0.3607 0.1656 0.1153 0.1068 0.0967 0.0618 0.0498 0.0278 0.0154]
누적 설명분산:   [0.3607 0.5263 0.6416 0.7485 0.8451 0.907  0.9568 0.9846 1.    ]
PCR CV RMSE: ['296,154', '292,911', '287,749', '287,519', '287,529', '275,200', '246,949', '246,740', '246,113']
M_opt(PCR) = 9,  CV RMSE = 246,113.3
모든 성분을 쓴 PCR 의 R^2 = 0.5932245909,  OLS = 0.5932245909,  차이 1.11e-16

from sklearn.cross_decomposition import PLSRegression

# 부분최소제곱(PLS): 성분을 찾을 때 y 와의 공분산까지 함께 본다.
# 그래서 같은 성분 수라면 대개 PCR 보다 낫지만, y 를 보고 방향을 정한
# 만큼 과적합의 여지도 생긴다.
pls_mse_scores = []
for M in range(1, X_scaled.shape[1] + 1):
    pls = PLSRegression(n_components=M)
    cv_scores = cross_val_score(
        pls, X_scaled, y,
        cv=kfold, scoring='neg_mean_squared_error'
    )
    pls_mse_scores.append(-cv_scores.mean())

M_opt_pls = np.argmin(pls_mse_scores) + 1
pls_cv_rmse = np.sqrt(pls_mse_scores[M_opt_pls - 1])

# (1) 의 첫 가중벡터를 직접 계산해 본다
w1 = X_scaled.T @ (y - y.mean())
w1 = w1 / np.linalg.norm(w1)
print("PLS 첫 성분의 가중 w1:")
for nm, v in sorted(zip(numeric_features, w1), key=lambda t: -abs(t[1])):
    print(f"  {nm:16s} {v:+.4f}")
print("PLS CV RMSE:", [f"{np.sqrt(m):,.0f}" for m in pls_mse_scores])
print(f"M_opt(PLS) = {M_opt_pls},  CV RMSE = {pls_cv_rmse:,.1f}")

PLS 첫 성분의 가중 w1:
  SqFtTotLiving    +0.5786
  BldgGrade        +0.5623
  Bathrooms        +0.4404
  Bedrooms         +0.2598
  SqFtFinBasement  +0.2478
  SqFtLot          +0.1142
  YrRenovated      +0.0897
  YrBuilt          +0.0684
  NbrLivingUnits   +0.0188
PLS CV RMSE: ['280,028', '251,752', '246,747', '246,317', '246,167', '246,123', '246,116', '246,113', '246,113']
M_opt(PLS) = 8,  CV RMSE = 246,113.2

from sklearn.linear_model import RidgeCV

# 능형회귀와도 견준다. 셋 다 상관된 변수를 다루는 방법이지만, 능형은
# 변수를 그대로 두고 계수만 줄이는 반면 PCR·PLS 는 변수를 적은 수의
# 성분으로 갈아 끼운다.
ridge_cv = RidgeCV(alphas=np.logspace(-2, 5, 100), cv=10)
ridge_cv.fit(X_scaled, y)
ridge_r2 = r2_score(y, ridge_cv.predict(X_scaled))
ridge_rmse = np.sqrt(mean_squared_error(y, ridge_cv.predict(X_scaled)))

print(f"능형: alpha = {ridge_cv.alpha_:,.2f},  R^2 = {ridge_r2:.6f},  RMSE = {ridge_rmse:,.1f}")
print(f"유효자유도 = {np.sum(d**2 / (d**2 + ridge_cv.alpha_)):.4f}  (p = 9)")
print(f"OLS RMSE {ols_rmse:,.1f} / 능형 {ridge_rmse:,.1f}  -> 차이 {ridge_rmse - ols_rmse:,.1f}")

능형: alpha = 335.16,  R^2 = 0.593056,  RMSE = 245,851.5
유효자유도 = 8.7301  (p = 9)
OLS RMSE 245,800.5 / 능형 245,851.5  -> 차이 51.0

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

np.random.seed(42)
n, p = 50, 5
X = np.random.randn(n, p)
beta_true = np.array([3, -1, 2, 0, 0])
y = X @ beta_true + np.random.randn(n) * 0.5

# 표준화. 벌점회귀에서는 필수다
X_s = StandardScaler().fit_transform(X)

# 특이값분해로 직접 구현한 주성분회귀
U, D, Vt = np.linalg.svd(X_s, full_matrices=False)
M = 3  # number of components
Z = U[:, :M] * D[:M]  # first M principal components
y_bar = y.mean()                                  # intercept
gamma = np.linalg.lstsq(Z, y - y_bar, rcond=None)[0]
beta_pcr_manual = Vt[:M].T @ gamma
y_pred_manual = X_s @ beta_pcr_manual + y_bar

# sklearn 으로 같은 일 하기
pca = PCA(n_components=M)
Z_sk = pca.fit_transform(X_s)
reg = LinearRegression().fit(Z_sk, y)
y_pred_sk = reg.predict(Z_sk)

print(f"Max prediction difference: {np.max(np.abs(y_pred_manual - y_pred_sk)):.2e}")

Max prediction difference: 9.77e-15