4장 모델 훈련

Outline

. 선형 회귀

. 경사 하강법

. 다항 회귀

. 학습 곡선

. 규제가 있는 선형 모델

. 로지스틱 회귀

4.0 설정

In [1]:
# 공통
import numpy as np
import os
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib import font_manager, rc
from scipy.linalg import lstsq
from scipy.linalg import pinv
from sklearn.linear_model import SGDRegressor
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn import datasets

from matplotlib.colors import ListedColormap
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso
from sklearn.linear_model import ElasticNet
from sklearn.linear_model import LogisticRegression

# 일관된 출력을 위해 유사난수 초기화
np.random.seed(42)

# 맷플롯립 설정
font_name = font_manager.FontProperties(fname="c:/Windows/Fonts/malgun.ttf").get_name()
rc('font', family=font_name)

plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
plt.rcParams['axes.unicode_minus'] = False

# 그림을 저장할 폴드
PROJECT_ROOT_DIR = "C:/Users/Admin/Desktop/ML/"
# PROJECT_ROOT_DIR = "C:/Users/sally/Desktop/ML/"
# PROJECT_ROOT_DIR = "C:/Users/User/Desktop/ML/"
# PROJECT_ROOT_DIR = "C:/Users/sally/Dropbox/2019-Fall-Semester/ML"

CHAPTER_ID = "training_linear_models"

IMAGES_PATH = os.path.join(PROJECT_ROOT_DIR, "images", CHAPTER_ID)

def save_fig(fig_id, tight_layout=True):
    path = os.path.join(IMAGES_PATH, fig_id + ".png")
    if tight_layout:
        plt.tight_layout()
    plt.savefig(path, format='png', dpi=300)

4.1 선형 회귀

. 선형 회귀 모델의 예측

.. $\hat{y}=h_\theta(x)=\theta'x$

... $h$ : 모델 파라미터 $\theta$를 사용한 가설 (hypothesis) 함수

... $\theta=(\theta_0,\theta_1,\ldots,\theta_n)'$

... $x=(1,x_1,x_2,\ldots,x_n)'$

. 선형 회귀 모델의 MSE 비용 함수

.. ${\rm MSE}(\theta)=\frac{1}{m}\sum_{i=1}^m(\theta'x^{(i)}-y^{(i)})^2$

4.1.1 정규 방정식

. 해석적인 방법 : $\hat\theta = (X'X)^{-1}X'y$ Why?

.. $X=(x^{(1)'},x^{(2)'},\ldots,x^{(m)'})': m\times (n+1)$ 특성 행렬

.. $y=(y^{(1)},y^{(2)},\ldots,y^{(m)})': m\times1$ 반응 벡터

In [2]:
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)
In [3]:
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([0, 2, 0, 15])
save_fig("generated_data_plot")
plt.show()
In [4]:
# 모든 샘플에 x0 = 1을 추가
X_b = np.c_[np.ones((100, 1)), X]
In [5]:
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
In [6]:
theta_best
Out[6]:
array([[4.21509616],
       [2.77011339]])
In [7]:
X_new = np.array([[0], [2]])
# 모든 샘플에 x0 = 1을 추가
X_new_b = np.c_[np.ones((2, 1)), X_new]  
In [8]:
y_predict = X_new_b.dot(theta_best)
y_predict
Out[8]:
array([[4.21509616],
       [9.75532293]])
In [9]:
plt.plot(X_new, y_predict, "r-", linewidth=2, label="Prediction")
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])
save_fig("linear_model_predictions")
plt.show()

. scikitlearn을 이용한 예측의 결과와 비교

. sklearn.linear_model.LinearRegression : 최소제곱법에 의한 선형 회귀

http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html

In [10]:
lin_reg = LinearRegression()
lin_reg.fit(X, y)
Out[10]:
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None,
         normalize=False)
In [11]:
# 회귀 계수 추정량
lin_reg.intercept_, lin_reg.coef_
Out[11]:
(array([4.21509616]), array([[2.77011339]]))
In [12]:
lin_reg.predict(X_new)
Out[12]:
array([[4.21509616],
       [9.75532293]])

. scipy.linalg : 선형대수 함수 모임

https://docs.scipy.org/doc/scipy/reference/linalg.html

. scipy.linalg.lstsq를 사용해서 직접 계산 : 방정식 $Ax = b$의 해를 직접 계산

https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html#scipy.linalg.lstsq

In [13]:
theta_best_svd, residuals, rank, s = lstsq(X_b, y)
# theta_best_svd, residuals, rank, s = np.linalg.lstsq(X_b, y, rcond=1e-6)
theta_best_svd
Out[13]:
array([[4.21509616],
       [2.77011339]])

. numpy.linalg.pinv를 사용해서 유사역행렬을 직접 계산

In [14]:
# pinv(X_b).dot(y)
np.linalg.pinv(X_b).dot(y)
Out[14]:
array([[4.21509616],
       [2.77011339]])

4.1.2 계산 복잡도 (computational complexity)

. $O(n^{2.4})\sim O(n^3)$

.. 특성 수가 두 배 늘어나면 $2^{2.4}\sim 2^3$ 배 증가

. 그러나 샘플 수는 두 배 늘어나도 선형적으로 증가 Why?

