로지스틱회귀 파라메터 추정
02 Jul 2017 | regression
이번 글에서는 로지스틱회귀 모델의 계수를 추정하는 방법을 살펴보도록 하겠습니다. 이번 글은 고려대 김성범 교수님 강의와 ‘밑바닥부터 시작하는 데이터과학(조엘 그루스 지음, 인사이트 펴냄)’을 정리하였음을 먼저 밝힙니다. 로지스틱회귀에 대한 기본적인 내용은 이곳을 참고하시면 좋을 것 같습니다. 그럼 시작하겠습니다.
최대우도추정법
로지스틱회귀 모델 계수를 추정하기 위해서는 우선 최대우도추정법(Maximum Likelihood Estimation) 개념을 먼저 짚고 넘어가야 합니다. 가령 데이터가 임의의 파라메터 $θ$에 의존하는 확률분포를 따르고 표본데이터 $v_1, v_2,…,v_n$이 주어졌다고 해보겠습니다.
우리는 데이터가 발생할 확률 $P(v_1,…,v_n$|$θ)$를 구하고 싶은데, $θ$를 모르는 상황입니다. 따라서 우리는 베이즈 규칙을 활용해 다음과 같이 $θ$가 발생할 우도(likelihood)로 바꿔서 생각할 수 있습니다. 다시 말해 조건절과 결과절을 뒤집어 원하는 값을 구하는 것입니다.
\(L(\theta |{ v }_{ 1 },{ v }_{ 2 },...,{ v }_{ n })\)
이 경우 가장 적절한 $θ$는 우도를 최대화해주는 값이 될 것입니다. 다시 말해 관측된 데이터가 발생할 경우를 가장 높게 만들어주는 값이라는 의미입니다.
전제 : 베르누이 분포
로지스틱 회귀는 베르누이 시행(Bernoulli trial)을 전제로 하는 모델입니다. 베르누이 시행이란 어떤 실험이 두 가지 결과만을 가지는 실험을 가리킵니다. 베르누이 시행의 결과에 따라 0(실패) 또는 1(성공)의 값을 대응시키는 확률변수(random variable)를 베르누이 확률변수라 합니다. 이 확률변수의 확률분포를 베르누이 분포라고 합니다. 베르누이 확률변수 Y의 분포는 아래 표와 같습니다.
$Y$
0
1
$P(Y=y_i)$
$1- p$
$p$
위 표를 수식으로 정리하면 아래와 같습니다.
\[P(Y={ y }_{ i })={ p }^{ { y }_{ i } }{ (1-p) }^{ 1-{ y }_{ i } }\quad ({ y }_{ i }=0,1)\]
베르누이 확률변수 Y에 관한 우도함수(likelihood function)는 아래와 같습니다.
\[L=\prod _{ i }^{ }{ { p }^{ { y }_{ i } }{ (1-p) }^{ 1-{ y }_{ i } }}\]
로지스틱회귀의 우도함수
학습데이터에 관측치가 $i$개가 있고, 정답 범주가 2개(1 또는 0)뿐인 이항로지스틱 모델의 파라메터 $β$가 주어졌다고 가정해 보겠습니다. 그러면 $i$번째 관측치의 종속변수 $y_i$는 $σ(β^Tx_i)$의 확률로 1, $1-σ(β^Tx_i)$의 확률로 0이 됩니다. 여기에서 $x_i$는 $i$번째 관측치의 독립변수, $σ$는 로지스틱 함수를 가리킵니다. 따라서 로지스틱회귀의 우도함수는 다음과 같이 쓸 수 있습니다.
\[L=\prod _{ i }^{ }{ { \sigma ({ \beta }^{ T }\overrightarrow { { x }_{ i } } ) }^{ { y }_{ i } }{ \left\{ 1-\sigma ({ \beta }^{ T }\overrightarrow { { x }_{ i } } ) \right\} }^{ 1-{ y }_{ i } } }\]
로지스틱 회귀의 파라메터 $β$는 앞서 언급한 MLE로 구합니다. 로그는 단조 증가함수이므로 로그 우도함수(log-likelihood function)를 최대로 하는 회귀계수 $β$는 동시에 우도를 최대화하는 $β$이며 그 역도 성립합니다. 로그 우도함수는 아래와 같이 정리할 수 있습니다.
\[\ln { L } =\sum _{ i }^{ }{ { y }_{ i }\ln { \left\{ \sigma ({ \beta }^{ T }{ \overrightarrow { x_{ i } } }) \right\} } } +\sum _{ i }^{ }{ \left( 1-{ y }_{ i } \right) \ln { \left\{ 1-\sigma ({ \beta }^{ T }{ \overrightarrow { x_{ i } } }) \right\} } }\]
다만 위 로그 우도함수는 추정 대상 파라메터인 회귀계수 $β$에 대해 비선형이기 때문에 선형회귀와 같이 명시적인 해가 존재하지 않습니다. 따라서 Stochastic Gradient Descent(SGD) 같은 반복적이고 점진적인 방식으로 해를 구하게 됩니다.
Stochastic Gradient Descent
SGD는 반복문을 돌 때마다 개별 데이터 포인트에 대한 그래디언트를 계산하고 이 그래디언트의 반대 방향으로 파라메터를 업데이트해 함수의 최소값을 구하는 기법입니다. 로지스틱 회귀의 파라메터는 명시적인 해가 존재하지 않기 때문에 이처럼 값을 조금씩 바꿔가면서 최적해를 찾아가는 것입니다.
그런데 우리는 로지스틱 회귀의 로그 우도를 최대화하는 파라메터 $β$를 찾으려는 것이므로 기존 SGD의 반대 방향으로 업데이트를 수행해야 할 것입니다. SGD 관련 메인 파이썬 코드는 다음과 같습니다.
def minimize_stochastic(target_fn, gradient_fn, x, y, theta_0, alpha_0=0.01):
# SGD 방식으로 gradient descent
# minimize_batch보다 훨씬 빠르다
data = zip(x, y)
# theta_0를 초기 theta로
theta = theta_0
# alpha_0를 초기 이동거리(step_size)로
alpha = alpha_0
# 시작할 때의 최소값
min_theta, min_value = None, float("inf")
iterations_with_no_improvement = 0
# 만약 100번 넘게 반복하는 동안 value가 더 작아지지 않으면 멈춤
while iterations_with_no_improvement < 100:
value = sum( target_fn(x_i, y_i, theta) for x_i, y_i in data )
# 새로운 최솟값을 찾았다면
if value < min_value:
# 이 값을 저장
min_theta, min_value = theta, value
# 100번 카운트도 초기화
iterations_with_no_improvement = 0
# 기본 이동거리로 돌아감
alpha = alpha_0
# 만약 최솟값이 줄어들지 않는다면
else:
# 이동거리 축소
alpha *= 0.9
# 100번 카운트에 1을 더함
iterations_with_no_improvement += 1
# 반복문이 돌 때마다 in_random_order를 호출하기 때문에
# 매 iter마다 그래디언트를 계산하는 순서가 달라짐
for x_i, y_i in in_random_order(data):
# 각 데이터 포인트에 대해 그래디언트를 계산
gradient_i = gradient_fn(x_i, y_i, theta)
# 기존 theta에서, 학습률(alpha)과 그래디언트를 뺀 것을 업데이트
theta = vector_subtract(theta, scalar_multiply(alpha, gradient_i))
return min_theta
# maximize_stochastic = -minimize_stochastic
def maximize_stochastic(target_fn, gradient_fn, x, y, theta_0, alpha_0=0.01):
return minimize_stochastic(negate(target_fn),
negate_all(gradient_fn),
x, y, theta_0, alpha_0)
# x를 입력하면 -f(x)를 반환해주는 함수
def negate(f):
return lambda *args, **kwargs: -f(*args, **kwargs)
# f가 여러 숫자를 반환할 때 모든 숫자를 음수로 변환
def negate_all(f):
return lambda *args, **kwargs: [-y for y in f(*args, **kwargs)]
그러면 이제는 target_fn과 gradient_fn, theta를 정의해야 합니다. target_fn은 목적함수로 로그 우도를 아웃풋으로 산출하는 함수를 의미합니다. 로지스틱회귀 모델의 로그 우도함수는 다음과 같습니다.
import math
# 로지스틱 함수
def logistic(x):
return 1.0 / (1 + math.exp(-x))
# 로지스틱 함수의 1차 도함수
def logistic_prime(x):
return logistic(x) * (1 - logistic(x))
# i번째 관측치에 대한 로그우도
def logistic_log_likelihood_i(x_i, y_i, beta):
# y_i의 레이블이 1이면
if y_i == 1:
return math.log(logistic(dot(x_i, beta)))
# y_i의 레이블이 0이면
else:
return math.log(1 - logistic(dot(x_i, beta)))
# 데이터 전체에 대한 로그우도 (전체 합)
def logistic_log_likelihood(x, y, beta):
return sum(logistic_log_likelihood_i(x_i, y_i, beta)
for x_i, y_i in zip(x, y))
‘gradient_fn’은 각 파라메터에 대한 목적함수의 그래디언트를 가리킵니다. 미분 공식을 활용해 로그우도 함수를 $β$에 대해 해석적으로 미분한 식을 코드로 표현한 것입니다.
def logistic_log_partial_ij(x_i, y_i, beta, j):
"""here i is the index of the data point,
j the index of the derivative"""
return (y_i - logistic(dot(x_i, beta))) * x_i[j]
def logistic_log_gradient_i(x_i, y_i, beta):
"""the gradient of the log likelihood
corresponding to the i-th data point"""
return [logistic_log_partial_ij(x_i, y_i, beta, j)
for j, _ in enumerate(beta)]
def logistic_log_gradient(x, y, beta):
return reduce(vector_add,
[logistic_log_gradient_i(x_i, y_i, beta)
for x_i, y_i in zip(x,y)])
마지막으로 ‘theta’는 해당 목적함수의 파라메터인데요, 우리 문제에선 $β$를 말합니다. 만약 독립변수가 2개라면 상수항을 포함해 3차원 벡터가 될 것입니다. 이밖에 위 코드 실행에 필요한 함수는 다음과 같습니다.
def in_random_order(data):
# Stochastic Gradient Descent 수행을 위한 함수로,
# 한번 반복문을 돌 때마다 임의의 순서로 데이터 포인트를 반환
# 데이터 포인트의 인덱스를 list로 생성
indexes = [i for i, _ in enumerate(data)]
# 이 인덱스를 랜덤하게 섞는다
random.shuffle(indexes)
# 이 순서대로 데이터를 반환한다
for i in indexes:
yield data[i]
def vector_subtract(v, w):
"""subtracts two vectors componentwise"""
return [v_i - w_i for v_i, w_i in zip(v,w)]
def scalar_multiply(c, v):
return [c * v_i for v_i in v]
def dot(v, w):
"""v_1 * w_1 + ... + v_n * w_n"""
return sum(v_i * w_i for v_i, w_i in zip(v, w))
분석 예시
분석 대상 데이터는 다음과 같습니다. 독립변수는 2개, 독립변수는 범주형 1개(0 또는 1)입니다.
data = [(0.7,48000,1),(1.9,48000,0),(2.5,60000,1),(4.2,63000,0),(6,76000,0),(6.5,69000,0),(7.5,76000,0),(8.1,88000,0),(8.7,83000,1),(10,83000,1),(0.8,43000,0),(1.8,60000,0),(10,79000,1),(6.1,76000,0),(1.4,50000,0),(9.1,92000,0),(5.8,75000,0),(5.2,69000,0),(1,56000,0),(6,67000,0),(4.9,74000,0),(6.4,63000,1),(6.2,82000,0),(3.3,58000,0),(9.3,90000,1),(5.5,57000,1),(9.1,102000,0),(2.4,54000,0),(8.2,65000,1),(5.3,82000,0),(9.8,107000,0),(1.8,64000,0),(0.6,46000,1),(0.8,48000,0),(8.6,84000,1),(0.6,45000,0),(0.5,30000,1),(7.3,89000,0),(2.5,48000,1),(5.6,76000,0),(7.4,77000,0),(2.7,56000,0),(0.7,48000,0),(1.2,42000,0),(0.2,32000,1),(4.7,56000,1),(2.8,44000,1),(7.6,78000,0),(1.1,63000,0),(8,79000,1),(2.7,56000,0),(6,52000,1),(4.6,56000,0),(2.5,51000,0),(5.7,71000,0),(2.9,65000,0),(1.1,33000,1),(3,62000,0),(4,71000,0),(2.4,61000,0),(7.5,75000,0),(9.7,81000,1),(3.2,62000,0),(7.9,88000,0),(4.7,44000,1),(2.5,55000,0),(1.6,41000,0),(6.7,64000,1),(6.9,66000,1),(7.9,78000,1),(8.1,102000,0),(5.3,48000,1),(8.5,66000,1),(0.2,56000,0),(6,69000,0),(7.5,77000,0),(8,86000,0),(4.4,68000,0),(4.9,75000,0),(1.5,60000,0),(2.2,50000,0),(3.4,49000,1),(4.2,70000,0),(7.7,98000,0),(8.2,85000,0),(5.4,88000,0),(0.1,46000,0),(1.5,37000,0),(6.3,86000,0),(3.7,57000,0),(8.4,85000,0),(2,42000,0),(5.8,69000,1),(2.7,64000,0),(3.1,63000,0),(1.9,48000,0),(10,72000,1),(0.2,45000,0),(8.6,95000,0),(1.5,64000,0),(9.8,95000,0),(5.3,65000,0),(7.5,80000,0),(9.9,91000,0),(9.7,50000,1),(2.8,68000,0),(3.6,58000,0),(3.9,74000,0),(4.4,76000,0),(2.5,49000,0),(7.2,81000,0),(5.2,60000,1),(2.4,62000,0),(8.9,94000,0),(2.4,63000,0),(6.8,69000,1),(6.5,77000,0),(7,86000,0),(9.4,94000,0),(7.8,72000,1),(0.2,53000,0),(10,97000,0),(5.5,65000,0),(7.7,71000,1),(8.1,66000,1),(9.8,91000,0),(8,84000,0),(2.7,55000,0),(2.8,62000,0),(9.4,79000,0),(2.5,57000,0),(7.4,70000,1),(2.1,47000,0),(5.3,62000,1),(6.3,79000,0),(6.8,58000,1),(5.7,80000,0),(2.2,61000,0),(4.8,62000,0),(3.7,64000,0),(4.1,85000,0),(2.3,51000,0),(3.5,58000,0),(0.9,43000,0),(0.9,54000,0),(4.5,74000,0),(6.5,55000,1),(4.1,41000,1),(7.1,73000,0),(1.1,66000,0),(9.1,81000,1),(8,69000,1),(7.3,72000,1),(3.3,50000,0),(3.9,58000,0),(2.6,49000,0),(1.6,78000,0),(0.7,56000,0),(2.1,36000,1),(7.5,90000,0),(4.8,59000,1),(8.9,95000,0),(6.2,72000,0),(6.3,63000,0),(9.1,100000,0),(7.3,61000,1),(5.6,74000,0),(0.5,66000,0),(1.1,59000,0),(5.1,61000,0),(6.2,70000,0),(6.6,56000,1),(6.3,76000,0),(6.5,78000,0),(5.1,59000,0),(9.5,74000,1),(4.5,64000,0),(2,54000,0),(1,52000,0),(4,69000,0),(6.5,76000,0),(3,60000,0),(4.5,63000,0),(7.8,70000,0),(3.9,60000,1),(0.8,51000,0),(4.2,78000,0),(1.1,54000,0),(6.2,60000,0),(2.9,59000,0),(2.1,52000,0),(8.2,87000,0),(4.8,73000,0),(2.2,42000,1),(9.1,98000,0),(6.5,84000,0),(6.9,73000,0),(5.1,72000,0),(9.1,69000,1),(9.8,79000,1),]
전체 코드를 구동하는 명령은 다음과 같습니다.
import random
data = map(list, data) # change tuples to lists
x = [[1] + row[:2] for row in data] # each element is [1, experience, salary]
y = [row[2] for row in data] # each element is paid_account
beta_0 = [1, 1, 1]
beta_hat = maximize_stochastic(logistic_log_likelihood_i,
logistic_log_gradient_i,
x, y, beta_0)
이번 글에서는 로지스틱회귀 모델의 계수를 추정하는 방법을 살펴보도록 하겠습니다. 이번 글은 고려대 김성범 교수님 강의와 ‘밑바닥부터 시작하는 데이터과학(조엘 그루스 지음, 인사이트 펴냄)’을 정리하였음을 먼저 밝힙니다. 로지스틱회귀에 대한 기본적인 내용은 이곳을 참고하시면 좋을 것 같습니다. 그럼 시작하겠습니다.
최대우도추정법
로지스틱회귀 모델 계수를 추정하기 위해서는 우선 최대우도추정법(Maximum Likelihood Estimation) 개념을 먼저 짚고 넘어가야 합니다. 가령 데이터가 임의의 파라메터 $θ$에 의존하는 확률분포를 따르고 표본데이터 $v_1, v_2,…,v_n$이 주어졌다고 해보겠습니다.
우리는 데이터가 발생할 확률 $P(v_1,…,v_n$|$θ)$를 구하고 싶은데, $θ$를 모르는 상황입니다. 따라서 우리는 베이즈 규칙을 활용해 다음과 같이 $θ$가 발생할 우도(likelihood)로 바꿔서 생각할 수 있습니다. 다시 말해 조건절과 결과절을 뒤집어 원하는 값을 구하는 것입니다. \(L(\theta |{ v }_{ 1 },{ v }_{ 2 },...,{ v }_{ n })\) 이 경우 가장 적절한 $θ$는 우도를 최대화해주는 값이 될 것입니다. 다시 말해 관측된 데이터가 발생할 경우를 가장 높게 만들어주는 값이라는 의미입니다.
전제 : 베르누이 분포
로지스틱 회귀는 베르누이 시행(Bernoulli trial)을 전제로 하는 모델입니다. 베르누이 시행이란 어떤 실험이 두 가지 결과만을 가지는 실험을 가리킵니다. 베르누이 시행의 결과에 따라 0(실패) 또는 1(성공)의 값을 대응시키는 확률변수(random variable)를 베르누이 확률변수라 합니다. 이 확률변수의 확률분포를 베르누이 분포라고 합니다. 베르누이 확률변수 Y의 분포는 아래 표와 같습니다.
| $Y$ | 0 | 1 |
|---|---|---|
| $P(Y=y_i)$ | $1- p$ | $p$ |
위 표를 수식으로 정리하면 아래와 같습니다.
\[P(Y={ y }_{ i })={ p }^{ { y }_{ i } }{ (1-p) }^{ 1-{ y }_{ i } }\quad ({ y }_{ i }=0,1)\]베르누이 확률변수 Y에 관한 우도함수(likelihood function)는 아래와 같습니다.
\[L=\prod _{ i }^{ }{ { p }^{ { y }_{ i } }{ (1-p) }^{ 1-{ y }_{ i } }}\]로지스틱회귀의 우도함수
학습데이터에 관측치가 $i$개가 있고, 정답 범주가 2개(1 또는 0)뿐인 이항로지스틱 모델의 파라메터 $β$가 주어졌다고 가정해 보겠습니다. 그러면 $i$번째 관측치의 종속변수 $y_i$는 $σ(β^Tx_i)$의 확률로 1, $1-σ(β^Tx_i)$의 확률로 0이 됩니다. 여기에서 $x_i$는 $i$번째 관측치의 독립변수, $σ$는 로지스틱 함수를 가리킵니다. 따라서 로지스틱회귀의 우도함수는 다음과 같이 쓸 수 있습니다.
\[L=\prod _{ i }^{ }{ { \sigma ({ \beta }^{ T }\overrightarrow { { x }_{ i } } ) }^{ { y }_{ i } }{ \left\{ 1-\sigma ({ \beta }^{ T }\overrightarrow { { x }_{ i } } ) \right\} }^{ 1-{ y }_{ i } } }\]로지스틱 회귀의 파라메터 $β$는 앞서 언급한 MLE로 구합니다. 로그는 단조 증가함수이므로 로그 우도함수(log-likelihood function)를 최대로 하는 회귀계수 $β$는 동시에 우도를 최대화하는 $β$이며 그 역도 성립합니다. 로그 우도함수는 아래와 같이 정리할 수 있습니다.
\[\ln { L } =\sum _{ i }^{ }{ { y }_{ i }\ln { \left\{ \sigma ({ \beta }^{ T }{ \overrightarrow { x_{ i } } }) \right\} } } +\sum _{ i }^{ }{ \left( 1-{ y }_{ i } \right) \ln { \left\{ 1-\sigma ({ \beta }^{ T }{ \overrightarrow { x_{ i } } }) \right\} } }\]다만 위 로그 우도함수는 추정 대상 파라메터인 회귀계수 $β$에 대해 비선형이기 때문에 선형회귀와 같이 명시적인 해가 존재하지 않습니다. 따라서 Stochastic Gradient Descent(SGD) 같은 반복적이고 점진적인 방식으로 해를 구하게 됩니다.
Stochastic Gradient Descent
SGD는 반복문을 돌 때마다 개별 데이터 포인트에 대한 그래디언트를 계산하고 이 그래디언트의 반대 방향으로 파라메터를 업데이트해 함수의 최소값을 구하는 기법입니다. 로지스틱 회귀의 파라메터는 명시적인 해가 존재하지 않기 때문에 이처럼 값을 조금씩 바꿔가면서 최적해를 찾아가는 것입니다.
그런데 우리는 로지스틱 회귀의 로그 우도를 최대화하는 파라메터 $β$를 찾으려는 것이므로 기존 SGD의 반대 방향으로 업데이트를 수행해야 할 것입니다. SGD 관련 메인 파이썬 코드는 다음과 같습니다.
def minimize_stochastic(target_fn, gradient_fn, x, y, theta_0, alpha_0=0.01):
# SGD 방식으로 gradient descent
# minimize_batch보다 훨씬 빠르다
data = zip(x, y)
# theta_0를 초기 theta로
theta = theta_0
# alpha_0를 초기 이동거리(step_size)로
alpha = alpha_0
# 시작할 때의 최소값
min_theta, min_value = None, float("inf")
iterations_with_no_improvement = 0
# 만약 100번 넘게 반복하는 동안 value가 더 작아지지 않으면 멈춤
while iterations_with_no_improvement < 100:
value = sum( target_fn(x_i, y_i, theta) for x_i, y_i in data )
# 새로운 최솟값을 찾았다면
if value < min_value:
# 이 값을 저장
min_theta, min_value = theta, value
# 100번 카운트도 초기화
iterations_with_no_improvement = 0
# 기본 이동거리로 돌아감
alpha = alpha_0
# 만약 최솟값이 줄어들지 않는다면
else:
# 이동거리 축소
alpha *= 0.9
# 100번 카운트에 1을 더함
iterations_with_no_improvement += 1
# 반복문이 돌 때마다 in_random_order를 호출하기 때문에
# 매 iter마다 그래디언트를 계산하는 순서가 달라짐
for x_i, y_i in in_random_order(data):
# 각 데이터 포인트에 대해 그래디언트를 계산
gradient_i = gradient_fn(x_i, y_i, theta)
# 기존 theta에서, 학습률(alpha)과 그래디언트를 뺀 것을 업데이트
theta = vector_subtract(theta, scalar_multiply(alpha, gradient_i))
return min_theta
# maximize_stochastic = -minimize_stochastic
def maximize_stochastic(target_fn, gradient_fn, x, y, theta_0, alpha_0=0.01):
return minimize_stochastic(negate(target_fn),
negate_all(gradient_fn),
x, y, theta_0, alpha_0)
# x를 입력하면 -f(x)를 반환해주는 함수
def negate(f):
return lambda *args, **kwargs: -f(*args, **kwargs)
# f가 여러 숫자를 반환할 때 모든 숫자를 음수로 변환
def negate_all(f):
return lambda *args, **kwargs: [-y for y in f(*args, **kwargs)]
그러면 이제는 target_fn과 gradient_fn, theta를 정의해야 합니다. target_fn은 목적함수로 로그 우도를 아웃풋으로 산출하는 함수를 의미합니다. 로지스틱회귀 모델의 로그 우도함수는 다음과 같습니다.
import math
# 로지스틱 함수
def logistic(x):
return 1.0 / (1 + math.exp(-x))
# 로지스틱 함수의 1차 도함수
def logistic_prime(x):
return logistic(x) * (1 - logistic(x))
# i번째 관측치에 대한 로그우도
def logistic_log_likelihood_i(x_i, y_i, beta):
# y_i의 레이블이 1이면
if y_i == 1:
return math.log(logistic(dot(x_i, beta)))
# y_i의 레이블이 0이면
else:
return math.log(1 - logistic(dot(x_i, beta)))
# 데이터 전체에 대한 로그우도 (전체 합)
def logistic_log_likelihood(x, y, beta):
return sum(logistic_log_likelihood_i(x_i, y_i, beta)
for x_i, y_i in zip(x, y))
‘gradient_fn’은 각 파라메터에 대한 목적함수의 그래디언트를 가리킵니다. 미분 공식을 활용해 로그우도 함수를 $β$에 대해 해석적으로 미분한 식을 코드로 표현한 것입니다.
def logistic_log_partial_ij(x_i, y_i, beta, j):
"""here i is the index of the data point,
j the index of the derivative"""
return (y_i - logistic(dot(x_i, beta))) * x_i[j]
def logistic_log_gradient_i(x_i, y_i, beta):
"""the gradient of the log likelihood
corresponding to the i-th data point"""
return [logistic_log_partial_ij(x_i, y_i, beta, j)
for j, _ in enumerate(beta)]
def logistic_log_gradient(x, y, beta):
return reduce(vector_add,
[logistic_log_gradient_i(x_i, y_i, beta)
for x_i, y_i in zip(x,y)])
마지막으로 ‘theta’는 해당 목적함수의 파라메터인데요, 우리 문제에선 $β$를 말합니다. 만약 독립변수가 2개라면 상수항을 포함해 3차원 벡터가 될 것입니다. 이밖에 위 코드 실행에 필요한 함수는 다음과 같습니다.
def in_random_order(data):
# Stochastic Gradient Descent 수행을 위한 함수로,
# 한번 반복문을 돌 때마다 임의의 순서로 데이터 포인트를 반환
# 데이터 포인트의 인덱스를 list로 생성
indexes = [i for i, _ in enumerate(data)]
# 이 인덱스를 랜덤하게 섞는다
random.shuffle(indexes)
# 이 순서대로 데이터를 반환한다
for i in indexes:
yield data[i]
def vector_subtract(v, w):
"""subtracts two vectors componentwise"""
return [v_i - w_i for v_i, w_i in zip(v,w)]
def scalar_multiply(c, v):
return [c * v_i for v_i in v]
def dot(v, w):
"""v_1 * w_1 + ... + v_n * w_n"""
return sum(v_i * w_i for v_i, w_i in zip(v, w))
분석 예시
분석 대상 데이터는 다음과 같습니다. 독립변수는 2개, 독립변수는 범주형 1개(0 또는 1)입니다.
data = [(0.7,48000,1),(1.9,48000,0),(2.5,60000,1),(4.2,63000,0),(6,76000,0),(6.5,69000,0),(7.5,76000,0),(8.1,88000,0),(8.7,83000,1),(10,83000,1),(0.8,43000,0),(1.8,60000,0),(10,79000,1),(6.1,76000,0),(1.4,50000,0),(9.1,92000,0),(5.8,75000,0),(5.2,69000,0),(1,56000,0),(6,67000,0),(4.9,74000,0),(6.4,63000,1),(6.2,82000,0),(3.3,58000,0),(9.3,90000,1),(5.5,57000,1),(9.1,102000,0),(2.4,54000,0),(8.2,65000,1),(5.3,82000,0),(9.8,107000,0),(1.8,64000,0),(0.6,46000,1),(0.8,48000,0),(8.6,84000,1),(0.6,45000,0),(0.5,30000,1),(7.3,89000,0),(2.5,48000,1),(5.6,76000,0),(7.4,77000,0),(2.7,56000,0),(0.7,48000,0),(1.2,42000,0),(0.2,32000,1),(4.7,56000,1),(2.8,44000,1),(7.6,78000,0),(1.1,63000,0),(8,79000,1),(2.7,56000,0),(6,52000,1),(4.6,56000,0),(2.5,51000,0),(5.7,71000,0),(2.9,65000,0),(1.1,33000,1),(3,62000,0),(4,71000,0),(2.4,61000,0),(7.5,75000,0),(9.7,81000,1),(3.2,62000,0),(7.9,88000,0),(4.7,44000,1),(2.5,55000,0),(1.6,41000,0),(6.7,64000,1),(6.9,66000,1),(7.9,78000,1),(8.1,102000,0),(5.3,48000,1),(8.5,66000,1),(0.2,56000,0),(6,69000,0),(7.5,77000,0),(8,86000,0),(4.4,68000,0),(4.9,75000,0),(1.5,60000,0),(2.2,50000,0),(3.4,49000,1),(4.2,70000,0),(7.7,98000,0),(8.2,85000,0),(5.4,88000,0),(0.1,46000,0),(1.5,37000,0),(6.3,86000,0),(3.7,57000,0),(8.4,85000,0),(2,42000,0),(5.8,69000,1),(2.7,64000,0),(3.1,63000,0),(1.9,48000,0),(10,72000,1),(0.2,45000,0),(8.6,95000,0),(1.5,64000,0),(9.8,95000,0),(5.3,65000,0),(7.5,80000,0),(9.9,91000,0),(9.7,50000,1),(2.8,68000,0),(3.6,58000,0),(3.9,74000,0),(4.4,76000,0),(2.5,49000,0),(7.2,81000,0),(5.2,60000,1),(2.4,62000,0),(8.9,94000,0),(2.4,63000,0),(6.8,69000,1),(6.5,77000,0),(7,86000,0),(9.4,94000,0),(7.8,72000,1),(0.2,53000,0),(10,97000,0),(5.5,65000,0),(7.7,71000,1),(8.1,66000,1),(9.8,91000,0),(8,84000,0),(2.7,55000,0),(2.8,62000,0),(9.4,79000,0),(2.5,57000,0),(7.4,70000,1),(2.1,47000,0),(5.3,62000,1),(6.3,79000,0),(6.8,58000,1),(5.7,80000,0),(2.2,61000,0),(4.8,62000,0),(3.7,64000,0),(4.1,85000,0),(2.3,51000,0),(3.5,58000,0),(0.9,43000,0),(0.9,54000,0),(4.5,74000,0),(6.5,55000,1),(4.1,41000,1),(7.1,73000,0),(1.1,66000,0),(9.1,81000,1),(8,69000,1),(7.3,72000,1),(3.3,50000,0),(3.9,58000,0),(2.6,49000,0),(1.6,78000,0),(0.7,56000,0),(2.1,36000,1),(7.5,90000,0),(4.8,59000,1),(8.9,95000,0),(6.2,72000,0),(6.3,63000,0),(9.1,100000,0),(7.3,61000,1),(5.6,74000,0),(0.5,66000,0),(1.1,59000,0),(5.1,61000,0),(6.2,70000,0),(6.6,56000,1),(6.3,76000,0),(6.5,78000,0),(5.1,59000,0),(9.5,74000,1),(4.5,64000,0),(2,54000,0),(1,52000,0),(4,69000,0),(6.5,76000,0),(3,60000,0),(4.5,63000,0),(7.8,70000,0),(3.9,60000,1),(0.8,51000,0),(4.2,78000,0),(1.1,54000,0),(6.2,60000,0),(2.9,59000,0),(2.1,52000,0),(8.2,87000,0),(4.8,73000,0),(2.2,42000,1),(9.1,98000,0),(6.5,84000,0),(6.9,73000,0),(5.1,72000,0),(9.1,69000,1),(9.8,79000,1),]
전체 코드를 구동하는 명령은 다음과 같습니다.
import random
data = map(list, data) # change tuples to lists
x = [[1] + row[:2] for row in data] # each element is [1, experience, salary]
y = [row[2] for row in data] # each element is paid_account
beta_0 = [1, 1, 1]
beta_hat = maximize_stochastic(logistic_log_likelihood_i,
logistic_log_gradient_i,
x, y, beta_0)
Comments