[모두의 딥러닝](keras) #2-2. 선형 회귀(Linear regression) 코딩으로 구현하기

'모두의 딥러닝' 개정 2판의 예제 코드를 이용하여 공부한 글입니다.

 

앞서서 선형 회귀의 이론에 대해 자세히 공부해보았습니다. 해당 글은 바로 아래에 링크해두었습니다.

2021.12.31 - [개인 공부/AI] - [모두의 딥러닝] #2. 선형 회귀(Linear regression) - 평균 제곱 오차(MSE), 경사 하강법(Gradient descent algorithm)

[모두의 딥러닝] #2. 선형 회귀(Linear regression) 이론 - 평균 제곱 오차(MSE), 경사 하강법(Gradient descent

 

이론을 공부해 보았으니, 실습을 통하여 확실히 이해해보는 시간을 가집시다.

 

아래의 세 가지 알고리즘을 코드로 구현해보겠습니다.

1. 평균 제곱 오차(MSE)

2. 최소제곱법
3. 경사 하강법 - 단순 선형 회귀

4. 경사 하강법 - 다중 선형 회귀

 

# 1. 평균 제곱 오차(MSE)

평균 제곱 오차를 python 코딩으로 구현해볼까요? 

import numpy as np

# 기울기 W와 y 절편 b -> 예측값(cost function)
fake_W_b = [3, 76]

# x, y의 데이터의 값 -> 실제값
data = [[2,81],[4,93],[6,91],[8,97]]
x=[i[0] for i in data]
y=[i[1] for i in data]

# y=Wx+b에 W와 b 값을 대입하여 결과를 출력하는 함수
def predict(x):
    return fake_W_b[0]*x+fake_W_b[1]

# MSE 함수
def mse(y, y_hat):
    return ((y-y_hat)**2).mean() # y: 실제값, y_hat: 예측값

# MSE 함수를 각 y 값에 대입하여 최종 값을 구하는 함수
def mse_val(y, predict_result):
    return mse(np.array(y), np.array(predict_result))

# 예측 값이 들어갈 빈 리스트
predict_result = []

# 모든 x 값을 한 번씩 대입하여
for i in range(len(x)):
    # 그 결과에 해당하는 predict_result 리스트를 완성
    predict_result.append(predict(x[i]))
    print("공부시간=%.f, 실제 점수=%.f, 예측 점수=%.f" %(x[i],y[i], predict(x[i])))
    
# 최종 MSE 출력
print("mse 최종값: "+str(mse_val(predict_result,y)))

fake_W_b 은 예측한 cost function의 기울기 W와 y 절편 b를 저장하는 배열입니다.

data에는 실제 x, y 값을 저장해둔 후에, x와 y 배열에 각각 분리하여 저장합니다.

앞에서 본 cost function을 식으로 구현한 것이 mse 함수입니다. 예측값과 실제값을 뺀 후에 제곱한 것의 평균을 코드로 나타냈죠. 이때 mse 함수의 인자로는 실제값과 예측값이 들어가는데 예측값을 계산하는 함수로 predict 함수를 사용합니다. predict 함수로 예측한 y값들을 predict_result 배열에 저장합니다. 

예측값과 실제값이 모두 준비되었으므로 mse_val 함수를 이용하여 mse를 계산하여 출력합니다. 참고로 이때 mse_val()의 두 인자의 순서가 바뀌어져 있는 것을 관찰할 수 있습니다. 그 이유는 mse()에서 두 인자의 차를 제곱하기 때문입니다.

 

# 2. 최소제곱오차(MSE) - 단순 선형 회귀

# -*- coding: utf-8 -*-
import numpy as np

# x 값과 y값
x=[2, 4, 6, 8]
y=[81, 93, 91, 97]

# x와 y의 평균값
mx = np.mean(x)
my = np.mean(y)
print("x의 평균값:", mx)
print("y의 평균값:", my)

# 기울기 공식의 분모
divisor = sum([(mx - i)**2 for i in x])

# 기울기 공식의 분자
def top(x, mx, y, my):
    d = 0
    for i in range(len(x)):
        d += (x[i] - mx) * (y[i] - my)
    return d
dividend = top(x, mx, y, my)

print("분모:", divisor)
print("분자:", dividend)

# 기울기와 y 절편 구하기
a = dividend / divisor
b = my - (mx*a)

# 출력으로 확인
print("기울기 a =", a)
print("y 절편 b =", b)

 

# 3. 경사 하강법 - 단순 선형 회귀

이론에서는 생각을 쉽게 하기 위하여 예측함수 H(x)를 y절편이 없는 일차함수로 단순화하여 생각했었습니다.