. 결국 특성이 매우 많고 훈련 샘플이 너무 많을 때는 불리한 방법

4.2 경사 하강법

. 원리 : $\theta$를 임의의 값(random initialization)으로 시작해서 한번에 조금씩 비용 함수가 감소되는 방향으로 진행하여 알고리즘이 최솟값에 도달할 때까지 점진적으로 향상시킴

. 비용 함수

.. 볼록(convex) 함수

.. 기울기가 갑자기 변하지 않음 (즉, Lipschitz 연속 https://en.wikipedia.org/wiki/Lipschitz_continuity)

... 전역 최솟값(global minimum)에 가깝게 접근 가능

.. 특성들의 스케일에 따라 비용함수 모양이 다양

Cap%202018-09-13%2011-04-07-679.png

Cap%202018-09-13%2011-29-06-709.png

Cap%202018-09-13%2011-41-02-234.jpg

. 스텝의 크기 : 학습률 하이퍼파라미터로 결정

Cap%202018-09-13%2011-17-17-858.png

Cap%202018-09-13%2011-17-33-696.png

4.2.1 뱃지 경사 하강법 (batch gradient descent)

. ${\rm MSE}(\theta)=\frac{1}{m}(y-X\theta)'(y-X\theta)=\frac{1}{m}(y'y-2\theta'(X'y)+\theta'(X'X)\theta)$

. $\frac{\partial {\rm MSE}(\theta)}{\partial\theta}=\nabla_\theta {\rm MSE}(\theta)=\frac{2}{m}X'(X\theta-y)$

. 경사 하강법의 스텝 : $\theta^{\rm new}=\theta^{\rm old}-\eta \nabla_\theta {\rm MSE}(\theta^{\rm old})$

.. $\eta$ : 학습률

. 특징 :

.. 매 스텝에서 훈련 세트 전체를 사용

.. 큰 훈련 세트에서는 느림

.. 그러나 특성 수에 민감하지 않음

. 경사하강법을 직접 구현해보면 ...

In [15]:
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
In [16]:
theta
Out[16]:
array([[4.21509616],
       [2.77011339]])

. 학습률 $\eta$가 회귀 추정에 미치는 영향을 그려보면 ...

In [17]:
def plot_gradient_descent(theta, eta, theta_path=None):
    m = len(X_b)
    plt.plot(X, y, "b.")
    n_iterations = 1000
    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
        if theta_path is not None:
            theta_path.append(theta)
    plt.xlabel("$x_1$", fontsize=18)
    plt.axis([0, 2, 0, 15])
    plt.title(r"$\eta = {}$".format(eta), fontsize=16)
In [18]:
theta_path_bgd = []
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)
save_fig("gradient_descent_plot")
plt.show()
In [19]:
theta #초기값
Out[19]:
array([[ 0.49671415],
       [-0.1382643 ]])

4.2.2 확률적 경사 하강법

. 특징 :

.. 매 스텝에서 한 개의 샘플을 무작위로 선택하고 그 샘플에 대한 경사를 계산

.. 알고리즘 수행 속도가 빠름

.. 그러나 배치 경사 하강법보다 알고리즘이 불안정

.. 최솟값에 안착하지 못할 수도!

. 해결책 :

.. 학습률 조정 : 시작할 때는 크게, 점차 작게

.. 매 스텝에서 학습 스케줄링

Cap%202018-09-18%2013-08-36-607.png

. 확률적경사하강법을 직접 구현한다면 ...

In [20]:
theta_path_sgd = []
m = len(X_b)
n_epochs = 50
t0, t1 = 5, 50  # 학습 스케줄 하이퍼파라미터 

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:                    
            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)
        eta = learning_schedule(epoch * m + i)
        theta = theta - eta * gradients
        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])                              
save_fig("sgd_plot")                                 
plt.show()                                           
In [21]:
theta
Out[21]:
array([[4.20742938],
       [2.74264448]])
In [22]:
# 기본 학습 스케줄링 
sgd_reg = SGDRegressor(max_iter=50, penalty=None, eta0=0.1, random_state=42)
# sgd_reg.fit(X, y.ravel())
sgd_reg.fit(X, y)
C:\ProgramData\Anaconda3\lib\site-packages\sklearn\linear_model\stochastic_gradient.py:183: FutureWarning: max_iter and tol parameters have been added in SGDRegressor in 0.19. If max_iter is set but tol is left unset, the default value for tol in 0.19 and 0.20 will be None (which is equivalent to -infinity, so it has no effect) but will change in 0.21 to 1e-3. Specify tol to silence this warning.
  FutureWarning)
C:\ProgramData\Anaconda3\lib\site-packages\sklearn\utils\validation.py:761: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().
  y = column_or_1d(y, warn=True)
Out[22]:
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)
In [23]:
sgd_reg.intercept_, sgd_reg.coef_
Out[23]:
(array([4.16782089]), array([2.72603052]))

4.2.3 미니뱃지 경사 하강법

. 임의의 작은 샘플 세트에 대해 경사를 계산

. SGD보다 덜 불규칙적임

. local minimum에서 빠져나오지 못할 수도!

In [24]:
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)
        eta = learning_schedule(t)
        theta = theta - eta * gradients
        theta_path_mgd.append(theta)
