로지스틱 회귀(Logistic Regression)

$p$개의 특징을 가지는 벡터 $x$와 가중치 $\omega$에 대한 다중 선형 회귀의 일반 식은 다음과 같이 쓸 수 있습니다.

 

$y=\omega_{0}+\omega_{1}x_{1}+\cdots+\omega_{p}x_{p}$

 

하지만 위와 같은 선형식은 범주형 종속변수를 표현하기에 적절하지 않습니다.

 

그래서 다음과 같이 좌변을 종속변수의 값이 1이 될 확률로 두고 식을 세워 보겠습니다.

 

$P(y=1|x;\omega) = \omega_{0}+\omega_{1}x_{1}+\cdots+\omega_{p}x_{p}$ = $\omega^{T}x$   $( x = [x_{1},x_{2},\cdots,x_{p}])$

                        

그런데, 좌변은 확률 값이기 때문에 [0,1] 구간의 값을 가지지만 우변은 $[-\infty,\infty]$구간의 값을 가지므로

 

식이 성립하지 않는 경우가 생깁니다.

 

이 문제를 해결하기 위해 아래와 같이 좌변에 확률 대신 승산(odds)을 사용해 봅시다.

 

$\frac{ P(y=1|x;\omega) }{ 1-P(y=1|x;\omega) }= \omega^{T}x$

 

그러나 좌변은 $[0,\infty]$구간의 값을 가지므로 여전히 스케일이 맞지 않습니다.

 

그래서 좌변에 자연로그를 취해봅니다. 이 과정을 로짓 변환이라고 합니다.

 

$log(\frac{ P(y=1|x;\omega) }{ 1-P(y=1|x;\omega) })=\omega^{T}x$

 

드디어 좌변과 우변 모두 $[-\infty,\infty]$구간의 값을 가지게 되었습니다.

 

 

 

이제는 아래와 같은 과정을 통해 식을 변형해 보겠습니다.

 

잠시 $P(y=1|x;\omega) = p$,   $\omega^{T}x=k$ 로 두겠습니다.

 

그러면 위의 식을 $log(\frac{ p }{ 1-p })=k$와 같이 쓸 수 있습니다.

 

그리고 다음 과정을 거쳐서...

$\frac{ p }{ 1-p }=e^{k}$

$p=e^{k}(1-p)=e^{k}-pe^{k}$

$p(1+e^{k})=e^{k}$

$p=\frac{ e^{k} }{ 1+e^{k} }=\frac{ 1 }{1+ e^{-k} }$

 

최종적으로 아래와 같은 식이 만들어집니다.

 

$P(y=1|x;\omega)=\frac{ 1 }{ 1+e^{-\omega^{T}x} }$

 

이 때 $\sigma(z)=\frac{ 1 }{ 1+e^{-z} }$와 같은 형태의 식이 관찰되는데, 이를 로지스틱 함수

 

혹은 시그모이드 함수라고 부르며 아래와 같은 그래프로 표현할 수 있습니다.

 

시그모이드 함수는 z값에 관계없이 항상 0과 1 사이의 값이 출력 되기때문에 확률밀도함수의 요건을 충족시킵니다.

 

 

계속 이어서, 다음 식은 y가 1일 확률을 나타냅니다.

 

$P(y=1|x;\omega)=\frac{ 1 }{ 1+e^{-\omega^{T}x} }=\sigma(\omega^{T}x)$

 

그리고 y는 0과 1 둘 중 하나의 값을 가지기 때문에 y가 0일 확률은 아래처럼 나타낼 수 있습니다.

 

$P(y=0|x;\omega)=1-P(y=1|x;\omega)$

 

따라서 다음이 성립합니다.

 

$P(y|x;\omega)=(\sigma(\omega^{T}x))^{y}(1-\sigma(\omega^{T}x))^{1-y}$

 

확률변수 y가 베르누이 분포를 따른다는 사실을 확인할 수 있습니다.

 

이제는 m개의 데이터에 대한 likelihood 함수를 구해보겠습니다.

 

$\mathcal{L}(\omega)$

 

=$P(y^{(1)},\cdots,y^{(m)}|x^{(1)},\cdots,x^{(m)};\omega)$

 

=$\prod_{ i=1 }^{ m } {P(y^{(i)}|x^{(i)};\omega  })$

 

=$\prod_{ i=1 }^{ m } {\sigma(\omega^{T}x^{(i)}))^{y^{(i)}}(1-\sigma(\omega^{T}x^{(i)}))^{1-y^{(i)}}}$

 

여기서 양변에 로그를 취해주면 각 확률 밀도 함수에 대한 곱을 합으로 표현해줄 수 있습니다.

 

