본문으로 바로가기



Optimization with R - 5. 뉴턴 방법을 사용한 로지스틱 회귀분석 모수 추정하기

category Statistics/Optimization 2017.07.16 22:49
Optimization with R - 5. 뉴턴 방법을 사용한 로지스틱 회귀분석 모수 추정하기

Multivariate Newton’s method with R

지난 포스팅에서 우리는 입력변수가 2개 이상인 함수에 대하여 뉴턴법을 적용하여 보았다. 오늘은 이에 관한 실제 R코드를 짜고 이를 가상의 데이터를 만들어 적용하여 보도록 하자. 좀 더 쉬운 코드를 위하여 우리가 사용하는 최적화 방법은 제약조건이 없는(unconstrained) 상황의 최적화를 가정한다.

우리가 정의할 함수의 이름은 newton_opt이라고 하고, 입력값으로는 grad, Hessian, x_now, 그리고 나머지 루프를 빠져나오기 위한 값들인 stop_criteria, max_iter, 마지막으로 stepsize를 설멍하였다.

newton_opt <- function(grad, Hessian, x_now,
                       stop_criteria = .Machine$double.eps,
                       max_iter = 10^4,
                       stepsize = 0.001)

뉴턴 방법 R로 구현하기

핵심 스텝 구현하기

지난 시간에 알아본 다변수 함수에 적용되는 뉴턴법의 스텝은 다음과 같았다.

\[ \mathbf {x} _{next}=\mathbf {x}_{now}-[\mathbf {H}(\mathbf {x} _{now})]^{-1}\nabla f(\mathbf {x} _{now}) \]

이 스텝을 구현함에 있어서, 헤시안(Hessian) 행렬의 역행렬을 구하는 대신 다음과 같은 두 가지 변형을 취했다.

# newton method
sterm <- solve(Hessian(x_now), grad(x_now))
x_next <- x_now - stepsize * sterm
  1. 역행렬의 계산 대신 선형방정식을 푸는 과정을 택한다. 이는 역행렬 계산보다 좀 더 안정적이다.

  2. stepsize를 사용하여 스텝이 너무 커서 발산하는 경우를 방지한다.

\[ \begin{align*} solve\Rightarrow & \mathbf{H}(\mathbf{x}_{now})y=\nabla f(\mathbf{x}_{now})\\ \Rightarrow & \mathbf{x}_{next}=\mathbf{x}_{now}-stepsize*y \end{align*} \]

루푸 탈출 방법 설정하기

정의할 newton_opt함수 안에 들어갈 반복문을 빠져나올 조건을 다음과 같이 설정하자.

  1. x_now와 x_next가 거의 같아서 움직이지 않을 경우
  2. gradient값이 0과 아주 가까울 경우
  3. 미리 설정한 반복 횟수를 초과했을 경우
while ((distance > stop_criteria) &
         (iters < max_iter) &
         (sqrt(sum(grad(x_now)^2)) > stop_criteria))

R 코드 전체

# Multi demesional newton's method in optimization
newton_opt <- function(grad, Hessian, x_now,
                       stop_criteria = .Machine$double.eps,
                       max_iter = 10^4,
                       stepsize = 0.01){
  # variable initialization
  iters <- 0
  x_next <- 0
  distance <- Inf
  
  # do while until
  # 1) x does not move
  # 2) exceed the maximun iteration
  # 3) the function value become close to zero
  while ((distance > stop_criteria) &
         (iters < max_iter) &
         #(min(x_now) > stop_criteria)
         (sqrt(sum(grad(x_now)^2)) > stop_criteria)){
    # newton method
    sterm <- solve(Hessian(x_now), grad(x_now))
    x_next <- x_now - stepsize * sterm
    
    distance <- sqrt(sum((x_now - x_next)^2))
    #distance <- min(x_now - x_next)
    iters <- iters + 1
    x_now <- x_next
  }
  
  # return the final point
  if (iters > max_iter){
    warning("\n", "Over the maximum number of iteration.", "\n")
  }
  x_now
}

