In [68]:
from __future__ import division
from collections import Counter
from functools import partial
from linear_algebra import dot, vector_add
from stats import median, standard_deviation, de_mean
from probability import normal_cdf
from gradient_descent import minimize_stochastic
#from simple_linear_regression import total_sum_of_squares
import math, random
In [69]:
def total_sum_of_squares(y):
"""the total squared variation of y_i's from their mean"""
return sum(v ** 2 for v in de_mean(y))
14장에서 다뤘던 모델
$y_i=\alpha+\beta x_i+\epsilon_i$
여기에 독립 변수를 추가하면
$\Rightarrow$시간(분)= $\alpha+\beta_1*$(친구 수)$+\beta_2*$(근무 시간)$+\beta_3$*(박사 학위 취득 여부)+$\epsilon$
각 입력값 $x_i$가 숫자 하나가 아니라 $k$개의 숫자인 $x_i1,\cdot\cdot\cdot,x_ik$라고 한다면,
다중 회귀 모델은 다음과 같은 형태를 띈다.
$y_i=\alpha+\beta_1 x_{i1}+\cdot\cdot\cdot+\beta_k x_{ik}+\epsilon_i$
다중 회귀 분석에서는 보통 파라미터 벡터를 $\beta$라고 부름
여기에 상수항 $\alpha$까지 덧붙이려면, 각 데이터 x_i의 앞부분에 1을 덧붙이면 된다.
beta = [alpha, beta_1,$\cdot\cdot\cdot$, x_ik]
그리고 각 데이터는 다음과 같이 된다.
x_i = [1, x_i1, $\cdot\cdot\cdot$, x_ik]
이렇게 하면 모델을 다음과 같이 나타낼 수 있다.
In [70]:
def predict(x_i, beta):
"""각 x_i의 첫 번째 항목은 1이라고 가정"""
return dot(x_i, beta)
독립 변수 x는 다음과 같은 벡터들의 열로 표현할 수 있다.
[1, # 상수항
49, # 친구의 수
4, # 하루 근무 시간
0] # 박사 학위 취득 여부오류 함수
In [71]:
def error(x_i, y_i, beta):
return y_i - predict(x_i, beta)
SGD를 사용하기 위한 오류 제곱 값
In [72]:
def squared_error(x_i, y_i, beta):
return error(x_i, y_i, beta) ** 2
만약 미적분을 알고 있다면, 오류를 직접 계산할 수도 있다.
In [73]:
def squared_error_gradient(x_i, y_i, beta):
"""i번째 오류 제곱 값의 beta에 대한 기울기"""
return [-2 * x_ij * error(x_i, y_i, beta)
for x_ij in x_i]
SGD를 사용해서 최적의 베타를 계산
In [74]:
def estimate_beta(x, y):
beta_initial = [random.random() for x_i in x[0]]
return minimize_stochastic(squared_error,
squared_error_gradient,
x, y,
beta_initial,
0.001)
실험 데이터 셋팅 (14장과 동일)
In [75]:
x = [[1,49,4,0],[1,41,9,0],[1,40,8,0],[1,25,6,0],[1,21,1,0],[1,21,0,0],[1,19,3,0],[1,19,0,0],[1,18,9,0],[1,18,8,0],[1,16,4,0],[1,15,3,0],[1,15,0,0],[1,15,2,0],[1,15,7,0],[1,14,0,0],[1,14,1,0],[1,13,1,0],[1,13,7,0],[1,13,4,0],[1,13,2,0],[1,12,5,0],[1,12,0,0],[1,11,9,0],[1,10,9,0],[1,10,1,0],[1,10,1,0],[1,10,7,0],[1,10,9,0],[1,10,1,0],[1,10,6,0],[1,10,6,0],[1,10,8,0],[1,10,10,0],[1,10,6,0],[1,10,0,0],[1,10,5,0],[1,10,3,0],[1,10,4,0],[1,9,9,0],[1,9,9,0],[1,9,0,0],[1,9,0,0],[1,9,6,0],[1,9,10,0],[1,9,8,0],[1,9,5,0],[1,9,2,0],[1,9,9,0],[1,9,10,0],[1,9,7,0],[1,9,2,0],[1,9,0,0],[1,9,4,0],[1,9,6,0],[1,9,4,0],[1,9,7,0],[1,8,3,0],[1,8,2,0],[1,8,4,0],[1,8,9,0],[1,8,2,0],[1,8,3,0],[1,8,5,0],[1,8,8,0],[1,8,0,0],[1,8,9,0],[1,8,10,0],[1,8,5,0],[1,8,5,0],[1,7,5,0],[1,7,5,0],[1,7,0,0],[1,7,2,0],[1,7,8,0],[1,7,10,0],[1,7,5,0],[1,7,3,0],[1,7,3,0],[1,7,6,0],[1,7,7,0],[1,7,7,0],[1,7,9,0],[1,7,3,0],[1,7,8,0],[1,6,4,0],[1,6,6,0],[1,6,4,0],[1,6,9,0],[1,6,0,0],[1,6,1,0],[1,6,4,0],[1,6,1,0],[1,6,0,0],[1,6,7,0],[1,6,0,0],[1,6,8,0],[1,6,4,0],[1,6,2,1],[1,6,1,1],[1,6,3,1],[1,6,6,1],[1,6,4,1],[1,6,4,1],[1,6,1,1],[1,6,3,1],[1,6,4,1],[1,5,1,1],[1,5,9,1],[1,5,4,1],[1,5,6,1],[1,5,4,1],[1,5,4,1],[1,5,10,1],[1,5,5,1],[1,5,2,1],[1,5,4,1],[1,5,4,1],[1,5,9,1],[1,5,3,1],[1,5,10,1],[1,5,2,1],[1,5,2,1],[1,5,9,1],[1,4,8,1],[1,4,6,1],[1,4,0,1],[1,4,10,1],[1,4,5,1],[1,4,10,1],[1,4,9,1],[1,4,1,1],[1,4,4,1],[1,4,4,1],[1,4,0,1],[1,4,3,1],[1,4,1,1],[1,4,3,1],[1,4,2,1],[1,4,4,1],[1,4,4,1],[1,4,8,1],[1,4,2,1],[1,4,4,1],[1,3,2,1],[1,3,6,1],[1,3,4,1],[1,3,7,1],[1,3,4,1],[1,3,1,1],[1,3,10,1],[1,3,3,1],[1,3,4,1],[1,3,7,1],[1,3,5,1],[1,3,6,1],[1,3,1,1],[1,3,6,1],[1,3,10,1],[1,3,2,1],[1,3,4,1],[1,3,2,1],[1,3,1,1],[1,3,5,1],[1,2,4,1],[1,2,2,1],[1,2,8,1],[1,2,3,1],[1,2,1,1],[1,2,9,1],[1,2,10,1],[1,2,9,1],[1,2,4,1],[1,2,5,1],[1,2,0,1],[1,2,9,1],[1,2,9,1],[1,2,0,1],[1,2,1,1],[1,2,1,1],[1,2,4,1],[1,1,0,1],[1,1,2,1],[1,1,2,1],[1,1,5,1],[1,1,3,1],[1,1,10,1],[1,1,6,1],[1,1,0,1],[1,1,8,1],[1,1,6,1],[1,1,4,1],[1,1,9,1],[1,1,9,1],[1,1,4,1],[1,1,2,1],[1,1,9,1],[1,1,0,1],[1,1,8,1],[1,1,6,1],[1,1,1,1],[1,1,1,1],[1,1,5,1]]
num_friends_good = [49,41,40,25,21,21,19,19,18,18,16,15,15,15,15,14,14,13,13,13,13,12,12,11,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,8,8,8,8,8,8,8,8,8,8,8,8,8,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
daily_minutes_good = [68.77,51.25,52.08,38.36,44.54,57.13,51.4,41.42,31.22,34.76,54.01,38.79,47.59,49.1,27.66,41.03,36.73,48.65,28.12,46.62,35.57,32.98,35,26.07,23.77,39.73,40.57,31.65,31.21,36.32,20.45,21.93,26.02,27.34,23.49,46.94,30.5,33.8,24.23,21.4,27.94,32.24,40.57,25.07,19.42,22.39,18.42,46.96,23.72,26.41,26.97,36.76,40.32,35.02,29.47,30.2,31,38.11,38.18,36.31,21.03,30.86,36.07,28.66,29.08,37.28,15.28,24.17,22.31,30.17,25.53,19.85,35.37,44.6,17.23,13.47,26.33,35.02,32.09,24.81,19.33,28.77,24.26,31.98,25.73,24.86,16.28,34.51,15.23,39.72,40.8,26.06,35.76,34.76,16.13,44.04,18.03,19.65,32.62,35.59,39.43,14.18,35.24,40.13,41.82,35.45,36.07,43.67,24.61,20.9,21.9,18.79,27.61,27.21,26.61,29.77,20.59,27.53,13.82,33.2,25,33.1,36.65,18.63,14.87,22.2,36.81,25.53,24.62,26.25,18.21,28.08,19.42,29.79,32.8,35.99,28.32,27.79,35.88,29.06,36.28,14.1,36.63,37.49,26.9,18.58,38.48,24.48,18.95,33.55,14.24,29.04,32.51,25.63,22.22,19,32.73,15.16,13.9,27.2,32.01,29.27,33,13.74,20.42,27.32,18.23,35.35,28.48,9.08,24.62,20.12,35.26,19.92,31.02,16.49,12.16,30.7,31.22,34.65,13.13,27.51,33.2,31.57,14.1,33.42,17.44,10.12,24.42,9.82,23.39,30.93,15.03,21.67,31.09,33.29,22.61,26.89,23.48,8.38,27.81,32.35,23.84]
In [76]:
random.seed(0)
beta = estimate_beta(x, daily_minutes_good)
In [77]:
beta # 분 = 30.63 + 0.972 친구 수 - 1.868 근무 시간 + 0.911 박사 학위 취득 여부
Out[77]:
모델의 계수는 해당 항목의 영향력을 나타낸다.
친구 수 증가 => 대략 1분 증가
근무 시간 증가 => 대략 2분 감소
박사 학위 취득 => 대략 1분 증가
이러한 해석은 변수간의 관계를 직접적으로 설명해 주지 못한다.
예를 들어, 친구의 수가 다른 사용자들의 근무 시간은 서로 다를 수 있다.
이 모델은 이러한 관계를 잡아내지 못한다.
이러한 문제는 친구 수와 근무 시간을 곱한 새로운 변수로 해결할 수 있다.
변수가 점점 추가되기 시작하면, 각 계수가 유의미한지 살펴봐야 한다.
변수끼리 곱한 값, 변수의 log값, 변수의 제곱 값 등 수 많은 변수를 추가할 수 있기 때문이다.
In [78]:
def multiple_r_squared(x, y, beta):
sum_of_squared_errors = sum(error(x_i, y_i, beta) ** 2
for x_i, y_i in zip(x, y))
return 1.0 - sum_of_squared_errors / total_sum_of_squares(y)
In [79]:
multiple_r_squared(x, daily_minutes_good, beta)
Out[79]:
회귀 분석 모델에 새로운 변수를 추가하면 R 제곱 값이 어쩔 수 없이 증가한다.
따라서 다중 회귀 분석 모델은 언제나 단순 회귀 분석 모델보다 작은 오류를 갖게 된다.
이러한 이유로 다중 회귀 분석 모델에서는 각 계수의 표준 오차를 살펴 봐야 한다.
계수의 표준 오차는 추정된 $\beta_1$의 계수가 얼마나 확실한지 알려준다.
표준오차 : 모집단 전체를 알 수 없는 상황에서, 모집단이 정상분포(정규분포)를 이루고 있다는 가정 하에, 여러번의 샘플링을 통해 각 표본 집단의 평균들로 이루어진 표준 평균 분포를 얻고, 이 표준 평균 분포의 표준 편차가 표준 오차가 된다.
표준오차는 모평균과 표본평균 사이의 오차를 알려주므로, 모집단의 표준편차가 클수록 표준오차도 커지고, 사례가 많을 수록 작아진다.
오차를 측정하기 위해서는 각 오류 $\epsilon_1$는 독립이며, 평균은 0이고 표준편차는 $\sigma$인 정규분포의 확률변수라는 가정이 필요하다.
표준오차가 클수록 해당 계수는 무의미해 진다.
In [80]:
def bootstrap_sample(data):
"""len(data)개의 항목을 중복을 허용한 무작위 추출"""
return [random.choice(data) for _ in data]
In [81]:
def bootstrap_statistic(data, stats_fn, num_samples):
"""num_samples개의 bootstrap 샘플에 대해 stats_fn을 적용"""
return [stats_fn(bootstrap_sample(data))
for _ in range(num_samples)]
예를 들어, 다음과 같은 두 가지 데이터를 살펴보자.
In [82]:
# 101개의 데이터가 모두 100에 인접
close_to_100 = [99.5 + random.random() for _ in range(101)]
# 101개의 데이터 중 50개는 0에 인접, 50개는 200에 인접
far_from_100 = ([99.5 + random.random()] +
[random.random() for _ in range(50)] +
[200 + random.random() for _ in range(50)])
In [83]:
print(close_to_100)
print(far_from_100)
만약 두 데이터의 중앙값을 계산해 보면 둘 다 대략 100에 가까운 것을 확인할 수 있다.
In [84]:
print(median(close_to_100))
print(median(far_from_100))
하지만 다음과 같이 bootstrap을 적용해 보면,
close_to_100은 100에 대부분 가깝지만,
far_from_100은 0 또는 200에 가까운 것을 확인할 수 있다.
In [85]:
print(bootstrap_statistic(close_to_100, median, 100))
In [86]:
print(bootstrap_statistic(far_from_100, median, 100))
첫번째는 표준편차가 0에 가깝지만, 두번째는 표준편차가 100에 가까운 것을 확인할 수 있다.
In [87]:
standard_deviation(bootstrap_statistic(close_to_100, median, 100))
Out[87]:
In [88]:
standard_deviation(bootstrap_statistic(far_from_100, median, 100))
Out[88]:
데이터가 이렇게 극단적인 경우에는 데이터를 직접 살펴보면 문제를 쉽게 파악할 수 있지만,
대부분의 경우 데이터만 살펴보는 것으로는 부족하다.
In [89]:
def estimate_sample_beta(sample):
x_sample, y_sample = zip(*sample) # magic unzipping trick
return estimate_beta(x_sample, y_sample)
In [90]:
random.seed(0) # 예시와 동일한 결과를 얻기 위해 설정
In [91]:
bootstrap_betas = bootstrap_statistic(list(zip(x, daily_minutes_good)),
estimate_sample_beta,
100)
In [92]:
bootstrap_betas
Out[92]:
그리고 각 계수의 표준 오차를 추정할 수 있다.
In [93]:
bootstrap_standard_errors = [
standard_deviation([beta[i] for beta in bootstrap_betas])
for i in range(4)
]
In [94]:
bootstrap_standard_errors
Out[94]:
이제 '과연 $\beta_1$는 0일까?' 같은 가설을 검증해 볼 수 있다.
p-value(유의 확률)은 귀무가설이 맞을 경우 대립가설 쪽의 값이 나올 확률을 나타내는 값. 확률 값이라고도 한다. 표본 평균이 귀무가설 값에서 멀수록 작아지게 된다.
In [95]:
def p_value(beta_hat_j, sigma_hat_j):
if beta_hat_j > 0:
return 2 * (1 - normal_cdf(beta_hat_j / sigma_hat_j))
else:
return 2 * normal_cdf(beta_hat_j / sigma_hat_j)
In [96]:
print(p_value(30.63, 1.174))
print(p_value(0.972, 0.079))
print(p_value(-1.868, 0.131))
print(p_value(0.911, 0.990))
대부분의 계수들은 0이 아닌 것으로 검증되었으나,
박사 학위 취득 여부에 대한 계수에 의미가 없을 수 있다는 것을 암시
In [97]:
# alpha는 패널티의 강도를 조절하는 하이퍼 파라미터
# 보통 "lamda"라고 표현하지만 파이썬에서는 이미 사용 중인 키워드이다.
def ridge_penalty(beta, alpha):
return alpha * dot(beta[1:], beta[1:])
def squared_error_ridge(x_i, y_i, beta, alpha):
"""beta를 사용할 때 오류와 패널티의 합을 추정"""
return error(x_i, y_i, beta) ** 2 + ridge_penalty(beta, alpha)
그리고 이전과 동일하게 경사 하강법을 적용할 수 있다.
In [98]:
def ridge_penalty_gradient(beta, alpha):
"""패널티의 기울기"""
return [0] + [2 * alpha * beta_j for beta_j in beta[1:]]
def squared_error_ridge_gradient(x_i, y_i, beta, alpha):
"""i번 오류 제곱 값과 패널티의 기울기"""
return vector_add(squared_error_gradient(x_i, y_i, beta),
ridge_penalty_gradient(beta, alpha))
def estimate_beta_ridge(x, y, alpha):
"""패널티가 alpha인 리지 회귀를 경사 하강법으로 학습"""
beta_initial = [random.random() for x_i in x[0]]
return minimize_stochastic(partial(squared_error_ridge, alpha=alpha),
partial(squared_error_ridge_gradient,
alpha=alpha),
x, y,
beta_initial,
0.001)
만약 alpha가 0이라면 패널티는 전혀 없으며, 이전과 동일한 모델이 학습될 것이다.
In [99]:
random.seed(0)
beta_0 = estimate_beta_ridge(x, daily_minutes_good, alpha=0.0)
beta_0
Out[99]:
In [100]:
dot(beta_0[1:], beta_0[1:])
Out[100]:
In [101]:
multiple_r_squared(x, daily_minutes_good, beta_0)
Out[101]:
그리고 alpha를 증가시킬수록 적합성은 감소하고, beta의 크기도 감소한다.
In [102]:
beta_0_01 = estimate_beta_ridge(x, daily_minutes_good, alpha=0.01)
print(beta_0_01)
print(dot(beta_0_01[1:], beta_0_01[1:]))
print(multiple_r_squared(x, daily_minutes_good, beta_0_01))
In [103]:
beta_0_1 = estimate_beta_ridge(x, daily_minutes_good, alpha=0.1)
print(beta_0_1)
print(dot(beta_0_1[1:], beta_0_1[1:]))
print(multiple_r_squared(x, daily_minutes_good, beta_0_1))
In [104]:
beta_1 = estimate_beta_ridge(x, daily_minutes_good, alpha=1)
print(beta_1)
print(dot(beta_1[1:], beta_1[1:]))
print(multiple_r_squared(x, daily_minutes_good, beta_1))
In [105]:
beta_10 = estimate_beta_ridge(x, daily_minutes_good, alpha=10)
print(beta_10)
print(dot(beta_10[1:], beta_10[1:]))
multiple_r_squared(x, daily_minutes_good, beta_10)
Out[105]:
패널티가 증가하면 박사 학위 취득 여부 변수는 사라진다.
다른 형태의 패널티를 사용하는 lasso regression도 있다.
In [106]:
def lasso_penalty(beta, alpha):
return alpha * sum(abs(beta_i) for beta_i in beta[1:])
리지 회귀의 패널티는 총 계수의 합을 줄여 주지만,
라쏘 회귀의 패널티는 모든 계수를 최대한 0으로 만들어 주며,
보다 희소한(sparse) 모델을 학습하게 해준다.
하지만 라쏘 회귀는 경사 하강법으로는 학습할 수 없다.