코딩으로 구현하면 컴퓨터가 대신 계산해줄 것이기 때문에 실제 이론대로 H(x)를 y절편이 존재하는 H(x)=ax+b로 두겠습니다. cost function을 기울기 a와 y절편 b에 대해 각각 편미분한 값이 0이 될 때까지 반복되도록 코드를 작성합니다.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 공부 시간 X와 성적 Y의 리스트 만들기
data = [[2, 81],[4, 93], [6, 91], [8, 97]]
x = [i[0] for i in data] # 리스트 안 for 문 
y = [i[1] for i in data] # [ 표현식 for 항목 in 리스트 or 튜플 if 조건문 ]

# 그래프로 나타내기
plt.figure(figsize=(8,5))
plt.scatter(x,y)
plt.show()

# 리스트로 되어 있는 x와 y값을 넘파이 배열로 바꾸기(인덱스를 주어 하나씩 불러와 계산이 가능하게 하기 위함)
x_data = np.array(x)
y_data = np.array(y)

# 기울기 a와 절편 b의 값 초기화
a = 0
b = 0

# 학습률 정하기
lr = 0.03

# 몇 번 반복될 지 설정
epochs = 2001

# 경사 하강법 시작
for i in range(epochs): # 에포크 수만큼 반복
    y_pred = a * x_data + b # 예측함수: y를 구하는 식 세우기
    error = y_data - y_pred # cost function: 오차를 구하는 식
    # 오차 함수를 a로 미분한 값
    a_diff = -(2/len(x_data))*sum(x_data * (error))
    # 오차 함수를 b로 미분한 값
    b_diff = -(2/len(x_data))*sum(error)
    
    # 각 parameter를 편미분 값의 반대 방향으로 학습률만큼 이동시킨다.
    a = a - lr*a_diff #학습률을 곱해 기존의 a값 업데이트
    b = b - lr*b_diff #학습률을 곱해 기존의 b값 업데이트
    
    if i % 100 == 0: # 100번 반복될 때마다 현재의 a값, b값 출력
        print("epoch=%.f, 기울기=%.04f, 절편=%.04f, a편미분=%.04f, b편미분=%.04f" %(i, a , b,a_diff,b_diff))

# 앞서 구한 기울기와 절편을 이용해 그래프 다시 그리기
y_pred = a * x_data + b
plt.scatter(x,y)
plt.plot([min(x_data), max(x_data)], [min(y_pred),max(y_pred)]) # 예측함수 그리기
plt.show()

 

위의 코드를 돌려보면 아래와 같은 output을 얻을 수 있습니다.

a는 2.3, b는 79로 수렴합니다.

epoch가 커질 수록, 즉 반복될 수록, a편미분 값이 0으로 수렴하는 것을 볼 수 있습니다. 마찬가지로 b편미분 값도 0으로 수렴하는 것을 볼 수가 있죠. 편미분 값이 0이 되기 때문에 자연스럽게 기울기 a와 y절편 b도 수렴하였습니다.

 

학습률을 키워볼까요? 0.3에서 0.4로 키워봅시다. 이동의 폭을 키우는 것이죠.

해석하기에 매우 곤란한 값이 나옵니다. 그래프조차 그릴 수 없는 숫자가 아닌 값입니다. 이론에서 공부한 바로 추측하면, 학습률이 너무 커서 발산한 것으로 예상할 수 있습니다. 이런 경우에는 다른 학습률로 다시 시도해봐야 합니다.

 

그럼 학습률을 줄어봅시다. 0.3에서 0.2로 줄여봅시다. 이동의 폭을 줄이는 것입니다.

다행히 수로 결과가 나오긴 합니다. b편미분 값은 0이 아닌 것으로 보아 b는 수렴하지 않은 것을 알 수 있습니다만 a편미분 값은 0.0000으로 수렴한 것 같기도 합니다.

하지만 꺼림칙 하죠? a편미분 값이 0.0000이라고 한데 정말 수렴한 것일까요? 서식지정자를 0.4f에서 0.8f로 수정하여 확인해보면 되죠.

그렇습니다. 역시나 수렴하지 않았네요. 수렴하기에는 학습률이 너무 작았습니다.

 

위의 시행을 통해서 우리는 최적의 학습률이 0.03이라는 것을 알 수 있습니다.

 

# 4. 경사 하강법 - 다중 선형 회귀

다중 선형 회귀라고 단순 선형 회귀의 논리와 다르지 않습니다.

x 데이터 값이 한 개가 아니라는 것만 빼면 모두 같습니다. 다만, 3D의 그래프가 필요하기에 코드는 다음과 같이 다르게 작성해주어야 합니다.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

