다운로드
작성자: admin 작성일시: 2016-10-06 21:55:10 조회수: 3159 다운로드: 155
카테고리: 머신 러닝 태그목록:

로지스틱 회귀 분석

로지스틱 회귀(Logistic Regression) 분석은 회귀 분석이라는 명칭을 가지고 있지만 분류(classsification) 방법의 일종이다.

로지스틱 회귀 모형에서는 베르누이 확률 변수(Bernoilli random variable)의 모수(parameter) $\theta$가 독립 변수 $x$에 의존한다고 가정한다.

$$ p(y \mid x, \theta) = \text{Ber} (y \mid \theta(x) )$$

여기에서 모수 $\theta$ 는 0과 1사이의 실수이며 다음과 같이 $x$의 값에 의존하는 함수이다.

$$ \theta = f(w^Tx) $$

시그모이드 함수

모수 $\theta$는 일반적인 회귀 분석의 종속 변수와 달리 0 부터 1까지의 실수값만 가질 수 있기 때문에 시그모이드 함수(sigmoid function)이라 불리는 특별한 형태의 함수 $f$를 사용해야 한다.

시그모이드 함수는 종속 변수의 모든 실수 값에 대해 유한한 구간 $(a,b)$ 사이의 한정된(bounded) 값과 양의 기울기를 가지는 함수를 말하며 다음과 같은 함수들이 주로 사용된다.

  • 로지스틱 함수 (Logistic Function)
$$ \text{logitstic}(z) = \dfrac{1}{1+\exp{(-z)}} $$
  • 오차 함수 (Error Function)
$$ \text{erf}(z) = \frac{2}{\sqrt\pi}\int_0^z e^{-t^2}\,dt $$
  • 하이퍼볼릭 탄젠트 함수 (Hyperbolic tangent)
$$ \tanh(z) = \frac{\sinh z}{\cosh z} = \frac {e^z - e^{-z}} {e^z + e^{-z}} $$
  • 역 탄젠트 함수 (Arc-tangent)
$$ \arctan(z) = \tan^{-1}(z) $$
In:
library(NORMT3)

x <- seq(-10, 10, length=1000)
plot(x, (1/(1 + exp(-x)))*2-1, type="l", xlim=c(-10, 10), ylim=c(-1.2, 1.2), ylab="", main="Sigmoids")
grid()
lines(x, erf(0.5 * sqrt(pi) * x), lty=4)
lines(x, tanh(x), lty=3)
legend("topleft", c("logistic", "erf", "tanh"), lty=c(1, 3, 4))
#plt.plot(, label="erf (scaled)")
#plt.plot(xx, np.tanh(xx), label="tanh")
#plt.ylim([-1.1, 1.1])
#plt.legend(loc=2)

로지스틱 함수

여러가지 시그모이드 중 로지스틱 함수는 다음과 같은 물리적인 의미를 부여할 수 있기 때문에 많이 사용된다.

우선 Bernoulli 시도에서 1이 나올 확률 $\theta$ 과 0이 나올 확률 $1-\theta$ 의 비(ratio)는 다음과 같은 수식이 되며 odds ratio 라고 한다.

$$ \text{odds ratio} = \dfrac{\theta}{1-\theta} $$

이 odds ratio 를 로그 변환한 것이 로지트 함수(Logit function)이다.

$$ z = \text{logit}(\text{odds ratio}) = \log \left(\dfrac{\theta}{1-\theta}\right) $$

로지스틱 함수(Logistic function) 는 이 로지트 함수의 역함수이다.

$$ \text{logitstic}(z) = \theta(z) = \dfrac{1}{1+\exp{(-z)}} $$

로지스틱 모형의 모수 추정

로지스틱 모형은 일종의 비선형 회귀 모형이지만 다음과 같이 MLE(Maximum Likelihood Estimation) 방법으로 모수 $w$를 추정할 수 있다.

여기에서는 종속 변수 $y$가 베르누이 확률 변수라고 가정한다.

$$ p(y \mid x, \theta) = \text{Ber} (y \mid \theta(x) )$$