In [25]:
theta #해
Out[25]:
array([[4.1984193 ],
       [2.76598192]])
In [26]:
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])
save_fig("gradient_descent_paths_plot")
plt.show()

4.3 다항 회귀

. 특성의 거듭제곱을 새로운 특성으로 추가한 모델

. PolynomialFeatures(degree=d)는 특성이 $n$ 개인 배열을 특성이 $n+d\choose d$ 개인 배열로 변환 Why?

.. eg, 특성아 $a, b$ 두 개이고 degree = 3일 때, 가능한 특성은 $5\choose3$ = $5\choose2$ = $10$ 개 즉, $1,a, b, a^2, b^2, ab,a^3, b^3, a^2b, ab^2$

.. $\begin{eqnarray*}

H(n,0)+H(n,1)+\cdots+H(n,d)&=&H(n+1,0)+H(n,1)+\cdots+H(n,d)\ &=&H(n+1,1)+H(n,2)+\cdots+H(n,d)\ &\cdots&\ &=&H(n+1,d-1)+H(n,d)\ &=&H(n+1,d) = {n+d\choose d} \end{eqnarray*}$

.. 점화식 $H(n+1,d)=H(n+1,d-1)+H(n,d)$ Why?

. 가상 자료 생성

In [27]:
m = 100
X = 6 * np.random.rand(m, 1) - 3
y = 0.5 * X**2 + X + 2 + np.random.randn(m, 1)
In [28]:
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([-3, 3, 0, 10])
save_fig("quadratic_data_plot")
plt.show()

. 2차 곡선 적합

In [29]:
poly_features = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly_features.fit_transform(X)
In [30]:
X_poly[:5]
Out[30]:
array([[-2.46208661,  6.06187049],
       [-0.97584553,  0.9522745 ],
       [-2.46179449,  6.06043209],
       [ 2.76926234,  7.66881389],
       [ 1.8435699 ,  3.39874996]])
In [31]:
lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)
Out[31]:
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None,
         normalize=False)
In [32]:
lin_reg.intercept_, lin_reg.coef_
Out[32]:
(array([1.7594576]), array([[1.01802176, 0.50874904]]))
In [33]:
X_new = np.linspace(-3, 3, 100).reshape(100, 1)
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="Prediction")
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])
save_fig("quadratic_data_plot")
plt.show()

. 1차, 2차, 300차 다항 회귀 비교

In [34]:
for style, width, degree in (("g-", 1, 300), ("b--", 2, 2), ("r-+", 2, 1)):
    polybig_features = PolynomialFeatures(degree=degree, include_bias=False)
    std_scaler = StandardScaler()
    lin_reg = LinearRegression()
    polynomial_regression = Pipeline([
            ("poly_features", polybig_features),
            ("std_scaler", std_scaler),
            ("lin_reg", lin_reg),
        ])
    polynomial_regression.fit(X, y)
    y_newbig = polynomial_regression.predict(X_new)
    plt.plot(X_new, y_newbig, style, label=str(degree), linewidth=width)

plt.plot(X, y, "b.", linewidth=3)
plt.legend(loc="upper left")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([-3, 3, 0, 10])
save_fig("high_degree_polynomials_plot")
plt.show()

4.4 학습 곡선

. 얼마나 복잡한 모델을 사용해야 하나?

. 어떻게 모델이 데이터에 과대적합 혹은 과소적합되었는지 알 수 있는가?

.. 교차 검증 : 훈련 세트에서는 좋은데 교차 검증이 좋지 않으면 과대적합, 둘 다 좋지 않으면 과소적합

.. 학습 곡선 : 훈련 세트와 검증 세트의 모델 성능을 훈련 세트의 크기로 나타낸 곡선

In [35]:
def plot_learning_curves(model, X, y):
    X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=10)
    train_errors, val_errors = [], []
    for m in range(1, len(X_train)):
        model.fit(X_train[:m], y_train[:m])
        y_train_predict = model.predict(X_train[:m])
        y_val_predict = model.predict(X_val)
        train_errors.append(mean_squared_error(y_train[:m], y_train_predict))
        val_errors.append(mean_squared_error(y_val, y_val_predict))

    plt.plot(np.sqrt(train_errors), "r-+", linewidth=2, label="훈련")
    plt.plot(np.sqrt(val_errors), "b-", linewidth=3, label="검증")
    plt.legend(loc="upper right", fontsize=14) 
    plt.xlabel("훈련 세트 크기", fontsize=14)    
    plt.ylabel("RMSE", fontsize=14)             

. 단순 선형 회귀 모델 가정하에서 ...

In [36]:
plot_learning_curves(lin_reg, X, y)

plt.axis([0, 80, 0, 3])
save_fig("underfitting_learning_curves_plot")  
plt.show()                      

. 10차 다항 회귀 모델 가정하에서

In [37]:
polynomial_regression = Pipeline([
        ("poly_features", PolynomialFeatures(degree = 10, include_bias=False)),
        ("lin_reg", LinearRegression()),
    ])

plot_learning_curves(polynomial_regression, X, y)
plt.axis([0, 80, 0, 3])     
save_fig("learning_curves_plot")  
plt.show()                      

편향과 분산 트레이드 오프

. 편향(bias) : 잘못된 가정으로 인한 오차. 편향이 크면 과소적합 위험

