С подробной реализацией на Python

Обобщенные линейные модели (GLM) – это регрессионные модели, в которых мы обобщаем линейное предположение обычной модели линейной регрессии. Из-за этой нелинейности оценка параметра регрессии будет не такой простой, как оценка параметра линейной регрессии.

Алгоритм Повторяющиеся взвешенные наименьшие квадраты (IRLS) или иногда также Повторяющиеся взвешенные наименьшие квадраты (IWLS) — это метод нахождения оценок максимального правдоподобия обобщенных линейных моделей. Это расширение метода взвешенных наименьших квадратов. Начнем с краткого предисловия.

Фон

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

Этот метод будет иметь ошибки со средними нулевыми и постоянной дисперсией.

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

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

Мы все еще не можем использовать это для GLM по определенной причине. Y-переменная GLM отличается от прогнозируемой переменной. Рассмотрим данные ниже (которые мы также будем использовать для кода).

Здесь мы видим, что переменная y — это переменная count, указывающая количество успехов. Мы сохраняем максимальное количество успехов равным одному для нашего случая, что делает это проблемой логистической регрессии. Биномиальная (и логистическая) регрессия предсказывает не количество успехов экземпляра, а вероятность успеха. Эта разница между переменной y и желаемой прогнозируемой переменной проблематична, потому что мы не можем поместить нашу переменную y во взвешенное уравнение наименьших квадратов. Чтобы решить эту проблему, мы вводим алгоритм IRLS. Но прежде чем мы перейдем к алгоритму, мне нужно объяснить основы GLM.

Обобщенные линейные модели (GLM)

Как следует из названия, GLM представляют собой обобщение линейной регрессии, где прогнозируемая переменная связана с линейной моделью через функцию связи, обозначенную буквой g.

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

Здесь μ — это вероятность предсказания. Но вы можете использовать это для моделирования любой другой переменной с любым другим ограничением предметной области с соответствующим выбором функции связи.

Примечание. Функция связи должна быть дифференцируемой.

Давайте поговорим о функции ссылки. В нашем случае нам нужно найти функцию, которая преобразует бесконечный домен в домен [0, 1]. Многие функции обладают этим свойством. Я покажу вам самую популярную функцию связи для этого преобразования.

Вышеупомянутая функция широко известна как логит-функция. Функция логит-связи работает с двоичными переменными y для прогнозирования вероятностей.

Примечание. GLM действительны только в том случае, если переменная y относится к экспоненциальному семейному распределению.

Экспоненциальная семья

Логистическая регрессия имеет бинарную переменную y со значениями либо 1, либо 0. Можно сказать, что эти данные получены из распределения Бернулли с непостоянным параметром вероятности. Именно эту вероятность мы и пытаемся смоделировать. Поскольку распределение Бернулли или биномиальное распределение относится к экспоненциальному семейству, мы можем смоделировать его с помощью GLM. Но что это за экспоненциальная семья?

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

Я покажу, как выразить биномиальное распределение в этой форме.

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

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

Среднее значение биномиального распределения равно np, а дисперсия - np(1-p). Поскольку мы знаем этоa(ϕ) = 1 из вышеприведенного вывода, я могу представить свою функцию дисперсии следующим образом.

Если мы заменим μ на np, мы получим дисперсию np(1-p). Предположим, что нет прямой связи между средним значением и дисперсией, как в примере с нормальным распределением. В этом случае вам придется найти функцию дисперсии вручную, используя формулу.

Повторяющиеся повторно взвешенные наименьшие квадраты

Мы начнем с определения «новой» проблемы регрессии,

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

Обратите внимание, что для GLM член a(ϕ) можно полностью игнорировать, потому что он аннулируется в следующем расчете.

Отсюда мы можем оценить параметры, используя взвешенное решение наименьших квадратов. Здесь a(ϕ) отменяется.

Используя текущую оценку параметра, находим новое значение μ,

Напомним, что мы определили z и Σявляются функциями μ, μ — функцией β, а β — функцией z и Σ. Мы оцениваем параметры, повторяя это рекурсивное отношение.

Алгоритм останавливается, если увеличение логарифмической вероятности больше не является значительным.

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

Давайте скомпилируем. Алгоритм повторяющихся перевзвешенных наименьших квадратов:

  1. Инициализируйте μ в пределах определенного домена. Я инициализирую массивом 0.5вероятностей.
  2. Учитывая текущее значение μ, рассчитайте z и Σ, используя уравнение 1 и eуравнение 2.
  3. Учитывая текущее значение z и Σ, вычислить β по формуле взвешенных наименьших квадратов; eвопрос 3.
  4. Учитывая текущее значение β, рассчитайте μ, используя уравнение 4.
  5. Если увеличение логарифмического правдоподобия (уравнение 5) меньше заданного эпсилон, остановитесь. В противном случае вернитесь к шагу 2.

Код Python

Давайте соберем все кусочки головоломки вместе. Во-первых, функция ссылки.

Далее наши кусочки алгоритма,

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

Обратите внимание, что для X в параметре требуется столбец перехвата. Запустив следующую строку кода,

b = IWLS(X, y)

Возвращает оценки параметров в переменную b.

IWLS completed on 8 iterations
array([1.10999575, 9.12479957, 2.18745957])

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

phat = ig(X @ b)
yhat = np.where(phat >= 0.5, 1, 0)

Здесь phat хранит прогнозы вероятности, а yhat хранит биномиальный прогноз.

Сравнение

Насколько мне известно, в библиотеке Python scikit-learn нет решателя IRLS, поэтому мы не можем сравнивать наши результаты с их результатами. Но, к счастью, пакет R faraway использует IRLS в качестве решателя GLM. Давайте сравним.

Бег

summary(model)

Вернулся бы,

Call:
glm(formula = Success ~ ., family = binomial, data = df)
Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.72517  -0.12458  -0.00007   0.48620   1.07683
Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)    1.110      1.283   0.865    0.387
Height         9.125      8.489   1.075    0.282
Weight         2.187      2.099   1.042    0.297
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 13.8629  on 9  degrees of freedom
Residual deviance:  5.6807  on 7  degrees of freedom
AIC: 11.681
Number of Fisher Scoring iterations: 8

Вспомним наш предыдущий вывод,

IWLS completed on 8 iterations
array([1.10999575, 9.12479957, 2.18745957])

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

Поздравления

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