데이터 표본이 $\{ x_i, y_i \}$일 경우 Log Likelihood $\text{LL}$ 를 구하면 다음과 같다.

$$ \begin{eqnarray} \text{LL} &=& \log \prod_{i=1}^N \theta_i(x_i)^{y_i} (1-\theta_i(x_i))^{1-y_i} \\ &=& \sum_{i=1}^N \left( y_i \log\theta_i(x_i) + (1-y_i)\log(1-\theta_i(x_i)) \right) \\ \end{eqnarray} $$

$\theta$가 로지스틱 함수 형태로 표현된다면

$$ \log \left(\dfrac{\theta(x)}{1-\theta(x)}\right) = w^T x $$$$ \theta(x) = \dfrac{1}{1 + \exp{(-w^Tx)}} $$

가 되고 이를 Log Likelihood 에 적용하면 다음과 같다.

$$ \begin{eqnarray} \text{LL} &=& \sum_{i=1}^N \left( y_i \log\theta_i(x_i) + (1-y_i)\log(1-\theta_i(x_i)) \right) \\ &=& \sum_{i=1}^N \left( y_i \log\left(\dfrac{1}{1 + \exp{(-w^Tx_i)}}\right) - (1-y_i)\log\left(\dfrac{\exp{(-w^Tx_i)}}{1 + \exp{(-w^Tx_i)}}\right) \right) \\ \end{eqnarray} $$

이 값의 최대화하는 값을 구하기 위해 chain rule를 사용하여 $w$로 미분해야 한다.

우선 $\theta$를 $w$로 미분하면

$$ \dfrac{\partial \theta}{\partial w} = \dfrac{\partial}{\partial w} \dfrac{1}{1 + \exp{(-w^Tx)}} \ = \dfrac{\exp{(-w^Tx)}}{(1 + \exp{(-w^Tx)})^2} x \ = \theta(1-\theta) x $$

chain rule를 적용하면

$$ \begin{eqnarray} \dfrac{\partial \text{LL}}{\partial w} &=& \sum_{i=1}^N \left( y_i \dfrac{1}{\theta_i(x_i;w)} - (1-y_i)\dfrac{1}{1-\theta_i(x_i;w)} \right) \dfrac{\partial \theta}{\partial w} \\ &=& \sum_{i=1}^N \big( y_i (1-\theta_i(x_i;w)) - (1-y_i)\theta_i(x_i;w) \big) x_i \\ &=& \sum_{i=1}^N \big( y_i - \theta_i(x_i;w) \big) x_i \\ \end{eqnarray} $$

이 값은 $w$에 대한 비선형 함수이므로 선형 모형과 같이 간단하게 그레디언트가 0이 되는 모수 $w$ 값에 대한 수식을 구할 수 없으며 수치적인 최적화 방법(numerical optimization)을 통해 최적 모수 $w$의 값을 구해야 한다.

수치적 최적화

단순한 Steepest Gradient 방법을 사용한다면 최적화 알고리즘은 다음과 같다.

그레디언트 벡터는 $$ g_k = \dfrac{d}{dw}(-LL) $$

이 방향으로 step size $\eta_k$ 만큼 움직이면 다음과 같이 반복적으로 최적 모수값을 구할 수 있다.

$$ \begin{eqnarray} w_{k+1} &=& w_{k} - \eta_k g_k \\ &=& w_{k} + \eta_k \sum_{i=1}^N \big( y_i - \theta_i(x_i) \big) x_i\\ \end{eqnarray} $$

R의 로지스틱 회귀

R은 로지스틱 회귀를 위한 glm 명령을 제공한다. glm은 일반 선형 회귀 모형(Generalized Linear Model)의 약자이다. 사실 로지스틱 회귀는 일반 선형 회귀 모형 중 link 함수가 logit 함수인 특수한 경우이다.

glm(formula, family=binomial(link="logit"), data)

예제 1: Michelin and Zagat 가이드 비교