. 분산(variance) : 훈련 데이터의 작은 변동에 모델이 과도하게 민감. 자유도가 높은 모델은 높은 분산을 갖기 쉬움. 과대적함 위험

. 줄일 수 없는 오차(irreducible error) : 데이터 자체의 노이즈로 인해 발생

. 모델의 복잡도가 커지면 편향은 줄지만 분산은 늘어나고, 복잡도가 줄어들면 편향은 커지고 분산은 작아짐

4.5 규제가 있는 선형 모델

. 모델의 가중치를 제한하거나 다항식의 차수를 감소

4.5.1 릿지(ridge) 회귀

. 비용 함수 :

.. $J(\theta)={\rm MSE}(\theta)+\alpha \frac{1}{2}\sum_{i=1}^n\theta_i^2$

.. intercept $(\theta_0)$는 규제항에 포함하지 않음

.. 입력 특성의 스케일에 의존함. 입력 특성의 표준화가 요구됨

.. $\alpha$ : 모델을 얼마나 규제할 지를 조절하는 하이퍼파라미터 :

... $\alpha=0$ : 선형 회귀, $\alpha$가 매우 크면 모든 가중치가 $0$에 가까워지고 데이터의 평균을 지나는 직선에 가까움

. 비용 함수의 경사 벡터 : $\frac{2}{m}X'(X\theta-y)+\alpha\theta$

. 정규방정식 해 : $\hat\theta=(X'X+\alpha A)^{-1}X'y$ :

.. $A$ : $(n+1)\times(n+1)$, $(0,0)$의 원소는 $0$인 단위 행렬

. numpy.random.rand : $U(0,1)$에서 난수 생성

https://docs.scipy.org/doc/numpy-1.14.1/reference/generated/numpy.random.rand.html

. 내장 함수

https://wikidocs.net/32#zip

. sklearn.linear_model.Ridge : $L_2$ 규제를 가진 선형 회귀 모델 적합

http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html

. sklearn.linear_model.SGDRegressor : $L_2$ 규제를 가진 회귀 모형 적합, SGD 방법을 사용해서

http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.SGDRegressor.html

. 가상 자료 생성

In [38]:
m = 20
X = 3 * np.random.rand(m, 1)
y = 1 + 0.5 * X + np.random.randn(m, 1) / 1.5
X_new = np.linspace(0, 3, 100).reshape(100, 1)
In [39]:
def plot_model(model_class, polynomial, alphas, **model_kargs):
    for alpha, style in zip(alphas, ("b-", "g--", "r:")):
        model = model_class(alpha, **model_kargs) if alpha > 0 else LinearRegression()
        if polynomial:
            model = Pipeline([
                    ("poly_features", PolynomialFeatures(degree=10, include_bias=False)),
                    ("std_scaler", StandardScaler()),
                    ("regul_reg", model),
                ])
        model.fit(X, y)
        y_new_regul = model.predict(X_new)
        lw = 2 if alpha > 0 else 1
        plt.plot(X_new, y_new_regul, style, linewidth=lw, label=r"$\alpha = {}$".format(alpha))
    plt.plot(X, y, "b.", linewidth=3)
    plt.legend(loc="upper left", fontsize=15)
    plt.xlabel("$x_1$", fontsize=18)
    plt.axis([0, 3, 0, 4])

plt.figure(figsize=(8,4))
plt.subplot(121)
plot_model(Ridge, polynomial=False, alphas=(0, 10, 100), random_state=42)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.subplot(122)
plot_model(Ridge, polynomial=True, alphas=(0, 10**-5, 1), random_state=42)
save_fig("ridge_regression_plot")
plt.show()

. Ridge 클래스와 SGD Regressor 클래스의 비교

In [40]:
# Chokesky decomposition of $(X'X+\alpha A)$ : LL', L : lower traiangular matrix
ridge_reg = Ridge(alpha=1, solver="cholesky", random_state=42)
ridge_reg.fit(X, y)
Out[40]:
Ridge(alpha=1, copy_X=True, fit_intercept=True, max_iter=None,
   normalize=False, random_state=42, solver='cholesky', tol=0.001)
In [41]:
ridge_reg.predict([[1.5]])
Out[41]:
array([[1.64129117]])
In [42]:
# SGDRegressor로 직접 계산

# sgd_reg = SGDRegressor(max_iter=5, penalty="l2", random_state=42)
sgd_reg = SGDRegressor(alpha=1,n_iter=1000, penalty="l2", random_state=42)
sgd_reg.fit(X, y.ravel())
C:\ProgramData\Anaconda3\lib\site-packages\sklearn\linear_model\stochastic_gradient.py:152: DeprecationWarning: n_iter parameter is deprecated in 0.19 and will be removed in 0.21. Use max_iter and tol instead.
  DeprecationWarning)
Out[42]:
SGDRegressor(alpha=1, average=False, early_stopping=False, epsilon=0.1,
       eta0=0.01, fit_intercept=True, l1_ratio=0.15,
       learning_rate='invscaling', loss='squared_loss', max_iter=None,
       n_iter=1000, n_iter_no_change=5, penalty='l2', power_t=0.25,
       random_state=42, shuffle=True, tol=None, validation_fraction=0.1,
       verbose=0, warm_start=False)
