Цель этой статьи — дать читателю четкое представление о методе логистической регрессии. Он предполагает знакомство с основами теории вероятностей, линейной алгебры и дифференциального исчисления. Его первая часть исследует обоснование логистической регрессии наряду с необходимыми математическими деталями. Вторая часть посвящена реализации, для которой я использую Python.
Зачем использовать логистическую регрессию?
Это соблазняет прибегнуть к старой знакомой линейной регрессии, даже если целевая переменная является дихотомической (также известной как двоичная), однако она окажется недействительной с точки зрения моделирования вероятностей и, следовательно, адекватных классификаций. Причина в том, что такая модель имеет неограниченный диапазон, то есть простирается от -inf до +inf. Было изобретено решение, включающее функцию, которая отлично подходит для моделирования вероятностей, а именно сигмовидную функцию:
Вы также можете иногда сталкиваться с его другой формой:
Таким образом, значение функции приближается к 1, когда x переходит в +inf, и к 0, когда x переходит в -inf. Логистическая регрессия использует это поведение, делая вероятностные предположения.
Теория
Как и в линейной регрессии, вход логистической регрессии представляет собой линейную комбинацию признаков и весов (коэффициентов), связанных с ними:
который мы собираемся переписать как:
Внешний вид w-ноль и 1 может показаться запутанным, однако это то же самое y-пересечение, которое представляет бета-ноль. В контексте весов w-naught часто называют смещением.
С входом, подключенным к сигмовидной функции, мы имеем вероятностную форму логистической регрессии:
Однако также часто можно увидеть логарифм шансов, который используется в так называемой логит-форме логистической регрессии:
но почему мы можем захотеть использовать эту странную меру? Причина в необходимости выразить линейную комбинацию признаков и весов.
Решение выразить линейную комбинацию через логарифм шансов имеет двоякую основу: во-первых, оно возвращает знакомую «прямолинейную» запись, а, во-вторых, коэффициенты становятся несколько более интерпретируемыми.
Таким образом, вопрос в том, как мы расположим логит-прямую, какой будет кривая вероятности. Возникает вопрос: как выбрать веса, которые будут давать наиболее подходящую кривую?
Максимальная вероятность
Метод OLS (Обычные наименьшие квадраты), с помощью которого мы оптимизируем параметры для линейной регрессии, не будет работать с логит-функцией, проблема в том, что независимо от линии расстояния до точек данных всегда будут бесконечными. . Вместо этого мы вводим понятие максимальной вероятности, где вероятность — это вероятность наблюдения наших данных при условии, что модель верна:
Обратите внимание, что MLE (оценка максимального правдоподобия) требует независимости наблюдений, что позволяет умножать вероятности:
Однако из-за проблемы округления чисел с плавающей запятой мы рискуем потерять значение при умножении вероятностей, что почти гарантированно приведет к нулевой вероятности, независимо от выбранных нами весов. Именно поэтому при выполнении MLE функция правдоподобия подвергается логарифмическому преобразованию, которое переводит произведение вероятностей в сумму их логов.
Градиентный спуск
Градиентный спуск — это метод численной оптимизации, который итеративно приближается к минимуму функции стоимости J(w). Это схема алгоритма:
- Выбираются скорость обучения η (число обычно между 0,001 и 0,1) и количество итераций (рекомендуется взять не менее 1000);
- Генерируются начальные веса, которые обычно представляют собой случайные числа;
- Градиент (производная) вычисляется относительно каждого задействованного параметра w_i;
- Новые веса рассчитываются по следующей формуле:
Примечание: позаботьтесь о том, какую скорость обучения вы выберете. Если вы выберете относительно большое значение, ваши веса будут бесконечно отскакивать без сходимости. Наоборот, если η слишком мало, вы, возможно, состаритесь до того, как значение градиента достигнет порога допуска (произвольно выбранный параметр, близкий к нулю, на котором алгоритм останавливается).
Чтобы соблюсти определенные соглашения, мы определяем нашу функцию стоимости J как -LL, поскольку в противном случае это был бы градиентный подъем, который упоминается гораздо реже. Ниже приведен полный вывод градиента J:
Давайте разложим задачу и сначала получим градиенты для членов в сумме:
затем свяжите все вместе и уберите избыточность:
Итак, наша цель — вывести градиенты для всех весов и приблизить их к нулю. Хотя теоретически к задаче можно подойти аналитически с помощью системы уравнений, по мере увеличения количества признаков все быстро становится безумным. Представьте, что вы обучаете компьютер находить аналитическое решение с использованием 10 000 признаков, и, следовательно, система из 10 001 производной (вспомните точку пересечения) обнуляется. Видеть? Численные методы спасают положение.
Многоклассовая логистическая регрессия
Хотя по своей природе цель логистической регрессии состоит в том, чтобы различать только два класса, ее можно использовать для мультиклассовой (n > 2) классификации. В этой статье рассматривается метод One-vs-Rest, который подразумевает подгонку n бинарных логистических моделей к n наборам данных, предварительно обработанным по-разному. Концепцию лучше всего понять на примере:
Учтите, что мы должны различать Setosa, Versicolor и Virginica из набора данных Iris, который должен быть хорошо знаком тем, кто знаком с машинным обучением. Мы должны построить 3 разные логистические модели:
- Setosa (1) против Versicolor и Virginica (0);
- Versicolor (1) против Setosa и Virginica (0);
- Вирджиния (1) против Сетозы и Версиколора (0),
а затем мы суммируем результаты, при этом окончательным прогнозом является класс, который имеет наибольшую вероятность наблюдения 1. Скажем, для наблюдения у нас есть вероятность получить Setosa, равную 0,035, Versicolor — 0,23, а Virginica — 0,035. 0006. Таким образом, Versicolor будет предсказанным классом.
Логистическая регрессия в Python
Поздравляю с освоением теории и переходом ко второй части статьи. Здесь мы собираемся построить логистическую регрессию в Python. Мы сделаем это в объектно-ориентированной манере, поэтому, если вы чувствуете, что плохо разбираетесь в ООП, считайте, что вы втиснулись в него. Для получения базовых знаний советую посетить миникурс Кори Шафера.
Для этой реализации вам также необходимо установить numpy.
Мы начнем с объявления класса и указания параметров, установленных при инициализации:
class LogisticRegression:
def __init__(self, eta=.01, n_iter=100000, tolerance=1e-5, random_state=42):
self._eta = eta
self._n_iter = n_iter
self._tolerance = tolerance
self.__random_state = random_state
Пара замечаний по параметрам:
- эта = η = скорость обучения;
- n_iter – количество итераций до остановки градиентного спуска;
- допуск – это порог, которого должны достичь все градиенты, чтобы алгоритм остановился;
- random_state — это случайное начальное число, на котором генерируются начальные веса. Указано для воспроизводимости.
Мы добавляем статический метод _sigmoidutility, который принимает матрицу признаков и массив весов и создает массив вероятностей:
@staticmethod
def _sigmoid(X, w):
denominator = 1 + np.exp(-(X * w).sum(axis=1))
return 1 / denominator
Теперь мы полностью сосредоточимся на методе fit, который я собираюсь объяснить шаг за шагом с помощью встроенных комментариев. Я включил параметр normalize, который дополнительно масштабирует функции, что делает конвергенцию градиентного спуска более плавной и быстрой:
def fit(self, X, y, normalize=False):
# Convert the input into numpy arrays
# Make y a 1-d array in case a matrix is supplied
self._X = np.array(X, dtype=np.dtype('float64'))
self._y = np.array(y, dtype=np.dtype('int64')).squeeze()
# We save normalize as the instance attribute;
# if set to True, test data in the predict method
# will also be normalized.
# We also save the data range as well as the minimum value.
self.__normalize = normalize
if self.__normalize:
self.__Xmin = self._X.min(axis=0)
self.__Xrange = self._X.max(axis=0) - self.__Xmin
self._X = (self._X - self.__Xmin) / self.__Xrange
# Check if the problem is multiclass:
self._classes = np.unique(self._y)
if len(self._classes) > 2:
# If we have more than 2 classes,
# we set the corresponding boolean to True
# and prepare a container for n binary models.
self.__multiclass = True
self.models_ = []
for class_ in self._classes:
# Setting 1 where an observation belongs to the
# class and 0 where it does not.
y = np.zeros(shape=self._y.shape)
y[self._y == class_] = 1
# Initialize and fit the model.
lr = LogisticRegression(
eta=self._eta,
n_iter=self._n_iter,
tolerance=self._tolerance,
random_state=self.__random_state
)
lr.fit(self._X, y)
# We initialize normalize as False regardless
# of whether or not the main model has True
# because, if it does, self._X is already normalized
# Instead, we set the necessary attributes after
# the model is fit.
if self.__normalize:
lr.__normalize = self.__normalize
lr.__Xmin = self.__Xmin
lr.__Xrange = self.__Xrange
self.models_.append(lr)
self.__fit = True
return self
else:
self.__multiclass = False
# We add the bias term to the data to fit the intercept.
self._X = np.concatenate(
[np.ones(shape=(len(X), 1)), self._X],
axis=1
)
# Generate the initial weights.
rs = np.random.RandomState(seed=self.__random_state)
self.w_ = rs.normal(size=(self._X.shape[1], ))
# Gradient descent
for _ in range(self._n_iter):
grad = ((self._sigmoid(self._X, self.w_) - self._y)[:, np.newaxis] * self._X).sum(axis=0)
self.w_ -= self._eta * grad
print(grad)
if all(np.absolute(grad) < self._tolerance):
break
self.intercept_ = self.w_[0]
self.coef_ = self.w_[1:]
self.__fit = True
return self
Давайте проверим, что у нас есть на данный момент:
from sklearn.model_selection import train_test_split from sklearn import datasets data = datasets.load_iris() X, y = data['data'][:, 3:4], data['target'] # Use one feature for simplicity (petal width). X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=42) lr = LogisticRegression(eta=.05) lr.fit(X_train, y_train, normalize=True) print(lr.models_[0].w_) print(lr.models_[1].w_) print(lr.models_[2].w_)
Теперь мы представляем метод predict, который принимает следующие параметры:
- X_тест. Проверка данных. Если модель была подобрана с использованием нормализованной версии обучающих данных, данные тестирования также будут нормализованы.
- вероятно. Возвращать ли вероятности. Они добавляются справа от массива предсказанных классов.
- порог. Правило принятия решения: если P ≥ порога, наблюдение классифицируется как 1 и 0 в противном случае. Следует манипулировать с осторожностью. Например, если жизнь пациента зависит от обнаружения болезни, порог следует снизить. Игнорируется, если модель мультиклассовая.
def predict(self, X_test, proba=False, threshold=.5):
if not self.__fit:
raise TypeError('The model has not been fit.')
if self.__multiclass:
# Extract the probabilities from each of the fitted
# models.
probas = map(lambda model: model.predict(X_test, proba=True)[:, [1]], self.models_)
# Tie them together.
self.proba_ = np.concatenate(list(probas), axis=1)
# Make the prediction: a class with the highest
# probability is chosen as the predicted class.
self.y_hat_ = np.array([self._classes[idx] for idx in self._classes[np.argmax(self.proba_, axis=1)]])
if proba:
return np.concatenate([self.y_hat_[:, np.newaxis], self.proba_], axis=1)
else:
return self.y_hat_
else:
self._X_test = np.array(X_test, dtype=np.dtype('float64'))
# Normalize the testing data.
if self.__normalize:
self._X_test = (self._X_test - self.__Xmin) / self.__Xrange
# Append the bias term.
self._X_test = np.concatenate([np.ones(shape=(self._X_test.shape[0], 1)), self._X_test], axis=1)
# Calculate the probabilities.
self.proba_ = self._sigmoid(self._X_test, self.w_)
self.y_hat_ = np.zeros(shape=(self.proba_.shape[0], ))
self.y_hat_[self.proba_ >= threshold] = 1
if proba:
return np.concatenate([self.y_hat_[:, np.newaxis], self.proba_[:, np.newaxis]], axis=1)
else:
return self.y_hat_
Давайте посмотрим на метод в действии:
print(lr.predict(X_test, proba=True))
Вывод:
[[0.00000 1.00000 0.26618 0.00000]
[2.00000 0.00000 0.35582 0.42736]
[1.00000 0.00000 0.33544 0.01199]
…
[0.00000 0.99901 0.28415 0.00000]
[2.00000 0.00000 0.40533 0.99991]
[1.00000 0.00000 0.34896 0.15900]]
Мы видим, как выглядит процесс принятия решений, но как будет работать такая модель?
from sklearn.metrics import accuracy_score print(accuracy_score(lr.predict(X_test), y_test))
Вывод:
0.8947368421052632
Что ж, он действительно неплохо справился со своей задачей, но можем ли мы сделать лучше? Я не собираюсь посвящать остальную часть этой статьи выбору функций, так как эта тема выходит за ее рамки, и поэтому мы просто подставляем все имеющиеся у нас данные.
X, y = data['data'], data['target'] X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=42) lr.fit(X_train, y_train, normalize=True) print(accuracy_score(lr.predict(X_test), y_test))
Вывод:
0.9473684210526315
Давайте сравним наш результат с реализацией scikit-learn:
from sklearn.linear_model import LogisticRegression as LinearRegression_ lr_sklearn = LogisticRegression_(penalty='none', tol=1e-5, random_state=42, max_iter=100000) # Same hyperparams lr_sklearn.fit(X_train, y_train) print(accuracy_score(lr_sklearn.predict(X_test), y_test))
Вывод:
0.8947368421052632
Не буду хвастаться этим незначительным преимуществом, так как оно, скорее всего, связано со случаем.
Кроме того, я полагаю, что можно сделать даже лучше, чем 94,7%, если мы удалим «шумные» функции.
Вывод
В этой статье я постарался быть одновременно кратким и всеобъемлющим. Я старался как можно яснее передать материал и сделать особый акцент на вещах, которые, на мой взгляд, были представлены либо поверхностно, либо чрезмерно запутанно в других источниках, которые мне приходилось перебирать во многих источниках. Я надеюсь, что вам понравилось это руководство, и если у вас остались какие-либо вопросы без ответа, свяжитесь со мной, и я отвечу на них в той мере, в какой это позволяет мой опыт.