Logistic regression 모수 추정 적용

위에서 정의한 newton_opt 함수를 사용하여 로지스틱 회귀분석의 모수 추정을 해보도록 하자.

데이터 셋 발생

먼저 가상의 데이터 셋을 발생시켰는데, 시험의 결과를 종속변수(y \(\in\) {0, 1}, 패스 1, 탈락 0)로 평균 공부시간(x1)과, 평균 수업 중 질문 횟수(x2)를 독립변수 데이터로 가지고 있다고 가정하자.

library(ggplot2)

set.seed(2017)
x1 <- c(rnorm(100, 5, 2), rnorm(100, 10, 2))
x2 <- c(rnorm(100, 3, 2), rnorm(100, 6, 2))
y <- rep(c(0, 1), c(100, 100))
plotdata <- data.frame(studytime = x1, nquestion = x2, pass = y)

p <- ggplot(data = plotdata) +
  geom_point(mapping =
               aes(x = studytime,
                   y = nquestion,
               color = factor(pass))) + theme_bw() +
  scale_colour_discrete(guide = guide_legend(reverse=TRUE),
                        name="Exam",
                        labels=c("Fail", "Pass"))
p

데이터를 보면 평균 공부시간이 길고, 수업시간에 질문을 많이 할수록 시험에 붙을 확률이 높아지는 것으로 보인다. 이러한 데이터를 바탕으로 로지스틱 모델을 설정해 보도록 하자.

로지스틱(logistic) 모델

로지스틱 모델은 위의 데이터에서 종속변수인 시험 결과가 베르누이(Bernoulli) 분포를 따른다고 가정한다.

\[ \begin{align*} y_{i} & \sim Bernoulli\left(\sigma\left(\beta_{0}+\beta_{1}x_{1i}+\beta_{2}x_{2i}\right)\right)\\ & \sim Bernoulli\left(\sigma\left(\mathbf{x}_{i}^{T}\beta\right)\right) \end{align*} \] 위의 데이터에서 i는 1에서 200까지의 값을 가질 수 있으며, \(\sigma\) 함수는 \(\sigma(x)=1/(1+exp(-x))\)로 정의하며, 학생들의 탈락 여부는 독립을 가정한다.

MLE 구하기

주어진 데이터를 바탕으로 likelihood 함수와 log likelihood 함수는 다음과 같이 구할 수 있다.

\[ \begin{align*} p\left(\underline{y}|\mathbf{X},\beta\right) & =\prod_{i=1}^{n}p\left(y_{i}|\mathbf{x}_{i},\beta\right)\\ & =\prod_{i=1}^{n}\sigma\left(\mathbf{x}_{i}^{T}\beta\right)^{y_{i}}\left(1-\sigma\left(\mathbf{x}_{i}^{T}\beta\right)\right)^{1-y_{i}}\\ & =\prod_{i=1}^{n}\pi_{i}^{y_{i}}\left(1-\pi_{i}\right)^{1-y_{i}},\,\mbox{ where }\pi_{i}:=\sigma\left(\mathbf{x}_{i}^{T}\beta\right) \end{align*} \]

MLE의 경우 이 likelihood 함수를 최대로 만드는 베타값을 구한다. 따라서 우리는 앞에서의 논조를 맞추기 위해 -log를 취한 값을 최소화 시키는 베타값을 찾도록 하겠다.

\[ \begin{align*} \ell\left(\beta\right) & =-log\left(p\left(\underline{y}|\mathbf{X},\beta\right)\right)\\ & =-\sum_{i=1}^{n}\left\{ y_{i}log\left(\pi_{i}\right)+\left(1-y_{i}\right)log\left(1-\pi_{i}\right)\right\} \\ & =-\sum_{i=1}^{n}\left\{ y_{i}log\left(\sigma\left(\mathbf{x}_{i}^{T}\beta\right)\right)+\left(1-y_{i}\right)log\left(1-\sigma\left(\mathbf{x}_{i}^{T}\beta\right)\right)\right\} \end{align*} \]

