Numpy로 직접 구현하는 Linear Regression

Regression

• 회귀 분석: 관찰된 연속형 변수들에 대해 두 변수 사이의 모형을 구한 뒤 적합도를 측정하는 분석 방법
• 변수: 값이 변하는 데이터 요소 또는 속성
  • 독립변수(설명변수) : 원인(x, input)
  • 종속변수(응답변수): 결과(y, output)

Linear Regression

• 독립 변수 하나 이상 + 종속 변수 하나 의 관계를 모델링하는 통계적 방법
• $y=mx+b$
• 가정
  • 선형성: 종속 변수와 독립 변수간의 관계가 선형적
  • 독립성: 관측값들은 서로 독립적 (잔차(residual, 오차의 추정치)들이 무작위 분포)
  • 등분산성: 오류의 분산이 일정 (잔차들이 일정한 분포를 보임)
  • 정규성: 오류가 정규분포를 따름

최소 제곱법(OLS) 방법

• 예측된 값 $\hat{y} = Xw$와 실제 값 y 사이의 차이(잔차)를 최소화하는 가중치 벡터 $w$를 찾는 closed-form solution이다.
• 선형회귀모델의 파라미터를 추정하는 방법 중 하나로, 다른 방식으로는 WLS(weighted least squares), GLS(generalized least squares), LASSO, Ridge 등이 있다.
  $$Cost(mb) = \sum(y_i - (mx_i +b))^2$$

$$L(w) = ||y-Xw||^2 $$
$$\frac{\partial L(w)}{\partial w} = -2X^T(y - Xw)$$
$$ X^T X w = X^T y$$
$$w = \left(X^T X\right)^{-1} X^T y$$
이렇게 마지막으로 구한 식은 정규방정식(Normal Equation)이라 불리며, 다음과 같은 형태로 잘 알려져 있다.
$$ \beta = (X^TX)^{-1}X^Ty $$ numpy 코드로 구현하면 다음과 같다.

import numpy as np
beta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)

모델 평가 지표

• 펑균 절대 오차(MAE)
  • 해석이 간단하고 단위가 종속 변수와 동일
    $$MAE = \frac{1}{n}\sum_{i=1}^{n}(y_i-\hat{y}_i)$$
• 평균 제곱 오차(MSE)
  • 제곱하므로 큰 오차에 더 큰 패널티 부여하여 큰 오차를 줄이는데 유용
    $$MSE = \frac{1}{n}\sum_{i=1}^{n}(y_i-\hat{y}_i)^2$$
• 제곱근 평균 제곱 오차(RMSE)
  • 해석이 용이하고 큰 오차에 민감하게 반응함과 동시에 실제 값과 같은 단위
    $$RMSE = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_i-\hat{y}_i)^2}$$
• 결정계수($R^2$, Coefficient of Determination)
  • 모델이 종속 변수의 변동성을 얼마나 설명하는지 나타냄
  • 0~1 사이 값을 가지며 1에 가까울수록 x와 y는 완변한 선형 상관관계를 이룸

$$ R^2 = \frac{SSR}{SST} = 1 - \frac{SSE}{SST}$$

(총 변동 SST = 설명 가능한 변동 SSR + 설명 불가능한 변동 SSE)

$$R^2 = 1 - \frac{\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}{\sum_{i=1}^{n}(y_i - \bar{y}_i)^2}$$

Numpy로 Linear Regression 코드 구현

Iris Data 불러오기

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
import numpy as np

iris = load_iris() # iris 데이터
X = iris.data[:, [0]] # 특성 데이터에서 첫 번째 특성만 사용
y = iris.data[:, 2] # 타겟 데이터로 세 번째 특성 사용

# 학습-테스트 데이터 분리
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

X_train = np.array(X_train).reshape(-1, 1)
y_train = np.array(y_train)
X_test = np.array(X_test).reshape(-1, 1)
y_test = np.array(y_test)

print(X_train.shape, y_train.shape)
print(X_test.shape, y_test.shape)

>>> (120, 1) (120,)
>>> (30, 1) (30,)

Linear Regression

class LinearRegression:
    def __init__(self):
        self.weights = None
        self.bias = None

    def fit(self, X, y):
        # X에 절편(bias) 추가
        X = np.hstack((np.ones((X.shape[0], 1)), X))
        # 최소제곱법(OLS)로 구한 정규 방정식(Normal Equation)
        beta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
        self.weights = beta[0]
        self.bias = beta[1]

    def predict(self, X):
        # X에 절편 추가
        X = np.hstack((np.ones((X.shape[0], 1)), X))
        # 학습한 가중치와 바이어스로 y = Wx+b 예측값 반환
        y_pred = X.dot(np.hstack((self.weights, self.bias)))
        return y_pred

model = LinearRegression()
model.fit(X_train, y_train)
predictions = model.predict(X_test)

평가지표

def mean_squared_error(y_true, y_pred):
    squared_errors = np.square(y_true - y_pred)
    mse = np.mean(squared_errors)

    return mse

def r_squared(y_true, y_pred):
    mean_y_true = np.sum(y_true) / len(y_true)
    tss = np.sum(np.square(y_true - mean_y_true))
    rss = np.sum(np.square(y_true - y_pred))
    r2 = 1 - (rss / tss)

    return r2

# 직접 구현
MSE = mean_squared_error(y_test, predictions)
r2 = r_squared(y_test, predictions)
MSE, r2

>>> (0.5960765879745361, 0.8181245472591383)

from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score

# sklearn metric 함수와 비교
MSE_sklearn = mean_squared_error(y_test, predictions)
r2_sklearn = r2_score(y_test, predictions)
MSE_sklearn,r2_sklearn

>>> (0.5960765879745361, 0.8181245472591383)

시각화

import matplotlib.pyplot as plt

plt.scatter(X_test, y_test, label='Data')
plt.plot(X_test, predictions, color="r", label="Regression Line")

for i in range(len(X_test)):
    plt.plot([X_test[i], X_test[i]], [y_test[i], predictions[i]], color='gray', linestyle='--', linewidth=0.5)

plt.title('Iris Linear Regression')
plt.legend()
plt.show()