Логистическая регрессия пытается решить класс проблем, которые кажутся более простыми, чем линейная регрессия. В последнем случае мы хотим найти линию (или плоскость, или гиперплоскость), которая наилучшим образом предсказывает среднее числовое значение для некоторых данных на основе наших наблюдений. В логистической регрессии мы просто пытаемся предсказать, в какую категорию попадает наблюдение, на основе зависимых переменных. Мужчина или женщина? Больной или нет?

В этом сообщении блога я показываю математику, которая лежит в основе логистической регрессии, и показываю, как вы можете реализовать ее на языке программирования R. Доступ к файлам можно получить в моем репозитории Github здесь.

Суть логистической регрессии

Когда мы используем логистическую регрессию, мы пытаемся определить вероятность того, что наблюдение будет в определенном классе. Я начну с случая двух классов (K = 2). Давайте определим наши переменные для классов A и B.

«Сырая» модель, с которой мы начинаем, представлена ​​ниже. Мы отвечаем на вопрос: каковы шансы того, что данные наблюдения i попадут в категорию A по сравнению с B с учетом набора параметров β? ?

Небольшое примечание: если мы просто попытаемся предсказать отношение шансов, мы попытаемся предсказать значение функции, которая сходится к бесконечности при p → ∞. Эту проблему можно обойти, вычислив натуральный логарифм отношения шансов. Использование этой стратегии не меняет положение максимума, но создает набор значений, которые ограничены между 0 и 1.

Тогда модель, для которой мы пытаемся найти параметры β, выглядит следующим образом:

На следующем шаге для решения p используется довольно простая алгебра:

Поиск параметров

В отличие от линейной регрессии, где параметры модели могут быть подобраны путем инвертирования матрицы, логистическая регрессия должна рассматриваться как проблема оптимизации.

Функция максимального правдоподобия

Функция правдоподобия - это некоторая мера, которую мы используем для определения вероятности наблюдаемых нами результатов с учетом имеющихся у нас данных. Принцип оценки максимального правдоподобия гласит, что мы должны определить некоторую функцию, которая измеряет эту вероятность, а затем искать параметры, которые ее оптимизируют.

Для нашей задачи бинарной классификации мы будем предполагать, что наша функция правдоподобия управляется биномиальной функцией распределения вероятностей:

Поскольку нас интересует только максимизация этой функции, нам не нужно беспокоиться о количестве возможных вариантов. Если это не учитывать, наша функция максимального правдоподобия будет выглядеть следующим образом:

Мы еще далеки от завершения! Нам нужно взять производную вышеупомянутой функции, что окажется трудным с учетом мультипликативного качества приведенного выше выражения. Чтобы упростить задачу, мы можем просто взять натуральный логарифм и получить выражение, которое будет намного легче вычислить:

Метод Ньютона-Рафсона

Выше я упоминал, что нам нужно подойти к этому как к проблеме численной оптимизации. Очень часто используется метод Ньютона-Рафсона.

В контексте нашей проблемы логистической регрессии мы ищем критическую точку (которая содержит глобальный максимум) для функции правдоподобия. Значит, мы должны найти точку, в которой градиент равен нулю. В этом контексте мы должны использовать градиент и гессиан функции f (x).

Моя реализация метода Ньютона-Рафсона в R:

# Simple implementation of Newton-Raphson:
# Inputs:
# [1] func - The function to be optimized by the algorithm.
# [2] initial_values - An initial guess to be passed to "func".
# [3] tol - The stopping criteria to the model.  For this implementation, the tolerance is compared to the 2-norm
#           of the gradient:
# Outputs:
# [1]  params - the output of the logistic regression.
newton <- function(func, initial_values,tol = 1e-16){
  params <- initial_values
  check <-1
    while(check > tol){
      func_eval <- func(params)
      params <- params -  solve(func_eval$ddf)%*%func_eval$df
      check <- sqrt(t(func_eval$df)%*%func_eval$df)
    }
  return(params)
}

Градиент