In [43]:
sgd_reg.predict([[1.5]])
Out[43]:
array([1.63288998])

4.5.2 라쏘(lasso) 회귀

. 비용 함수 : $J(\theta)={\rm MSE}(\theta)+\alpha\sum_{i=1}^n|\theta_i|$

. 특징 :

.. 덜 중요한 특성의 가중치를 완전히 제거 (가중치를 $0$)

.. 희소 모델을 만듦

. 서브그래디언트(sub-gradient) 벡터 : 미분 불가능한 점이 있을 때 ($\theta_j=0$에서 미분 불가능)

.. $g(\theta, J)=\frac{2}{m}X'(X\theta-y)+\alpha({\rm sign}(\theta_1),\ldots,{\rm sign}(\theta_n))'$

.. ${\rm sign}(t)=1$ if $t>0;$ $0$ if $t=0;$ $-1$ if $t<0$

. sklearn.linear_model.Lasso : $L_1$ 규제를 가진 선형 회귀 모델 적합

http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html

In [44]:
plt.figure(figsize=(8,4))
plt.subplot(121)
plot_model(Lasso, polynomial=False, alphas=(0, 0.1, 1), random_state=42)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.subplot(122)
plot_model(Lasso, polynomial=True, alphas=(0, 10**-7, 1), tol=1, random_state=42)
save_fig("lasso_regression_plot")
plt.show()

. Lasso 클래스와 SGD Regressor 클래스의 비교

In [45]:
lasso_reg = Lasso(alpha=0.1)
lasso_reg.fit(X, y)
Out[45]:
Lasso(alpha=0.1, copy_X=True, fit_intercept=True, max_iter=1000,
   normalize=False, positive=False, precompute=False, random_state=None,
   selection='cyclic', tol=0.0001, warm_start=False)
In [46]:
lasso_reg.predict([[1.5]])
Out[46]:
array([1.63550962])
In [47]:
# SGDRegressor로 직접 계산

# sgd_reg = SGDRegressor(max_iter=1000, penalty="l1", random_state=42)
sgd_reg = SGDRegressor(n_iter=1000, alpha=0.1, penalty="l1", random_state=42)
sgd_reg.fit(X, y.ravel())
C:\ProgramData\Anaconda3\lib\site-packages\sklearn\linear_model\stochastic_gradient.py:152: DeprecationWarning: n_iter parameter is deprecated in 0.19 and will be removed in 0.21. Use max_iter and tol instead.
  DeprecationWarning)
Out[47]:
SGDRegressor(alpha=0.1, average=False, early_stopping=False, epsilon=0.1,
       eta0=0.01, fit_intercept=True, l1_ratio=0.15,
       learning_rate='invscaling', loss='squared_loss', max_iter=None,
       n_iter=1000, n_iter_no_change=5, penalty='l1', power_t=0.25,
       random_state=42, shuffle=True, tol=None, validation_fraction=0.1,
       verbose=0, warm_start=False)
In [48]:
sgd_reg.predict([[1.5]])
Out[48]:
array([1.63189741])

4.5.3. 엘라스틱넷(elastic net)

. 릿지 회귀와 라쏘 회귀를 절충한 모델

. $J(\theta)={\rm MSE}(\theta)+r\alpha\sum_{i=1}^n|\theta_i|+(1-r)\alpha \frac{1}{2}\sum_{i=1}^n\theta_i^2$

. 모델 선택 기준 (선형 회귀, 릿지, 라쏘, 엘레스틱넷)

.. 선형 회귀(규제 없음)보다 릿지를 선호

.. 유용한 특성이 소수이면 라쏘나 엘레스틱넷

.. $m<n$이거나 특성들이 강하게 연관되어 있으면 라쏘보다 엘레스틱넷을 선호

. sklearn.linear_model.ElasticNet : $L_1$ 규제와 $L_2$ 규제를 가진 선형 회귀 모델 적합

Linear regression with combined L1 and L2 priors as regularizer.
http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html

. ElasticNet 클래스와 SGD Regressor 클래스의 비교

In [49]:
elastic_net = ElasticNet(alpha=0.1, l1_ratio=0.5, random_state=42)
elastic_net.fit(X, y)
Out[49]:
ElasticNet(alpha=0.1, copy_X=True, fit_intercept=True, l1_ratio=0.5,
      max_iter=1000, normalize=False, positive=False, precompute=False,
      random_state=42, selection='cyclic', tol=0.0001, warm_start=False)
In [50]:
elastic_net.predict([[1.5]])
Out[50]:
array([1.63813146])
In [51]:
# SGDRegressor로 직접 계산

# sgd_reg = SGDRegressor(max_iter=1000, penalty="l1", random_state=42)
sgd_reg = SGDRegressor(n_iter=1000, penalty="elasticnet", alpha=0.1, l1_ratio=0.5, random_state=42)
sgd_reg.fit(X, y.ravel())
C:\ProgramData\Anaconda3\lib\site-packages\sklearn\linear_model\stochastic_gradient.py:152: DeprecationWarning: n_iter parameter is deprecated in 0.19 and will be removed in 0.21. Use max_iter and tol instead.
  DeprecationWarning)
