[머신러닝] Multiple Linear Regression - 연비(MPG) 예측

2022.07.05 - [Studying/Machine Learning] - [머신러닝] python 지구 온도 변화 분석

[머신러닝] Linear Regression - 지구 온도 변화 분석

저번 포스트가 feature가 하나인 single linear regression 이였으니

이번에는 multiple linear regression을 해보겠다.

 

여러 feature가 나와서 처리하게 되는 프로세스들이 추가된다.

---

이번에도 데이터를 먼저 다운로드해보겠다.

 

import pandas

import seaborn	#시각화를 위해 사용
seaborn.set()

from urllib.request import urlretrieve

URL = 'https://go.gwu.edu/engcomp6data3'
urlretrieve(URL, 'auto_mpg.csv')

받은 파일을 열어서 데이터를 확인해보겠다.

 

mpg_data = pandas.read_csv('/content/auto_mpg.csv')
mpg_data.info()

# <class 'pandas.core.frame.DataFrame'>
# RangeIndex: 392 entries, 0 to 391
# Data columns (total 9 columns):
#  #   Column        Non-Null Count  Dtype  
# ---  ------        --------------  -----  
#  0   mpg           392 non-null    float64
#  1   cylinders     392 non-null    int64  
#  2   displacement  392 non-null    float64
#  3   horsepower    392 non-null    float64
#  4   weight        392 non-null    float64
#  5   acceleration  392 non-null    float64
#  6   model year    392 non-null    int64  
#  7   origin        392 non-null    int64  
#  8   car name      392 non-null    object 
# dtypes: float64(5), int64(3), object(1)
# memory usage: 27.7+ KB

총 392개의 데이터가 있고 9개의 정보들이 있다.

 

이 데이터 중에 연비에 영향이 가지 않는 origin, car_name은 제외하고 linear regression을 해보려고 한다.

 

y_col = 'mpg'
x_cols = mpg_data.columns.drop(['car name', 'origin', 'mpg'])

print(x_cols)
# Index(['cylinders', 'displacement', 'horsepower', 'weight', 'acceleration', 'model year'],
#       dtype='object')

car_name과 origin과 예측할 mpg 컬럼을 drop 시킨 결과이다.

 

---

linear regression을 하기 전에 자동차의 스펙들과 연비의 상관관계를 시각화해서 알아보자

 

seaborn.pairplot(data=mpg_data, height=5, aspect=1,
             x_vars=x_cols,
             y_vars=y_col);

데이터는 이렇게 보인다.

 

acceleration과 model year과 양의 상관관계, 나머지는 음의 상관관계에 있는 것을 볼 수 있다.

 

이러한 상관관계이기 때문에 연비를 예측하는데 linear regression이 충분하다는 것을 알 수 있다.

 

---

입력변수가 d개일 때 우리가 찾아야 하는 모델은 다음과 같다.

 

$$\hat{y} = w_0 + w_1x_1 + w_2x_2 + ... + w_dx_d$$

 

여기서 y hat은 우리가 예측한 값이다.

 

여기서 연산의 편의 상 x_0=1을 추가해서 위의 식을 정리하면 다음과 같다.

 

$$\hat{y} = \sum_{i=0}^d w_ix_i = \mathbf{x}^\top \mathbf{w}$$

 

위에서 392개의 데이터를 받았는데 이를 N개라 하고 식을 나열해보면 다음과 같다.

 

 

$$\hat{y}^{(1)} = w_0x_0^{(1)} + w_1x_1^{(1)} + w_2x_2^{(1)} + ... + w_dx_d^{(1)}$$ $$\hat{y}^{(2)} = w_0x_0^{(2)} + w_1x_1^{(2)} + w_2x_2^{(2)} + ... + w_dx_d^{(2)}$$ $$\vdots$$ $$\hat{y}^{(N)} = w_0x_0^{(N)} + w_1x_1^{(N)} + w_2x_2^{(N)} + ... + w_dx_d^{(N)}$$

 

위의 식을 한 번에 행렬로 표현하면 다음과 같다.

 

$$\hat{\mathbf{y}} = \begin{bmatrix} \hat{y}^{(1)} \\ \hat{y}^{(2)}\\ \vdots \\ \hat{y}^{(N)} \end{bmatrix} = \begin{bmatrix} x_0^{(1)} & x_1^{(1)} & \cdots & x_d^{(1)} \\ x_0^{(2)} & x_1^{(2)} & \cdots & x_d^{(2)} \\ \vdots & \vdots & \ddots & \vdots \\ x_0^{(N)} & x_1^{(N)} & \cdots & x_d^{(N)} \end{bmatrix} \begin{bmatrix} w_0 \\ w_1\\ \vdots \\ w_d \end{bmatrix} = \mathbf{X}\mathbf{w}$$

 

