UPDATE: 2022-09-04 10:28:51

はじめに

このノートではロジスティク回帰モデルのパラメタを最尤法で計算する方法についてまとめておく。この記事は以前ブログに書いたものを転記しただけで追記も何もしていない。

ロジスティク回帰モデル

数理部分は、下記を参照しました。パラメタを推定するところは、optim()を使いました。

library(tidyverse)
set.seed(1989)
df <-  tibble(y = sample(c(0, 1), 10, replace = TRUE),
              x1 = rnorm(10),
              x2 = rnorm(10))

X <- df %>% select(starts_with("x"))
y <- df %>% select(y)


# Sigmoid function
sigmoid <- function(x) {
  1.0 / (1.0 + exp(-x))
}

# Cost function
cost_func <- function(theta, X, y){
  h <- sigmoid(X %*% theta)
  cost <- -1*(t(y) %*% log(h) + t(1-y) %*% log(1-h)) 
  cost
}

# Gradient function
gradient <- function(theta, X, y){
  h <- sigmoid(X %*% theta)
  grad <- (t(X) %*% (h - y)) 
  grad
}

log_reg <- function(X, y){
  # Add intercept and convert matrix
  X <- X %>% mutate(alpha = 1) %>% select(alpha, everything()) %>% as.matrix()
  y <- y %>% as.matrix()
  # Initialize parms
  theta <- matrix(rep(0, ncol(X)), nrow = ncol(X))
  # Optimize by BFGS(準ニュートン法)
  res_opt <- optim(par = theta, fn = cost_func, gr = gradient, X = X, y = y, method = "BFGS")
  # Get parms
  return(res_opt$par)
}

実行結果はこちら。

log_reg(X = X, y = y)
##            [,1]
## [1,] -1.1322717
## [2,] -0.4637470
## [3,] -0.7010454

glm()の結果と一緒みたいなのでOKでしょう。

g <- glm(y ~ x1 + x2, data = df, family = binomial(link = "logit"))
g$coefficients %>% as.matrix()
##                   [,1]
## (Intercept) -1.1322713
## x1          -0.4637478
## x2          -0.7010468