Out[51]:
SGDRegressor(alpha=0.1, average=False, early_stopping=False, epsilon=0.1,
       eta0=0.01, fit_intercept=True, l1_ratio=0.5,
       learning_rate='invscaling', loss='squared_loss', max_iter=None,
       n_iter=1000, n_iter_no_change=5, penalty='elasticnet', power_t=0.25,
       random_state=42, shuffle=True, tol=None, validation_fraction=0.1,
       verbose=0, warm_start=False)
In [52]:
sgd_reg.predict([[1.5]])
Out[52]:
array([1.63503003])

4.6 로지스틱 회귀

4.6.1 확률 추정

. 특징 :

.. 선형 회귀처럼 입력 특성의 가중값을 기준으로!

.. 먼저 로지스틱 값 (확률)을 출력

. 확률 추정 :

.. $\hat{p}=h_\theta(x)=\sigma(\theta'x)$

.. 로지스틱 (로짓) :

... 시그모이드 함수(sigmoid function, S자 모양)

... (0,1)의 값을 가짐

... $\sigma(t)=\frac{1}{1+\exp(-t)}$

. 예측 :

.. $\hat{y}=\left\{\begin{array}{ll}0,&\hat{p}\lt0.5\\1,&\hat{p}\ge0.5\end{array}\right.$

%E1%84%89%E1%85%B3%E1%84%8F%E1%85%B3%E1%84%85%E1%85%B5%E1%86%AB%E1%84%89%E1%85%A3%E1%86%BA%202018-09-24%20%E1%84%8B%E1%85%A9%E1%84%92%E1%85%AE%2012.37.31.png

. numpy.linspace : 명시한 구간을 등간격으로 나눔

https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.linspace.html

In [53]:
t = np.linspace(-10, 10, 100)
sig = 1 / (1 + np.exp(-t))
plt.figure(figsize=(9, 3))
plt.plot([-10, 10], [0, 0], "k-")
plt.plot([-10, 10], [0.5, 0.5], "k:")
plt.plot([-10, 10], [1, 1], "k:")
plt.plot([0, 0], [-1.1, 1.1], "k-")
plt.plot(t, sig, "b-", linewidth=2, label=r"$\sigma(t) = \frac{1}{1 + e^{-t}}$")
plt.xlabel("t")
plt.legend(loc="upper left", fontsize=20)
plt.axis([-10, 10, -0.1, 1.1])
plt.show()

4.6.2 훈련과 비용 함수

. 목적 : $y=1$에 대해서는 높은 확률, $y=0$에 대해서는 낮은 확률로 추정하는 $\theta$를 찾음

. 비용 함수 :

.. $J(\theta)=-\frac{1}{m}\sum_{i=1}^m\{y^{(i)}\log{\hat p}^{(i)}+(1-y^{(i)})\log(1-{\hat p}^{(i)})\}$

.. 훈련 샘플의 비용 함수 : $\left\{\begin{array}{ll}-\log\hat p,&y=1\\ -\log(1-\hat p),&y=0\end{array}\right.$

. 편도함수 :

.. $\frac{\partial}{\partial\theta_j}J(\theta)=\frac{1}{m}\sum_{i=1}^m\left(\sigma(\theta'x^{(i)})-y^{(i)}\right)x_j^{(i)}$

.. note that $\frac{\partial}{\partial\theta_j}{\hat p}^{(i)}={\hat p}^{(i)}(1-{\hat p}^{(i)})x_j$

4.6.3 결정 경계

. 붓꽃 자료 :

.. 세 품종 : Setosa, Versicolor, Virginica

.. 150 개 자료

.. 꽃잎과 꽃받침의 너비와 길이

. 양쪽의 확률이 0.5가 되는 점으로 결정, 즉 $\theta'x=0$을 만족하는 $x$ 값으로 결정

.. 특성이 하나일 때, $x_1=-\theta_0/\theta_1$

.. 특성이 두 개일 때, $x_1$가 주어졌다고 가정하면 $x_2=-(\theta_0+\theta_1x_1)/\theta_2$

%E1%84%89%E1%85%B3%E1%84%8F%E1%85%B3%E1%84%85%E1%85%B5%E1%86%AB%E1%84%89%E1%85%A3%E1%86%BA%202018-09-24%20%E1%84%8B%E1%85%A9%E1%84%92%E1%85%AE%201.39.41.png

. sklearn.datasets : 토이 데이터셋

http://scikit-learn.org/stable/datasets/index.html

. sklearn.datasets.load_iris : 붓꽃 자료 로드

http://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_iris.html

In [54]:
iris = datasets.load_iris()
In [55]:
list(iris.keys())
Out[55]:
['data', 'target', 'target_names', 'DESCR', 'feature_names', 'filename']
In [56]:
iris["target"]
Out[56]:
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])
In [57]:
iris["target_names"]
Out[57]:
array(['setosa', 'versicolor', 'virginica'], dtype='<U10')
In [58]:
iris["feature_names"]
Out[58]:
['sepal length (cm)',
 'sepal width (cm)',
 'petal length (cm)',
 'petal width (cm)']
In [59]:
# 꽃잎의 너비
X = iris["data"][:, 3:]
X.shape
Out[59]:
(150, 1)

. 이진 클래스로 축소해서

In [60]:
# Iris-Virginica이면 1 아니면 0
y = (iris["target"] == 2).astype(np.int)
y.shape
Out[60]:
(150,)

. sklearn.linear_model.LogisticRegression : 로지스틱 회귀 분류기

http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html

. feature 개수 : 1개

In [61]:
log_reg = LogisticRegression(random_state=42)
log_reg.fit(X, y)
C:\ProgramData\Anaconda3\lib\site-packages\sklearn\linear_model\logistic.py:433: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)
Out[61]:
LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
          intercept_scaling=1, max_iter=100, multi_class='warn',
          n_jobs=None, penalty='l2', random_state=42, solver='warn',
          tol=0.0001, verbose=0, warm_start=False)
In [62]:
# 꽃잎의 너비가 0~3cm인 꽃에 대해 모델의 추정 확률을 계산
X_new = np.linspace(0, 3, 1000).reshape(-1, 1)
y_proba = log_reg.predict_proba(X_new)
print(y_proba[:5,:])
[[0.98554411 0.01445589]
 [0.98543168 0.01456832]
 [0.98531838 0.01468162]
 [0.98520422 0.01479578]
 [0.98508919 0.01491081]]
In [63]:
# 추정 확률이 0.5 이상이 되는 최소의 꽃잎의 너비 결정
decision_boundary = X_new[y_proba[:, 1] >= 0.5][0]
In [64]:
plt.figure(figsize=(8, 3))
plt.plot(X[y==0], y[y==0], "bs")
plt.plot(X[y==1], y[y==1], "g^")
plt.plot([decision_boundary, decision_boundary], [0, 1], "k:", linewidth=2)
plt.plot(X_new, y_proba[:, 1], "g-", linewidth=2, label="Iris-Virginica")
plt.plot(X_new, y_proba[:, 0], "b--", linewidth=2, label="Not Iris-Virginica")
plt.text(decision_boundary+0.02, 0.15, "결정 경계", fontsize=14, color="k", ha="center")
plt.arrow(decision_boundary, 0.08, -0.3, 0, head_width=0.05, head_length=0.1, fc='b', ec='b')
plt.arrow(decision_boundary, 0.92, 0.3, 0, head_width=0.05, head_length=0.1, fc='g', ec='g')
plt.xlabel("꽃잎의 너비 (cm)", fontsize=14)
plt.ylabel("확률", fontsize=14)
plt.legend(loc="center left", fontsize=14)
plt.axis([0, 3, -0.02, 1.02])
save_fig("logistic_regression_plot")
plt.show()
In [65]:
decision_boundary
Out[65]:
array([1.61561562])

. 결정 경계를 직접 계산하면 ...

In [66]:
-log_reg.intercept_/log_reg.coef_
Out[66]:
array([[1.61278266]])
In [67]:
log_reg.predict([[1.7], [1.5]])
Out[67]:
array([1, 0])

. feature 개수 : 2개

In [68]:
X = iris["data"][:, (2, 3)]  # petal length, petal width
y = (iris["target"] == 2).astype(np.int)
In [69]:
# 모델 규제 강도 하이퍼파라미터 : C (클수록 규제 감소)
log_reg = LogisticRegression(C=10**10, random_state=42)
log_reg.fit(X, y)
C:\ProgramData\Anaconda3\lib\site-packages\sklearn\linear_model\logistic.py:433: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)
Out[69]:
LogisticRegression(C=10000000000, class_weight=None, dual=False,
          fit_intercept=True, intercept_scaling=1, max_iter=100,
          multi_class='warn', n_jobs=None, penalty='l2', random_state=42,
          solver='warn', tol=0.0001, verbose=0, warm_start=False)