#공부시간 X와 성적 Y의 리스트를 만듭니다.
data = [[2, 0, 81], [4, 4, 93], [6, 2, 91], [8, 3, 97]]
x1 = [i[0] for i in data]
x2 = [i[1] for i in data]
y = [i[2] for i in data]

#그래프로 확인해 봅니다.
ax = plt.axes(projection='3d')
ax.set_xlabel('study_hours')
ax.set_ylabel('private_class')
ax.set_zlabel('Score')
ax.dist = 11 
ax.scatter(x1, x2, y)
plt.show()

 

 

#리스트로 되어 있는 x와 y값을 넘파이 배열로 바꾸어 줍니다.(인덱스를 주어 하나씩 불러와 계산이 가능해 지도록 하기 위함입니다.)
x1_data = np.array(x1)
x2_data = np.array(x2)
y_data = np.array(y)

# 기울기 a와 절편 b의 값을 초기화 합니다.
a1 = 0
a2 = 0
b = 0

#학습률을 정합니다.
lr = 0.02 

#몇 번 반복될지를 설정합니다.(0부터 세므로 원하는 반복 횟수에 +1을 해 주어야 합니다.)
epochs = 2001 

#경사 하강법을 시작합니다.
for i in range(epochs): # epoch 수 만큼 반복
    y_pred = a1 * x1_data + a2 * x2_data + b  #y를 구하는 식을 세웁니다
    error = y_data - y_pred  #오차를 구하는 식입니다.
    a1_diff = -(2/len(x1_data)) * sum(x1_data * (error)) # 오차함수를 a1로 미분한 값입니다. 
    a2_diff = -(2/len(x2_data)) * sum(x2_data * (error)) # 오차함수를 a2로 미분한 값입니다. 
    b_new = -(2/len(x1_data)) * sum(y_data - y_pred)  # 오차함수를 b로 미분한 값입니다. 
    a1 = a1 - lr * a1_diff  # 학습률을 곱해 기존의 a1값을 업데이트합니다.
    a2 = a2 - lr * a2_diff  # 학습률을 곱해 기존의 a2값을 업데이트합니다.
    b = b - lr * b_new  # 학습률을 곱해 기존의 b값을 업데이트합니다.
    if i % 100 == 0:    # 100번 반복될 때마다 현재의 a1, a2, b값을 출력합니다.
        print("epoch=%.f, 기울기1=%.04f, 기울기2=%.04f, 절편=%.04f, a1편미분=%.04f, a2편미분=%.04f, b편미분=%.04f" % (i, a1, a2, b,a1_diff,a2_diff,b_new))

단순 선형 회귀와 마찬가지로 세 parameter의 편미분 값이 0으로 향하는 것을 볼 수 있습니다. 이를 통해 학습률의 최적화가 잘 됐다는 것을 알 수 있죠.

#참고 자료, 다중 선형회귀 '예측 평면' 3D로 보기

import statsmodels.api as statm
import statsmodels.formula.api as statfa
#from matplotlib.pyplot import figure

X = [i[0:2] for i in data]
y = [i[2] for i in data]

X_1=statm.add_constant(X)
results=statm.OLS(y,X_1).fit()

hour_class=pd.DataFrame(X,columns=['study_hours','private_class'])
hour_class['Score']=pd.Series(y)

model = statfa.ols(formula='Score ~ study_hours + private_class', data=hour_class)

results_formula = model.fit()

a, b = np.meshgrid(np.linspace(hour_class.study_hours.min(),hour_class.study_hours.max(),100),
                   np.linspace(hour_class.private_class.min(),hour_class.private_class.max(),100))

X_ax = pd.DataFrame({'study_hours': a.ravel(), 'private_class': b.ravel()})
fittedY=results_formula.predict(exog=X_ax)

fig = plt.figure()
graph = fig.add_subplot(111, projection='3d')

graph.scatter(hour_class['study_hours'],hour_class['private_class'],hour_class['Score'],
              c='blue',marker='o', alpha=1)
graph.plot_surface(a,b,fittedY.values.reshape(a.shape),
                   rstride=1, cstride=1, color='none', alpha=0.4)
graph.set_xlabel('study hours')
graph.set_ylabel('private class')
graph.set_zlabel('Score')
graph.dist = 11

plt.show()

1차원 예측 직선이 3차원 예측 평면으로 표현된 것을 볼 수 있습니다. 다중 선형 회귀도 위의 단순 선형 회귀의 경우와 같이 학습률을 변경해보며 최적화할 수 있습니다.