다음 데이터는 뉴욕시의 레스토랑에 대한 두 개의 가이드북에서 발취한 것이다.

  • Food: Zagat Survey 2006 의 고객 평가 점수
  • InMichelin: 해당 고객 평가 점수를 받은 레스토랑 중 2006 Michelin Guide New York City 에 실린 레스토랑의 수
  • NotInMichelin: 해당 고객 평가 점수를 받은 레스토랑 중 2006 Michelin Guide New York City 에 실리지 않은 레스토랑의 수
  • mi: 해당 고객 평가 점수를 받은 레스토랑의 수
  • proportion: 해당 고객 평가 점수를 받은 레스토랑 중 2006 Michelin Guide New York City 에 실린 레스토랑의 비율
In:
df1 <- read.table("~/data/sheather/MichelinFood.txt", header=T)
df1
FoodInMichelinNotInMichelinmiproportion
15 0 1 1 0.00
16 0 1 1 0.00
17 0 8 8 0.00
18 2 13 15 0.13
19 5 13 18 0.28
20 8 25 33 0.24
21 15 11 26 0.58
22 4 8 12 0.33
23 12 6 18 0.67
24 6 1 7 0.86
25 11 1 12 0.92
26 1 1 2 0.50
27 6 1 7 0.86
28 4 0 4 1.00
In:
plot(proportion ~ Food, data=df1, cex=3, xlim=c(10, 30))
In:
model1 <- glm(proportion ~ Food, binomial(link="logit"), data=df1)
summary(model1)
Warning message in eval(expr, envir, enclos):
“non-integer #successes in a binomial glm!”
Call:
glm(formula = proportion ~ Food, family = binomial(link = "logit"), 
    data = df)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.84164  -0.29190  -0.02041   0.32627   0.43466  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -10.3600     5.2111  -1.988   0.0468 *
Food          0.4671     0.2345   1.991   0.0464 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 8.7727  on 13  degrees of freedom
Residual deviance: 1.8342  on 12  degrees of freedom
AIC: 15.388

Number of Fisher Scoring iterations: 5
In:
plot(proportion ~ Food, data=df1, cex=3, xlim=c(10, 30))
Xnew <- data.frame(Food=seq(10, 30, length=100))
lines(Xnew$Food, predict(model1, Xnew, type="response"), lwd=2)

예제 2: Michelin 가이드 예측

다음 데이터는 뉴욕시의 개별 레스토랑의 고객 평가 점수와 Michelin 가이드 수록 여부를 보인 것이다.

  • InMichelin: Michelin 가이드 수록 여부
  • Restaurant Name: 레스토랑 이름
  • Food: 식사에 대한 고객 평가 점수 (1~30)
  • Decor: 인테리어에 대한 고객 평가 점수 (1~30)
  • Service: 서비스에 대한 고객 평가 점수 (1~30)
  • Price: 저녁 식사 가격 (US$)
In:
df2 <- read.csv("~/data/sheather/MichelinNY.csv", header=T)
head(df2)
InMichelinRestaurant.NameFoodDecorServicePrice
0 14 Wall Street19 20 19 50
0 212 17 17 16 43
0 26 Seats 23 17 21 35
1 44 19 23 16 52
0 A 23 12 19 24
0 A.O.C. 18 17 17 36
In:
plot(jitter(InMichelin, 1, 0.03) ~ jitter(Food, 1, 0.1), data=df2, 
     xlim=c(10, 30), xlab="Food", ylab="InMichelin")