In [70]:
x0, x1 = np.meshgrid(
        np.linspace(2.9, 7, 500).reshape(-1, 1),
        np.linspace(0.8, 2.7, 200).reshape(-1, 1),)
x0.shape,x1.shape
Out[70]:
((200, 500), (200, 500))
In [71]:
print(x0[0:2,:5])
print(x1[:5,0:2])
[[2.9        2.90821643 2.91643287 2.9246493  2.93286573]
 [2.9        2.90821643 2.91643287 2.9246493  2.93286573]]
[[0.8        0.8       ]
 [0.80954774 0.80954774]
 [0.81909548 0.81909548]
 [0.82864322 0.82864322]
 [0.83819095 0.83819095]]
In [72]:
X_new = np.c_[x0.ravel(), x1.ravel()]
print(X_new.shape)
print(X_new[:5,:])
(100000, 2)
[[2.9        0.8       ]
 [2.90821643 0.8       ]
 [2.91643287 0.8       ]
 [2.9246493  0.8       ]
 [2.93286573 0.8       ]]
In [73]:
y_proba = log_reg.predict_proba(X_new)
print(y_proba.shape)
(100000, 2)
In [74]:
plt.figure(figsize=(10, 4))
plt.plot(X[y==0, 0], X[y==0, 1], "bs") # Not Virginica
plt.plot(X[y==1, 0], X[y==1, 1], "g^") #Virginica

zz = y_proba[:, 1].reshape(x0.shape)
contour = plt.contour(x0, x1, zz, cmap=plt.cm.brg)

print(log_reg.coef_)
print("")
print(log_reg.intercept_)
print("")