주어진 negative loglikehood 함수에 뉴턴 방법을 적용하기 위해서는 기울기와 헤시안 행렬을 구해야 한다.

기울기(Gradient)

\[ \begin{align*} \frac{\partial\ell\left(\beta\right)}{\partial\beta} & =-\sum_{i=1}^{n}\left\{ y_{i}\frac{\partial}{\partial\beta}\left(log\left(\sigma\left(\mathbf{x}_{i}^{T}\beta\right)\right)\right)+\left(1-y_{i}\right)\frac{\partial}{\partial\beta}\left(log\left(1-\sigma\left(\mathbf{x}_{i}^{T}\beta\right)\right)\right)\right\} \\ \tag{1} & =-\sum_{i=1}^{n}\left\{ y_{i}\frac{\sigma'\left(\mathbf{x}_{i}^{T}\beta\right)}{\sigma\left(\mathbf{x}_{i}^{T}\beta\right)}\mathbf{x}_{i}+\left(1-y_{i}\right)\frac{-\sigma'\left(\mathbf{x}_{i}^{T}\beta\right)}{\left(1-\sigma\left(\mathbf{x}_{i}^{T}\beta\right)\right)}\mathbf{x}_{i}\right\} \end{align*} \]

뭔가 엄청 복잡해 보인다. 하지만 \(\sigma\) 함수의 좋은 특성 때문에 위의 식을 간략하게 나타낼 수 있다. 먼저 \(\sigma\) 함수를 미분해보자.

