Логистическая регрессия пытается решить класс проблем, которые кажутся более простыми, чем линейная регрессия. В последнем случае мы хотим найти линию (или плоскость, или гиперплоскость), которая наилучшим образом предсказывает среднее числовое значение для некоторых данных на основе наших наблюдений. В логистической регрессии мы просто пытаемся предсказать, в какую категорию попадает наблюдение, на основе зависимых переменных. Мужчина или женщина? Больной или нет?
В этом сообщении блога я показываю математику, которая лежит в основе логистической регрессии, и показываю, как вы можете реализовать ее на языке программирования 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] Логистическая регрессия; Доху, И. Мартин, В. Стрюнн, Х.
Особая благодарность моему другу Тайлеру Уилкоксу, который прислал мне несколько источников, которые действительно помогли сделать этот пост более полным.