$l(\omega)=log\mathcal{L}(\omega)=\sum_{ i=1 }^{ m } { y^{(i)}log\sigma(\omega^{T}x^{(i)})+(1-y^{(i)})log(1-\sigma(\omega^{T}x^{(i)})) }$

 

이때, 로그 가능도가 최대가 되는 ω를 찾아야 하므로 우리의 목표는 아래의 최적화 문제를 푸는 것으로 귀결됩니다.

 

$\hat{\omega}=arg\max\limits_{\omega} {l(\omega)}$

 

이제, 경사 하강법을 이용하여 문제를 해결해 보겠습니다.

 

우선 아래와 같은 gradient vector를 구해야 합니다. 

 

$\nabla l(\omega)=\begin{bmatrix} {\frac{ \partial l(\omega) }{ \partial\omega_{1} }} \\ \vdots \\{\frac{ \partial l(\omega) }{ \partial\omega_{n} }}\end{bmatrix}$

 

그러기 위해서 데이터 하나에 대한 편미분을 먼저 해 봅시다.

 

 

 

참고로, 유도과정 중에 시그모이드 함수의 도함수가 다음과 같음이 이용되었습니다. 

 

 

 

 

따라서 m개 데이터에 대한 편미분은 아래와 같고

 

 

 

 

이를 행렬 계산을 위해 벡터로 표현해 주면 다음과 같이 정리됩니다.

 

 

 

 

 

 

그리고 경사 하강법을 이용하여 ω벡터를 구해주면 됩니다.

 

 

 

 

 

이제, 위의 과정을 코드로 옮겨서 문제를 풀어보겠습니다.

# n 개의 표본 생성

n = 1000
m = 3
X = np.hstack([ 40+70*np.random.rand(n,1), 80+80*np.random.rand(n,1)])
y = np.zeros(n)
for i in range(n):
  if 2*X[i][0] + X[i][1] > 300:
    y[i] = 1
  else:
    y[i] = 0

dfx = pd.DataFrame(X)
dfy = pd.DataFrame(y)
df = pd.concat([dfy,dfx],axis=1)
df.columns = ['심장병 유무','체중', '평균혈압']
df=df.astype(int)
df

이것은 제가 임의로 만든 데이터이며, 체중과 혈압을 통해서 심장병의 유무를 예측하는 문제를 풀어보겠습니다.

 

# 데이터 표준화

scaler = StandardScaler()
X = scaler.fit_transform(X)
X = np.hstack([np.ones([n,1]),X])
y = y.reshape(n,1)

우선 데이터 표준화를 거치고, bias항을 추가해줍니다.

 

#sigmoid 함수 정의

def sigmoid(z):
  return 1/(1 + np.exp(-z))

sigmoid 함수를 정의합니다.

# 경사하강법 적용

def gradient_descent(X):
  alpha = 0.001
  w = np.zeros([m,1])
  for i in range(10000):
    del_lw = X.T @ (y - sigmoid(X @ w))
    w = w - alpha * (-del_lw)
  return w

 

편미분 자체를 scipy를 이용하거나 반복문을 이용하여 numerical 하게 해 줄 수도 있지만, 행렬 수준에서 수학적인 유도로 구해주었기 때문에 Gradient descent의 코드가 간결해졌습니다. 함수의 최댓값을 찾아야 하기 때문에 -부호를 꼭 붙여줍니다.

 

w = gradient_descent(X)

최적화된 가중치 벡터를 리턴 받았습니다.

 

y_predict = np.zeros(n)

for i in range(n):
  if sigmoid(-X @ w)[i] > 0.5:
    y_predict[i] = 1
  else:
     y_predict[i] = 0

데이터 행렬 X와 리턴 받은 가중치 벡터의 행렬곱 결과를 시그모이드 함수에 대입하여 결과를 예측(분류) 합니다.

 

count = 0
for i in range(n):
  if y[i]==y_predict[i]:
    count=count+1

accuracy = (count / n) *100

print('accuracy: '+str(accuracy)+'%')

accuracy: 99.8%

 

 

 

 

 

 

 

추후에 다음 내용도 업데이트하도록 하겠습니다. 

 

• Likelihood와 MLE(최대 우도 추정법)에 대한 설명
• 시그모이드 함수의 미분 결과에 대한 증명
• accuracy를 더 엄밀하고 정확하게 구하는 방법
• 경사 하강법의 적절한 하이퍼 파라미터(학습률, 사이클 횟수 등)를 정하는 방법
• 사용된 함수들을 매소드로 가지는 클래스 구현
• 사이킷런 라이브러리와의 결과 비교
• 데이터 시각화

 

 

머신러닝을 공부한 지 얼마 되지 않아 내용에 미흡한 점이 많겠지만, 고쳐야 할 내용이 보인다면

지적해 주시면 감사하겠습니다.