%load_ext watermark
%watermark -v -p numpy,scipy,sklearn,pandas,matplotlib
# 파이썬 2와 파이썬 3 지원
from __future__ import division, print_function, unicode_literals
# 공통
import numpy as np
import os
# 일관된 출력을 위해 유사난수 초기화
np.random.seed(42)
# 맷플롯립 설정
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
# 한글출력
matplotlib.rc('font', family='AppleGothic')
plt.rcParams['axes.unicode_minus'] = False
The watermark extension is already loaded. To reload it, use: %reload_ext watermark CPython 3.5.4 IPython 6.2.1 numpy 1.14.0 scipy 1.0.0 sklearn 0.20.0 pandas 0.22.0 matplotlib 2.1.2
다음은 1장에서 본 삶의 만족도에 대한 간단한 선형 회귀 모델
삶의 만족도 = $\theta$0 + $\theta$1 * 1인당_GDP
일반적으로 선형 모델은 [식 4-1]에서처럼 입력 특성의 가중치 합과 편향bias(또는 절편intercept)이라는 상수를 더해 예측을 만듦
이 식은 [식 4-2]처럼 벡터 형태로 더 간단하게 쓸 수 있음
이를 위해 먼저 모델이 훈련 데이터에 잘 들어맞는지 측정해야 함.
성능 측정 지표로는 평균 제곱 오차Mean Square Error(MSE)를 최소화
훈련 세트 X에 대한 선형 회귀 가설 h$\theta$의 MSE는 [식 4-3]처럼 계산함
정규방정식Normal Equation : 비용 함수를 최소화하는 $\theta$ 값을 얻을 수 있는 수학 공식(식 4-4)
이 공식을 테스트하기 위해 선형처럼 보이는 데이터를 생성
import numpy as np
X = 2 * np.random.rand(100, 1) #균일 분포[0,1) , 100행 1열
y = 4 + 3 * X + np.random.randn(100, 1) #정규분포
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([0, 2, 0, 15])
plt.show()
정규방정식을 사용해 $\hat{\theta}$을 계산
X_b = np.c_[np.ones((100, 1)), X] # 모든 샘플에 x0 = 1을 추가합니다. 두 개의 1차원 배열을 칼럼으로 세로로 붙여서 2차원 배열 만듦
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y) #inv():역행렬 계산, dot():행렬 곱셈
theta_best
array([[4.21509616], [2.77011339]])
$\hat{\theta}$을 사용한 예측
X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new] # 모든 샘플에 x0 = 1을 추가합니다.
y_predict = X_new_b.dot(theta_best)
y_predict
array([[4.21509616], [9.75532293]])
아래 그래프는 모델의 예측을 나타냄
plt.plot(X_new, y_predict, "r-", linewidth=2, label="예측") # x범위, y범위, 스타일=빨간색 선
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.legend(loc="upper left", fontsize=14)
plt.axis([0, 2, 0, 15])
plt.show()
아래는 같은 작업을 하는 사이킷런 코드
from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X, y)
lin_reg.intercept_, lin_reg.coef_
(array([4.21509616]), array([[2.77011339]]))
lin_reg.predict(X_new)
array([[4.21509616], [9.75532293]])
특성 수가 두 배로 늘어나면 계산 시간이 대략 5.3 ~ 8배 증가
하지만 이 공식의 복잡도가 훈련 세트의 샘플 수에는 선형적으로 증가(즉, O(m))
$\theta$를 임의의 값으로 시작해서 한번에 조금씩 비용 함수가 감소되는 방향으로 진행하여 알고리즘이 최솟값에 수렴할 때까지 점진적으로 향상 (그림 4-3)
[그림 4-6]은 경사 하강법의 두 가지 문제점을 보여줌
[그림 4-7]에서 왼쪽의 경사 하강법 알고리즘이 최솟값으로 곧장 진행하고 있어 빠르게 도달하지만 오른쪽 그래프는 돌아서 나가기 때문에 시간이 오래 걸림
편도함수를 각각 계산하는 대신 [식 4-6]을 사용하여 한꺼번에 계산
그래디언트 벡터 $\nabla$$\theta$MSE($\theta$)는 비용 함수의 (모델 파라미터마다 한 개씩인) 편도함수를 모두 담고 있음
아래는 알고리즘을 구현한 코드
theta_path_bgd = [] # 알고리즘이 훈련 과정 동안 파라미터 공간에서 움직인 경로
np.random.seed(42) # 초기 난수 생성값 지정
eta = 0.1 #학습률
n_iterations = 1000
m = 100
theta = np.random.randn(2,1)
for iteration in range(n_iterations):
gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
theta = theta - eta * gradients
theta_path_bgd.append(theta)
theta
array([[4.21509616], [2.77011339]])
[그림 4-8]은 세 가지 다른 학습률을 사용하여 진행한 경사 하강법의 스텝 처음 10개를 보여줌(점선은 시작점)
def plot_gradient_descent(theta, eta, theta_path=None):
plt.plot(X, y, "b.")
for iteration in range(n_iterations):
if iteration < 10:
y_predict = X_new_b.dot(theta)
style = "b-" if iteration > 0 else "r--"
plt.plot(X_new, y_predict, style)
gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
theta = theta - eta * gradients
plt.xlabel("$x_1$", fontsize=18)
plt.axis([0, 2, 0, 15])
plt.title(r"$\eta = {}$".format(eta), fontsize=16)
np.random.seed(42) # 초기 난수 생성값 지정
theta = np.random.randn(2,1) # random initialization
plt.figure(figsize=(10,4))
plt.subplot(131); plot_gradient_descent(theta, eta=0.02)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.subplot(132); plot_gradient_descent(theta, eta=0.1, theta_path=theta_path_bgd)
plt.subplot(133); plot_gradient_descent(theta, eta=0.5)
plt.show()
반복 횟수를 아주 크게 지정하고 그래디언트 벡터가 아주 작아지면, 즉 벡터의 노름이 어떤 값 $\varepsilon$(허용오차tolerance)보다 작아지면 경사 하강법이 (거의) 최솟값에 도달한 것이므로 알고리즘을 중지함
비용 함수가 최솟값에 다다를 때까지 부드럽게 감소하지 않고 위아래로 요동치면서 평균적으로 감소함. 시간이 지나면 최솟값에 매우 근접하겠지만 최적치는 아님 (그림 4-9 참조)
시작할 때는 학습률을 크게 하고(수렴 빠름, 지역 최솟값에 빠지지 않음) 점차 작게 줄여서 알고리즘이 전역 최솟값에 도달하게 함
아래은 확률적 경사 하강법의 구현 코드
theta_path_sgd = [] # 알고리즘이 훈련 과정 동안 파라미터 공간에서 움직인 경로
m = len(X_b) # 훈련 세트에 있는 샘플 수
np.random.seed(42) # 초기 난수 생성값 지정
n_epochs = 50
t0, t1 = 5, 50 # 학습 스케줄 하이퍼파라미터 learning schedule hyperparameters
def learning_schedule(t):
return t0 / (t + t1)
theta = np.random.randn(2,1) # 무작위 초기화
for epoch in range(n_epochs):
for i in range(m):
if epoch == 0 and i < 20: # 훈련 스텝의 첫 20개를 그림
y_predict = X_new_b.dot(theta)
style = "b-" if i > 0 else "r--"
plt.plot(X_new, y_predict, style)
random_index = np.random.randint(m)
xi = X_b[random_index:random_index+1]
yi = y[random_index:random_index+1]
gradients = 2 * xi.T.dot(xi.dot(theta) - yi) # 하나의 샘플에 대한 그래디언트를 계산하므로 식 4-6에서 /m 이 없음
eta = learning_schedule(epoch * m + i)
theta = theta - eta * gradients # 식 4-7
theta_path_sgd.append(theta)
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([0, 2, 0, 15])
plt.show()
theta
array([[4.21076011], [2.74856079]])
사이킷런에서 SGD 방식으로 선형 회귀를 사용하려면 SGDRegressor 클래스를 사용
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter=50, penalty=None, eta0=0.1, random_state=42) # 학습률 0.1, 에포크 50번
sgd_reg.fit(X, y.ravel()) # 다차원 배열(array)을 1차원 배열로 평평하게 펴줌, NumPy의 reshape()함수와 반대 기능
SGDRegressor(alpha=0.0001, average=False, early_stopping=False, epsilon=0.1, eta0=0.1, fit_intercept=True, l1_ratio=0.15, learning_rate='invscaling', loss='squared_loss', max_iter=50, n_iter=None, n_iter_no_change=5, penalty=None, power_t=0.25, random_state=42, shuffle=True, tol=None, validation_fraction=0.1, verbose=0, warm_start=False)
sgd_reg.intercept_, sgd_reg.coef_
(array([4.16782089]), array([2.72603052]))
theta_path_mgd = [] # 알고리즘이 훈련 과정 동안 파라미터 공간에서 움직인 경로
n_iterations = 50
minibatch_size = 20
np.random.seed(42) # 초기 난수 생성값 지정
theta = np.random.randn(2,1) # 무작위 초기화
t0, t1 = 200, 1000
def learning_schedule(t):
return t0 / (t + t1)
t = 0
for epoch in range(n_iterations):
shuffled_indices = np.random.permutation(m)
X_b_shuffled = X_b[shuffled_indices]
y_shuffled = y[shuffled_indices]
for i in range(0, m, minibatch_size):
t += 1
xi = X_b_shuffled[i:i+minibatch_size]
yi = y_shuffled[i:i+minibatch_size]
gradients = 2/minibatch_size * xi.T.dot(xi.dot(theta) - yi) # 식 4-6
eta = learning_schedule(t)
theta = theta - eta * gradients # 식 4-7
theta_path_mgd.append(theta)
theta
array([[4.25214635], [2.7896408 ]])
아래 그림은 세 가지 경사 하강법 알고리즘이 훈련 과정 동안 파라미터 공간에서 움직인 경로를 나타냄
theta_path_bgd = np.array(theta_path_bgd)
theta_path_sgd = np.array(theta_path_sgd)
theta_path_mgd = np.array(theta_path_mgd)
plt.figure(figsize=(7,4))
plt.plot(theta_path_sgd[:, 0], theta_path_sgd[:, 1], "r-s", linewidth=1, label="SGD")
plt.plot(theta_path_mgd[:, 0], theta_path_mgd[:, 1], "g-+", linewidth=2, label="미니배치")
plt.plot(theta_path_bgd[:, 0], theta_path_bgd[:, 1], "b-o", linewidth=3, label="배치")
plt.legend(loc="upper left", fontsize=16)
plt.xlabel(r"$\theta_0$", fontsize=20)
plt.ylabel(r"$\theta_1$ ", fontsize=20, rotation=0)
plt.axis([2.5, 4.5, 2.3, 3.9])
plt.show()
지금까지 논의한 알고리즘을 선형 회귀를 사용해 비교함(m은 훈련 샘플 수, n은 특성 수). [표 4-1] 참조
먼저 간단한 2차방정식quadratic equation으로 비선형 데이터를 생성(약간의 노이즈 포함)
np.random.seed(42)
m = 100
X = 6 * np.random.rand(m, 1) - 3
y = 0.5 * X**2 + X + 2 + np.random.randn(m, 1)
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([-3, 3, 0, 10])
plt.show()
사이킷런의 PolynomialFeatures를 사용해 훈련 데이터를 변환
from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly_features.fit_transform(X)
X[0]
array([-0.75275929])
X_poly[0]
array([-0.75275929, 0.56664654])
이 확장된 훈련 데이터에 Linear Regression을 적용
lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)
lin_reg.intercept_, lin_reg.coef_
(array([1.78134581]), array([[0.93366893, 0.56456263]]))
X_new=np.linspace(-3, 3, 100).reshape(100, 1) # 시작, 끝점을 균일간격으로 100개 생성
X_new_poly = poly_features.transform(X_new)
y_new = lin_reg.predict(X_new_poly)
plt.plot(X, y, "b.")
plt.plot(X_new, y_new, "r-", linewidth=2, label="예측")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.legend(loc="upper left", fontsize=14)
plt.axis([-3, 3, 0, 10])
plt.show()