grid()
In:
model2 <- glm(InMichelin ~ Food, binomial(link="logit"), data=df2)
summary(model2)
Call:
glm(formula = InMichelin ~ Food, family = binomial(link = "logit"), 
    data = df2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3484  -0.8555  -0.4329   0.9028   1.9847  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -10.84154    1.86234  -5.821 5.83e-09 ***
Food          0.50124    0.08767   5.717 1.08e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 225.79  on 163  degrees of freedom
Residual deviance: 175.73  on 162  degrees of freedom
AIC: 179.73

Number of Fisher Scoring iterations: 4
In:
plot(jitter(InMichelin, 1, 0.03) ~ jitter(Food, 1, 0.1), data=df2, 
     xlim=c(10, 30), xlab="Food", ylab="InMichelin")
grid()
Xnew <- data.frame(Food=seq(10, 30, length=100))
ynew <- predict(model2, Xnew, type="response")
lines(Xnew$Food, ynew, lwd=2)

c <- coef(model2)
threshold <- -c[[1]]/c[[2]]

abline(h=0.5, col=rgb(0.8, 0.2, 0.2, 0.5), lwd=3)
abline(v=threshold, col=rgb(0.2, 0.7, 0.2, 0.5), lwd=3)
text(threshold + 2.4, 0.3, sprintf("Food = %4.2f", threshold))

질문/덧글

질문이 있습니다. 박사님. 오타가 나지 않은가 싶습니다. sh.h*** 2018년 9월 4일 1:15 오후

로지스틱 모형의 모수 추정의 Log Likelihood 함수에 관한 질문입니다.

LL = log∏θ*(χ)y*{1-θ(χ)}1-y 가 아닌

LL = log∏{θ*χ}y*{1-θ(χ)}1-y 이 옳은 것이 아닌가 싶어 질문을 남깁니다.

그림으로 파일을 올릴 수 없어
y승과 (1-y)승을 윗첨자로 적지 못한 점 양해바랍니다.

항상 좋은 자료 감사하게 보고있습니다.
박사님.

답변: 질문이 있습니다. 박사님. 오타가 나지 않은가 싶습니다. 관리자 2018년 9월 4일 3:17 오후

다음 식을 말씀하시는 것인가요?
$$
\begin{eqnarray}
\text{LL}
&=& \log \prod_{i=1}^N \theta_i(x_i)^{y_i} (1-\theta_i(x_i))^{1-y_i} \\
&=& \sum_{i=1}^N \left( y_i \log\theta_i(x_i) + (1-y_i)\log(1-\theta_i(x_i)) \right) \\
\end{eqnarray}
$$
덧글에도 `$$ 수식 $$` 의 형태로 LaTeX 명령을 쓸 수 있으니 수식을 사용해주시면 감사하겠습니다.

sh.h*** 2018년 9월 4일 6:42 오후

네 그렇습니다. 첫번째 줄 수식과 관련된 질문입니다.
파이 안의 수식이
$$ \theta _{i} (x _{i} ) ^{y _{i}}(1- \theta _{i} (x _{i} )) ^{1-y _{i}} $$이 아닌
$$( \theta _{i} x _{i} ) ^{y _{i}} (1- \theta _{i} (x _{i} )) ^{1-y _{i}} $$이렇게 되어야

두번째 줄의 수식[특히 시그마 안의 수식의 앞 부분]
$$ log( \theta _{i} x _{i} ) ^{y _{i}} $$ 에서
$$ y _{i} $$가 로그 앞으로 나와서
$$ y _{i} log \theta_{i} (x _{i} ) $$ 이렇게 되는것이 아닌가요?

답변: Latex 사용이 익숙치 않아서, 가독성이 좋지 않은 점 양해부탁드립니다. sh.h*** 2018년 9월 4일 7:25 오후

데이터 사이언스를 조금씩 공부중인데, 수학과 관련된 기초가 약해서 제가 이해를 잘못하고 있는 부분이 있을 수도 있을 것 같습니다. 답변 남겨주시면 감사하겠습니다. 박사님

답변: 답변: Latex 사용이 익숙치 않아서, 가독성이 좋지 않은 점 양해부탁드립니다. 관리자 2018년 9월 4일 7:33 오후

$\theta$가 $x$의 함수이기 때문에 $\theta(x)$라고 써야 합니다. 첨자 $i$는 $i$번째 데이터를 뜻합니다.

답변: 답변: 답변: Latex 사용이 익숙치 않아서, 가독성이 좋지 않은 점 양해부탁드립니다. sh.h*** 2018년 9월 5일 5:08 오후

아하!! 감사합니다. 박사님!!!