Чтобы выразить это в матричной записи, пусть Y будет вектором из N элементов, представляющим результат класса для данных, и пусть P будет вектором длины N, составленным из элементов п_и. Тогда формула градиента для функции логарифмического правдоподобия:

Гессен

Следуя соглашению из [2], пусть W будет диагональной матрицей N на N с ненулевыми элементами W_ {i, i} = p_i (1-p_i). Имея это в виду, определите матрицу Гессе как:

Моя реализация R функции логарифма правдоподобия показана ниже:

# Implementation of binary logistic regression:
# Inputs:
# [1]  y - a column vector of N elements containing the actual classification results.
# [2]  X - a (p+1) by N design matrix.
# [3]  beta - a p-length column vector of coefficients.
# Outputs:
# [1]  f - The log-likelihood function evaluated using the inputs.
# [2]  df - p-length column vector representing the log-likelihood's gradient.
# [3]  ddf - a p by p Hessian matrix for the log-likelihood function.
L_L <- function(y,beta,X){
  # Calculate the vector of probabilities:
  pi <- exp(X%*%beta)/(1+exp(X%*%beta))
  f <- prod(ifelse(y==1,pi,1-pi))
  # First derivative:
  df <- t(t(y- pi)%*%X )
  # Second derivative:
  W <- diag(as.vector(pi*(1-pi)))
  ddf <- -t(X)%*%W%*%X
  # Output:
  output_list <- list(f = f, df = df,ddf = ddf)
  return(output_list)
}

Санитарная проверка

Чтобы убедиться, что моя реализация работает, я тестирую ее на нескольких разных примерах. Я решил опубликовать тот, что ниже, потому что он работает с реальным набором данных разумного размера.

Я сравнил свою реализацию со знаменитой функцией glm. Ожидание слишком многого с точки зрения минимизации градиента заставило меня столкнуться с некоторыми проблемами, связанными с сингулярностью матрицы Гессе, поэтому я сохранил допуск на уровне 1e-13, который был достаточно хорош для получения результатов, очень похожих на функцию glm:

Фактические данные

# Peter Caya
# 2017-03-19
# Simple script demonstrating binary classification using the standard glm function,
# my simpler L_L function for binary classification, and my multi-class 
# implementation.
rm(list = ls())
source("/home/peter/Dropbox/Data Science/ML Scripts/Log_Functions.R")
set.seed(13)
mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
admission <- as.matrix(mydata$admit)
beta <- as.matrix(rep(0,3))
vars <- as.matrix(mydata[,2:dim(mydata)[2]])
my_coefs <- newton(function(input){L_L(y = admission,X = vars,beta = input)},initial_values = beta,tol = 1e-13)
fromGLM <- glm(formula = admission~vars-1,family = binomial())
GLMcoefficients <- fromGLM$coefficients
my_results <- L_L(y,my_coefs,x)
glm_results <- L_L(y,GLMcoefficients,x)
L_L_Difference <- my_results$f - glm_results$f
DL_L_Difference <- t(my_results$df - glm_results$df)%*%(my_results$df - glm_results$df)
print(data.frame(matrix(GLMcoefficients),my_coefs))

Заключение

Мы можем видеть, что, хотя математическая основа для логистической регрессии немного сложнее, чем линейная регрессия, как только мы поймем ее, этот метод очень прост в реализации как проблема численной оптимизации.

Однако что, если мы хотим классифицировать объекты двух или более разных классов? Основные идеи остаются прежними, но готовый продукт немного сложнее. Первоначально я собирался добавить в эту статью свою реализацию логистической регрессии для нескольких переменных, но мне показалось, что она становится слишком длинной, так что это будет задокументировано в следующем сообщении в блоге.

Источники

[1] Элементы статистического обучения,

[2] Оценка максимального правдоподобия моделей логистической регрессии: теория и реализация

[3] Логистическая регрессия; Доху, И. Мартин, В. Стрюнн, Х.

Особая благодарность моему другу Тайлеру Уилкоксу, который прислал мне несколько источников, которые действительно помогли сделать этот пост более полным.