left_right = np.array([2.9, 7])
boundary = -(log_reg.coef_[0][0] * left_right + log_reg.intercept_[0]) / log_reg.coef_[0][1]
print(boundary)

# Virginica라고 판단할 확률이 15%(왼쪽 아래)부터 90%(오른쪽 위) 변하는 6개의 등고선 그림
plt.clabel(contour, inline=1, fontsize=12)
plt.plot(left_right, boundary, "k--", linewidth=3)
plt.text(3.5, 1.5, "Not Iris-Virginica ", fontsize=14, color="b", ha="center")
plt.text(6.5, 2.3, "Iris-Virginica", fontsize=14, color="g", ha="center")
plt.xlabel("꽃앞의 길이", fontsize=14)
plt.ylabel("꽃잎의 너비", fontsize=14)
plt.axis([2.9, 7, 0.8, 2.7])
save_fig("logistic_regression_contour_plot")
plt.show()
[[ 5.75286199 10.44454566]]

[-45.26057652]

[2.73609573 0.47781328]

4.6.4 소프트맥스(softmax) 회귀

. 다중 클래스를 가진 로지스틱 회귀

. 다항 로지스틱 회귀라고 함

. 원리 : 샘플 $x$가 주어지면 소프트맥스 함수 계산

.. 클래스 $k$에 대한 소프트맥스 함수 ; $s_k(x)=(\theta^{(k)})'x$

.. 클래스 $k$에 속할 확률 : ${\hat p}_k=\sigma(s(x))_k=\frac{\exp(s_k(x))}{\sum_{k=1}^K\exp(s_k(x))}$

. 예측 : ${\hat y}={\rm argmax}_k\sigma(s(x))_k={\rm argmax}_ks_k(x)={\rm argmax}_k((\theta^{(k)})'x)$

. Cross-entropy 비용 함수 :

.. $J(\Theta)=-\frac{1}{m}\sum_{i=1}^m\sum_{k=1}^Ky_k^{(i)}\log{\hat p}_k^{(i)}$

.. $i$ 번째 샘플에 대한 target 클래스가 $k$일 때, $y_k^{(i)}=1$이고 그 외에는 $0$

. 클래스 $k$에 대한 cross-entropy의 경사 벡터 : $\nabla_{\theta^{(k)}}J(\Theta)=\frac{1}{m}\sum_{i=1}^m ({\hat p}_k^{(i)}-y_k^{(i)})x^{(i)}$

. sklearn에서는 "multi_class = mutinomial", "solver = lbfgs"로 조정

. 다중(3개) 클래스로 하여

In [75]:
# 꽃잎 길이, 꽃잎 넓이
X = iris["data"][:, (2, 3)]  
y = iris["target"]
In [76]:
softmax_reg = LogisticRegression(multi_class="multinomial",solver="lbfgs", C=10, random_state=42)
softmax_reg.fit(X, y)
Out[76]:
LogisticRegression(C=10, class_weight=None, dual=False, fit_intercept=True,
          intercept_scaling=1, max_iter=100, multi_class='multinomial',
          n_jobs=None, penalty='l2', random_state=42, solver='lbfgs',
          tol=0.0001, verbose=0, warm_start=False)

. matplotlib.pyplot.contourf

Draw filled contours.
https://matplotlib.org/2.1.0/api/_as_gen/matplotlib.pyplot.contourf.html

In [77]:
x0, x1 = np.meshgrid(
        np.linspace(0, 8, 500).reshape(-1, 1),
        np.linspace(0, 3.5, 200).reshape(-1, 1),)
X_new = np.c_[x0.ravel(), x1.ravel()]
X_new[:5,]
Out[77]:
array([[0.        , 0.        ],
       [0.01603206, 0.        ],
       [0.03206413, 0.        ],
       [0.04809619, 0.        ],
       [0.06412826, 0.        ]])
In [78]:
y_proba = softmax_reg.predict_proba(X_new)
y_predict = softmax_reg.predict(X_new)

print(y_proba[:, 1].shape)
print(y_predict.shape)
(100000,)
(100000,)
In [79]:
zz1 = y_proba[:, 1].reshape(x0.shape) #ver
zz = y_predict.reshape(x0.shape)

plt.figure(figsize=(10, 4))
plt.plot(X[y==2, 0], X[y==2, 1], "g^", label="Iris-Virginica")
plt.plot(X[y==1, 0], X[y==1, 1], "bs", label="Iris-Versicolor")
plt.plot(X[y==0, 0], X[y==0, 1], "yo", label="Iris-Setosa")

custom_cmap = ListedColormap(['#fafab0','#9898ff','#a0faa0'])

plt.contourf(x0, x1, zz, cmap=custom_cmap)
contour = plt.contour(x0, x1, zz1, cmap=plt.cm.brg)
plt.clabel(contour, inline=1, fontsize=12)
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.legend(loc="center left", fontsize=14)
plt.axis([0, 7, 0, 3.5])
save_fig("softmax_regression_contour_plot")
plt.show()
In [80]:
softmax_reg.predict([[5, 2]])
Out[80]:
array([2])
In [81]:
softmax_reg.predict_proba([[5, 2]])
Out[81]:
array([[6.38014896e-07, 5.74929995e-02, 9.42506362e-01]])

The End of Chapter 4