\[ \begin{align*} \sigma\left(x\right) & =\frac{1}{1+exp\left(-x\right)}=\left(1+exp\left(-x\right)\right)^{-1}\\ \sigma'\left(x\right) & =\frac{exp\left(-x\right)}{\left(1+exp\left(-x\right)\right)^{2}} \end{align*} \]

따라서 식 (1)에 나타나는 시그마 함수의 분수꼴은 다음과 같이 간단하게 나타낼 수 있다.

\[ \begin{align*} \frac{\sigma'\left(x\right)}{\sigma\left(x\right)} & =\frac{\frac{exp\left(-x\right)}{\left(1+exp\left(-x\right)\right)^{2}}}{\frac{1}{1+exp\left(-x\right)}}=\frac{exp\left(-x\right)}{1+exp\left(-x\right)}\\ \frac{\sigma'\left(x\right)}{1-\sigma\left(x\right)} & =\frac{\frac{exp\left(-x\right)}{\left(1+exp\left(-x\right)\right)^{2}}}{\frac{exp\left(-x\right)}{1+exp\left(-x\right)}}=\frac{1}{1+exp\left(-x\right)} \end{align*} \]

따라서 식 (1)에 위의 결과를 집어넣어보면 다음과 같다.

\[ \begin{align*} \frac{\partial\ell\left(\beta\right)}{\partial\beta} & =-\sum_{i=1}^{n}\left\{ y_{i}\frac{\sigma'\left(\mathbf{x}_{i}^{T}\beta\right)}{\sigma\left(\mathbf{x}_{i}^{T}\beta\right)}\mathbf{x}_{i}+\left(1-y_{i}\right)\frac{-\sigma'\left(\mathbf{x}_{i}^{T}\beta\right)}{\left(1-\sigma\left(\mathbf{x}_{i}^{T}\beta\right)\right)}\mathbf{x}_{i}\right\} \\ & =-\sum_{i=1}^{n}\left\{ y_{i}\left(\frac{exp\left(-\mathbf{x}_{i}^{T}\beta\right)}{1+exp\left(-\mathbf{x}_{i}^{T}\beta\right)}\right)\mathbf{x}_{i}-\left(1-y_{i}\right)\left(\frac{1}{1+exp\left(-\mathbf{x}_{i}^{T}\beta\right)}\right)\mathbf{x}_{i}\right\} \\ & =-\sum_{i=1}^{n}\left\{ y_{i}\left(1-\pi_{i}\right)\mathbf{x}_{i}-\left(1-y_{i}\right)\pi_{i}\mathbf{x}_{i}\right\} \\ & =\sum_{i=1}^{n}\left\{ \left(\pi_{i}-y_{i}\pi_{i}-y_{i}+y_{i}\pi_{i}\right)\mathbf{x}_{i}\right\} \\ & =\sum_{i=1}^{n}\left\{ \left(\pi_{i}-y_{i}\right)\mathbf{x}_{i}\right\} \\ & =\left(\pi_{1}-y_{1}\right)\left(\begin{array}{c} 1\\ x_{11}\\ x_{21} \end{array}\right)+...+\left(\pi_{n}-y_{n}\right)\left(\begin{array}{c} 1\\ x_{1n}\\ x_{2n} \end{array}\right)\\ & =\left(\begin{array}{ccc} 1 & ... & 1\\ x_{11} & ... & x_{1n}\\ x_{21} & ... & x_{2n} \end{array}\right)\left(\begin{array}{c} \pi_{1}-y_{1}\\ \vdots\\ \pi_{n}-y_{n} \end{array}\right)=X^{T}\left(\underline{\pi}-\underline{y}\right) \end{align*} \]

헤시안(Hessian) 행렬

다음은 헤시안 행렬을 구할 차례이다. 헤시안 행렬은 위에서 구한 기울기 벡터를 베타 트랜스 포즈로 다시 미분을 해줘야 한다.

\[ \begin{align*} \frac{\partial^{2}\ell\left(\beta\right)}{\partial\beta^{T}\beta} & =\frac{\partial}{\partial\beta^{T}}\left(X^{T}\underline{\pi}-X^{T}\underline{y}\right)\\ & =\frac{\partial}{\partial\beta^{T}}\left(X^{T}\underline{\pi}\right)\\ & =\sum_{i=1}^{n}\frac{\partial}{\partial\beta^{T}}\left(\pi_{i}\mathbf{x_{i}}\right)\\ & =\sum_{i=1}^{n}\frac{\partial}{\partial\beta^{T}}\left(\sigma\left(\mathbf{x}_{i}^{T}\beta\right)\mathbf{x_{i}}\right)\\ & =\sum_{i=1}^{n}\left(\sigma'\left(\mathbf{x}_{i}^{T}\beta\right)\mathbf{x_{i}}\mathbf{x_{i}}^{T}\right)\\ & =X^{T}DX \end{align*} \]

여기서 D 행렬은 다음과 같은 대각행렬이 되며, 이것을 앞에서 정의한 \(\pi\)를 사용하여 좀 더 간단히 나타낼 수 있다.

\[ \begin{align*} D: & =Diag\left(\sigma'\left(\mathbf{x}_{1}^{T}\beta\right),...,\sigma'\left(\mathbf{x}_{n}^{T}\beta\right)\right)\\ & =Diag\left(\pi_{1}\left(1-\pi_{1}\right),...,\pi_{n}\left(1-\pi_{n}\right)\right) \end{align*} \]

이것으로 바탕으로 R함수를 작성해 보도록 하자.

# negative loglikelihood
nloglik <- function(beta, X, y){
  -sum(y * X %*% beta - log(1 + exp(X %*% beta)))
}

# gradient function
dnloglik <- function(beta, X, y) {
  inner <- X %*% beta
  pi <- exp(inner) / (1 + exp(inner))
  t(t(pi - y) %*% X)
}

Hnloglik <- function(beta, X, y) {
  inner <- X %*% beta
  pi <- exp(inner) / (1 + exp(inner))
  D <- diag(as.vector(pi * (1 - pi)))
  # sweep(X, 1, pi * (1-pi), "*")
  # = D %*% X
  t(X) %*% sweep(X, 1, pi * (1 - pi), "*")
}

자, 이것으로 길고 길었던 준비과정이 끝났다. 우리가 구현한 함수 newton_opt 함수에 위에서 정의한 기울기 함수와 헤시안 행렬 함수를 집어넣고 결과를 확인해보자.

X <- cbind(Intercept = 1, x1, x2)
gradf <- function(x){dnloglik(x, X, y)}
Hessianf <- function(x){Hnloglik(x, X, y)}
x_init <- c(0,0,0)
result <- newton_opt(gradf, Hessianf, x_init)
result
##                 [,1]
## Intercept -16.034103
## x1          1.320595
## x2          1.252586

우리가 짠 코드가 맞는지 stat 패키지에 들어있는 glm 함수를 사용하여 답을 확인해보자.

result2 <- glm(y~x1+x2, family="binomial")
result2$coefficients
## (Intercept)          x1          x2 
##  -16.034103    1.320595    1.252586

위의 결과에서 우리가 짠 코드가 정상적으로 작동하였으며, 정상적으로 작동하고 있다는 것을 확인 할 수 있다. 몇몇 독자들은 정의된 stepsize를 변경하고, 시작점을 변경 시켜가면서 시도해볼 수 있겠다. 시작점과 stepsize에 따라서 뉴턴법은 수렴하기도 하고, 그렇지 않기도 하다. 즉, 올바른 시작점과 stepsize는 알고리즘이 수렴을 하는데에 영향을 미친다는 것을 깨달을 수 있을 것이다. 이제 결과가 나왔으니 이것을 기준으로 판단하는 선을 그을 수 있다. 피팅이 잘 되었는지 그림으로 확인해 보도록 하자.

addline <- function(a, b, c){
  my_slope <- -a / b
  my_intercept <- -c / b
  geom_abline(slope = my_slope,
              intercept = my_intercept)
}
p + addline(result[2], result[3], result[1])

주어진 데이터에 잘 적합되었음을 알 수 있다. 마지막으로 새로운 데이터 셋에는 어느정도 잘 들어맞는지 확인하면서 포스팅을 마치도록 하겠다.

set.seed(2020)
x1 <- c(rnorm(100, 5, 2), rnorm(100, 10, 2))
x2 <- c(rnorm(100, 3, 2), rnorm(100, 6, 2))
y <- rep(c(0, 1), c(100, 100))
plotdata <- data.frame(studytime = x1, nquestion = x2, pass = y)

p <- ggplot(data = plotdata) +
  geom_point(mapping = aes(x = studytime,
                           y = nquestion,
                           color = factor(pass))) +
  scale_colour_discrete(guide = guide_legend(reverse=TRUE),
                        name="Exam",
                        labels=c("Fail", "Pass")) + theme_bw()

X <- cbind(Intercept = 1, x1, x2)
t_result <- newton_opt(gradf, Hessianf, x_init)
p + addline(t_result[2], t_result[3], t_result[1])


Reference

[1] Wikipedia Newton’s method page: https://en.wikipedia.org/wiki/Newton%27s_method

[2] Wikipedia Newton’s method in optimization page: https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization

[3] Monahan, John F. A primer on linear models. CRC Press, 2008.

[4] U of Iowa STAT7400 - Luke Tierney’s lecture note: http://homepage.divms.uiowa.edu/~luke/classes/STAT7400/notes.pdf

[5] Petersen, Kaare Brandt, and Michael Syskind Pedersen. “The matrix cookbook.” Technical University of Denmark 7 (2008): 15.

[6] Cookbook for R(Plot Legends): http://www.cookbook-r.com/Graphs/Legends_(ggplot2)/


SHARE TO



티스토리 툴바