이제 이것들을 코드로 표현해보자

 

---

from autograd import numpy
from autograd import grad

X = mpg_data[x_cols].values
X = numpy.hstack((numpy.ones((X.shape[0], 1)), X))  #연산의 편의를 위해 x_0=1
y = mpg_data[y_col].values

print("X.shape = {}, y.shape = {}".format(X.shape, y.shape))
# X.shape = (392, 7), y.shape = (392,)

 

이제 cost function을 정의해 보자

mean squared error(MSE)를 사용해보겠다.

 

$$L(\mathbf{y}, \hat{\mathbf{y}}) = \frac{1}{N} \sum_{i=1}^N (y^{(i)} - \hat{y}^{(i)})^2$$

위의 식을 코드로 나타내면 다음과 같다.

numpy.mean( numpy.sum((y-y_pred)**2) )

그러면 linear regression과 cost function을 코드로 나타내 보겟다.

def linear_regression(params, X):
    return numpy.dot(X, params)

def cost_function(params, model, X, y):
    y_pred = model(params, X)
    return numpy.mean( numpy.sum((y-y_pred)**2) )

---

이제 gradient descent로 cost function을 최소로 하는 계수를 찾아보겠다.

autograd.grad함수를 사용해서 기울기를 구하고 랜덤한 값을 넣어서 기울기 값이 잘 구해지는지 알아보자

 

gradient = grad(cost_function)

gradient(numpy.random.rand(X.shape[1]), linear_regression, X, y)

#	array([1.89156671e+06, 1.12205854e+07, 4.23341815e+08, 2.16679469e+08,
#	       6.10689089e+09, 2.87049061e+07, 1.43051715e+08])

기울기 값이 너무 크게 나오는 걸 볼 수 있다.

이것은 입력변수들 중에 특정 값들이 너무 커서 발생한 것이다.

 

이를 보정하기 위해 모든 값들을 0~1사이로 scaling해주겠다.

 

사용할 방법은 min-max scaling으로 scikit-learn패키지에서 가져오겠다.

 

이 normalization의 변환식은 다음과 같다.

 

$$x' = \frac{x - \min(x)}{\max(x)-\min(x)}$$

 

이제 앞의 포스트와 같이 gradient descent를 진행해보겠다.

 

max_iter = 1000
alpha = 0.001
params = numpy.zeros(X.shape[1])

for i in range(max_iter):
    descent = gradient(params, linear_regression, X_scaled, y)
    params = params - descent * alpha
    loss = cost_function(params, linear_regression, X_scaled, y)
    if i%100 == 0:
        print("iteration {}, loss = {}".format(i, loss))

학습이 완료된 후의 param은 다음과 같다.

 

params
#	array([ 27.56342896,  -1.24000064,   1.70706502,  -0.2853682 ,
#	       -23.08358717,   1.06664226,   8.99271696])

아까 위의 식 처럼 우리는 예측 값을 params와 X의 곱으로 나타낼 수 있다.

y_pred_gd = X_scaled @ params

 

---

이제 이 모델이 얼마나 정확한지 알아보겠다.

regression 문제에서는 주로 두 개의 지표가 있다.

 

MAE(mean absolute error), RMSE(root mean squared error)

두 지표의 식은 다음과 같다.

 

 

$$\text{MAE}(\mathbf{y}, \hat{\mathbf{y}}) = \frac{1}{N} \sum_{i=1}^N \left\vert y^{(i)} - \hat{y}^{(i)}\right\vert$$ $$\text{RMSE}(\mathbf{y}, \hat{\mathbf{y}}) =\sqrt{ \frac{1}{N} \sum_{i=1}^N (y^{(i)} - \hat{y}^{(i)})^2}$$

 

이 지표들도 scikit-learn 패키지를 통해 사용해보겠다.

 

from sklearn.metrics import mean_absolute_error, mean_squared_error

mae = mean_absolute_error(y, y_pred_gd)
rmse = mean_squared_error(y, y_pred_gd, squared=False)
print("mae  = {}".format(mae))
print("rmse = {}".format(rmse))

#	mae  = 2.613991601156043
#	rmse = 3.40552056741184

 

예측이 어느 정도 잘 된 것을 확인할 수 있다.

 

---

지금 하는 것은 이해를 위해 기능들을 하나하나 직접 구현하는 것인데

이 부분을 이해한 후에 나중에는 코드 3줄 이하로 이 모든 것을 할 수 있게 된다.