Введение в анализ данных

Домашнее задание 3, сложная часть. Линейная и логистическая регрессии.

Правила, прочитайте внимательно:

• Выполненную работу нужно отправить телеграм-боту @thetahat_ds25_bot. Для начала работы с ботом каждый раз отправляйте /start. Дождитесь подтверждения от бота, что он принял файл. Если подтверждения нет, то что-то не так. Работы, присланные иным способом, не принимаются.
• Дедлайн см. в боте. После дедлайна работы не принимаются кроме случаев наличия уважительной причины.
• Прислать нужно ноутбук в формате ipynb. Если вы строите интерактивные графики, их стоит прислать в формате html.
• Следите за размером файлов. Бот не может принимать файлы весом более 20 Мб. Если файл получается больше, заранее разделите его на несколько.
• Выполнять задание необходимо полностью самостоятельно. При обнаружении списывания всем участникам списывания дается штраф -2 балла к итоговой оценке за семестр.
• Решения, размещенные на каких-либо интернет-ресурсах, не принимаются. Кроме того, публикация решения в открытом доступе может быть приравнена к предоставлении возможности списать.
• Обратите внимание на правила использования ИИ-инструментов при решении домашнего задания.
• Код из рассказанных на занятиях ноутбуков можно использовать без ограничений.
• Для выполнения задания используйте этот ноутбук в качестве основы, ничего не удаляя из него. Можно добавлять необходимое количество ячеек.
• Комментарии к решению пишите в markdown-ячейках.
• Выполнение задания (ход решения, выводы и пр.) должно быть осуществлено на русском языке.
• Решение проверяется системой ИИ-проверки ThetaGrader. Результат проверки валидируется и исправляется человеком, после чего комментарии отправляются студентам.
• Если код будет не понятен проверяющему, оценка может быть снижена.
• Никакой код из данного задания при проверке запускаться не будет. Если код студента не выполнен, недописан и т.д., то он не оценивается.
• Код из рассказанных на занятиях ноутбуков можно использовать без ограничений.

Правила оформления теоретических задач:

• Решения необходимо оформить в виде $\LaTeX$ в markdown-ячейках. Иные способы (в т.ч. фотографии) не принимаются.
• Если вы не знаете $\LaTeX$, используйте ИИ-инструменты для оформления черновика решения. Примеры были показаны на лекции 2 по ИИ-инструментам.
• В решениях поясняйте, чем вы пользуетесь, хотя бы кратко.
• Решение, в котором есть только ответ, и отсутствуют вычисления, оценивается в 0 баллов.

Важно!!! Правила заполнения ноутбука:

• Запрещается удалять имеющиеся в ноутбуке ячейки, менять местами положения задач.
• Сохраняйте естественный линейный порядок повествования в ноутбуке сверху-вниз.
• Отвечайте на вопросы, а также добавляйте новые ячейки в предложенных местах, которые обозначены <...>.
• В markdown-ячейка, содержащих описание задачи, находятся специальные отметки, которые запрещается модифицировать.
• При нарушении данных правил работа может получить 0 баллов.

Баллы за задание:

Легкая часть (достаточно на "хор"):

• Задачи 1-4: скачайте первый ноутбук с условием задания со страницы курса.

Сложная часть (необходимо на "отл"):

• Задача 5 — 80 баллов;
• Задача 6 — 70 баллов;
• Задача 7 — 50 баллов.

Баллы учитываются в обязательной части курса и не влияют на оценку по факультативной части.

# Bot check

# HW_ID: fpmi_ad3_part2
# Бот проверит этот ID и предупредит, если случайно сдать что-то не то.

# Status: not final
# Перед отправкой в финальном решении удали "not" в строчке выше.
# Так бот проверит, что ты отправляешь финальную версию, а не промежуточную.
# Никакие значения в этой ячейке не влияют на факт сдачи работы.

import numpy as np
import pandas as pd
import seaborn as sns
from typing import Literal

sns.set(style="whitegrid", palette="Set2")

При решении задания используйте sklearn. Пропишите сюда необходимые импорты

from sklearn.base import BaseEstimator
from sklearn.preprocessing import StandardScaler

...

---

Легкая часть

Задачи 1-4: скачайте первый ноутбук с условием задания со страницы курса.

---

Сложная часть

---

Задача 5.

1. Реализуйте логистическую регрессию для двух вариантов поиска оценки параметров:

• простой градиентный спуск;
• стохастический градиентный спуск с batch_size элементами на каждой итерации.

Останавливайте итерации при выполнении хотя бы одного из двух условий:

• количество итераций превзошло число max_iter;
• оптимизируемый функционал изменился за итерацию не более чем на tol.

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

Замечания.

1. Для чистоты эксперимента время шага внутри цикла нужно замерять от конца предыдущего шага до конца текущего, а не от начала текущего шага. Время измеряйте с помощью from time import time.

2. Иногда при подсчете сигмоиды и оптимизируемого функционала могут возникать вычислительные ошибки. Для их избежания существуют специальные трюки.

   • How to Evaluate the Logistic Loss and not NaN trying
   • Exp-normalize trick
     

3. Трюки не обязательно реализовывать самостоятельно, можете воспользоваться функциями для них из numpy или scipy:

   • numpy.logaddexp;
   • scipy.special.logsumexp.

4. Обратите внимание, что класс LogisticRegression — наследник класса BaseEstimator, это с легкостью позволит использовать наш класс в различных пайплайнах библиотеки sklearn.

5. Следите за качеством кода, комментируйте логические этапы кода. Несоблюдение этого требования может привести к потере баллов.

# При реализации класса запрещено пользоваться ИИ-инструментами.


class LogisticRegression(BaseEstimator):
    """Модель логистической регрессии.

    Параметры:
    method (Literal['gd', 'sgd']): Метод оптимизации ('gd' - градиентный спуск, 
        'sgd' - стохастический градиентный спуск).
    learning_rate (float): Константа скорости обучения, на которую домножаем градиент при обучении
    tol (float): Допустимое изменение функционала между итерациями.
    max_iter (int): Максимальное число итераций.
    batch_size (int): Размер выборки для оценки градиента (используется только при 'sgd').
    fit_intercept (bool): Добавлять ли константу в признаки.
    save_history (bool): Сохранять ли историю обучения.
    """

    def __init__(
        self,
        method: Literal["gd", "sgd"] = "gd",
        learning_rate: float = 0.5,
        tol: float = 1e-3,
        max_iter: int = int(1e4),
        batch_size: int = 64,
        fit_intercept: bool = True,
        save_history: bool = True,
    ):
        """Создает модель и инициализирует параметры."""
        self.method = method
        self.learning_rate = learning_rate
        self.tol = tol
        self.max_iter = max_iter
        self.batch_size = batch_size
        self.fit_intercept = fit_intercept
        self.save_history = save_history
        self.history = []  # История обучения

    @staticmethod
    def _sigmoid(x: np.ndarray) -> np.ndarray:
        """Вычисляет сигмоидную функцию."""
        return 1 / (1 + np.exp(-x))

    def _add_intercept(self, X: np.ndarray) -> np.ndarray:
        """Добавляет свободный коэффициент к матрице признаков.

        Параметры: X (np.ndarray): Исходная матрица признаков.

        Возвращает: np.ndarray: Матрица X с добавленным свободным
        коэффициентом.
        """
        X_copy = np.full((X.shape[0], X.shape[1] + 1), fill_value=1)
        X_copy[:, :-1] = X
        return X_copy

    def fit(self, X: np.ndarray, Y: np.ndarray) -> "LogisticRegression":
        """Обучает модель логистической регрессии.

        Также, в случае self.save_history=True, добавляет в self.history
        текущее значение оптимизируемого функционала и затраченное время.

        Параметры:
        X (np.ndarray): Матрица признаков.
        Y (np.ndarray): Вектор истинных меток.

        Возвращает:
        LogisticRegression: Обученная модель.
        """
        if X.shape[0] != Y.shape[0]:
            raise ValueError("Количество строк в X и Y должно совпадать")

        if self.fit_intercept:
            X_copy = self._add_intercept(X)
        else:
            X_copy = X.copy()

        ...

        self.coef_ = ...  # Коэффициенты модели
        self.intercept_ = ...  # Свободный коэффициент
        self.n_iter_ = ...  # Число итераций

        return self

    def predict(self, X: np.ndarray) -> np.ndarray:
        """Возвращает предсказанные классы.

        Параметры: X (np.ndarray): Матрица признаков.

        Возвращает: np.ndarray: Предсказанные классы.
        """
        if self.fit_intercept:
            X_copy = self._add_intercept(X)
        else:
            X_copy = X.copy()

        if X_copy.shape[1] != self.coef_.shape[0]:
            raise ValueError("Число признаков в X не соответствует числу коэффициентов модели")

        ...
        return predictions

    def predict_proba(self, X: np.ndarray) -> np.ndarray:
        """Возвращает вероятности классов 0 и 1.

        Параметры: X (np.ndarray): Матрица признаков.

        Возвращает: np.ndarray: Матрица вероятностей классов (n_samples,
        2).
        """
        if self.fit_intercept:
            X_copy = self._add_intercept(X)
        else:
            X_copy = X.copy()

        if X_copy.shape[1] != self.coef_.shape[0]:
            raise ValueError("Число признаков в X не соответствует числу коэффициентов модели")

        ...
        return prob_predictions

Рассмотрим датасет Diabetes Health Indicators.

Для данного задания будем рассматривать версию датасета diabetes_binary_5050split_health_indicators_BRFSS2015.csv

Этот датасет содержит статистику здравоохранения и информацию об образе жизни, полученную в результате опросов вместе с меткой наличия/отсутствия диабета у участников. Среди признаков есть демографические данные, результаты лабораторных тестов и ответы на вопросы анкеты. Целевая переменная Diabetes_binary определяет статус пациента: есть ли у него диабет или предиабет (1), или он здоров (0).

Рассмотрим некоторые признаки, представленные в датасете.

Показатели здоровья

• HighBP: Высокое кровяное давление (1 = да, 0 = нет).

• HighChol: Высокий уровень холестерина (1 = да, 0 = нет).

• CholCheck: Проверка уровня холестерина за последние 5 лет (1 = да, 0 = нет).

• BMI: Индекс массы тела (рассчитывается как вес (кг) / рост² (м²)).

• GenHlth: Общая оценка здоровья (1 = отличное, 2 = очень хорошее, ..., 5 = плохое).

Образ жизни

• Smoker: Статус курения (1 = выкурил ≥100 сигарет за жизнь, 0 = нет).

• PhysActivity: Физическая активность вне работы (1 = да, 0 = нет).

• Fruits: Регулярное употребление фруктов (1 = не менее 1 раз в день, 0 = реже).

Доступ к медицине

• AnyHealthcare: Наличие медицинской страховки (1 = да, 0 = нет).

• NoDocbcCost: Отказ от визита к врачу из-за стоимости (1 = да, 0 = нет).

Скачайте файл и прочитайте его с помощью pandas.

dataset = pd.read_csv("diabets_health_indicators.csv")
dataset.head()

Разделите выборку на обучающую и тестовую и выполните преобразование категориальных признаков.

Для интерпретации коэффициентов необходимо нормализовать данные. Воспользуемся для этого классом StandardScaler из библиотеки sklearn.

scaler = StandardScaler()
...

2. Обучите две модели логистической регрессии с помощью методов

• простой градиентный спуск;
• стохастический градиентный спуск.

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

Замечания:

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

Сделайте выводы. Что будет при обучении на датасете, если увеличить количество объектов, а число признаков оставить прежним?

...

3. Исследуйте влияние размер шага (learning_rate) на качество модели для двух режимов обучения (простой и стохастический градиентный спуск). Для каждого размера шага получите качество модели при использовании простого и стохастического градиентного спуска. Сравните качество полученных моделей по метрике accuracy.

learning_rate_list = np.logspace(-5, 3, 8)

Сделайте выводы

...

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

Какой learning_rate стоит выбирать в зависимости от способа обучения модели? Чем плохи маленькие и большие learning_rate?

...

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

Как можно проинтерпретировать полученные результаты относительно решаемой задачи?

...

5. Сравните данную модель с бейзлайном, который в качестве предсказания выдает самый частый класс на обучающей выборке.

Насколько хорошее получилось качество обученной модели?

...

6. В исходной выборке оставьте два вещественных признака, которые имеют наибольшее влияние на предсказание в предыдущем пункте. Обучите на них модель на 10000 итерациях. Визуализируйте предсказание класса $1$ для нескольких промежуточных итераций.

iters = [10, 20, 500, 1000, 5000, 10000]

Вывод:

...

---

Задача 6.

В этой задаче вам предлагается реализовать регрессию Хьюбера, а также применить ее к данным с выбросами. Для начала реализуйте класс по шаблону снизу. Обратите внимание, что класс HuberRegression — наследник класса BaseEstimator, это с легкостью позволит использовать наш класс в различных пайплайнах библиотеки sklearn.

1. Задача оптимизации для регрессии Хьюбера выглядит следующим образом: $$\sum_{i=1}^n R\left(Y_i - x_i^T\theta\right) \rightarrow \min_\theta,$$ где $R(x)$ — функция потерь Хьюбера, определяемая как $$R(x) = \frac{x^2}{2} I\left\{|x| < c\right\} + c \left(|x| - \frac{c}{2}\right)I\left\{|x| > c\right\}.$$

Выпишите формулы для градиентного и стохастического градиентного спусков. В чем польза такой функции потерь?

...

2. Реализуем теперь класс. При реализации класса запрещено пользоваться ИИ-инструментами.

# При реализации класса запрещено пользоваться ИИ-инструментами.


class HuberRegression(BaseEstimator):
    """Класс, реализующий линейную регрессию с функцией потерь Хьюбера."""

    def __init__(self, c: float = 1.0, fit_intercept: bool = True, max_iter: int = 1000) -> None:
        """Инициализирует модель.

        Параметры: c (float): Константа из функции потерь Хьюбера.
        fit_intercept (bool): Добавлять ли константный признак. max_iter
        (int): Максимальное число итераций оптимизации.
        """
        self.c = c
        self.fit_intercept = fit_intercept
        self.max_iter = max_iter

    def fit(self, X: np.ndarray, y: np.ndarray) -> "HuberRegression":
        """Обучает модель.

        Параметры:
        X (np.ndarray): Матрица признаков.
        y (np.ndarray): Вектор целевой переменной.

        Возвращает:
        HuberRegression: Обученная модель.
        """
        if X.shape[0] != y.shape[0]:
            raise ValueError("Количество строк в X и y должно совпадать")

        ...

        self.coef_ = ...  # Коэффициенты модели
        self.intercept_ = ...  # Свободный коэффициент
        self.n_iter_ = ...  # Число итераций

        return self

    def predict(self, X: np.ndarray) -> np.ndarray:
        """Делает предсказание на новых данных.

        Параметры: X (np.ndarray): Матрица признаков.

        Возвращает: np.ndarray: Вектор предсказанных значений.
        """

        if X_copy.shape[1] != self.coef_.shape[0]:
            raise ValueError("Число признаков в X не соответствует числу коэффициентов модели")

        ...
        return pred

3. Загрузите данные из файлов train.csv, test.csv. Не забудьте, что всю аналитику, а также процесс обучения и подбор гиперпараметров необходимо выполнять на обучающей выборке.

...

Посмотрите на зависимость целевой переменной от каждого признака.

...

Что можно сказать о наличии возможных выбросов? Какое влияние они могут оказать?

...

4. Обучите простую линейную регрессию и посчитайте качество на тестовой выборке по метрике MSE.

...

Что можно сказать о качестве нашей модели?

...

5. Теперь обучите линейную регресcию Хьюбера и посчитайте качество на тестовой части по метрикe MSE.

...

Что изменилось?

...

6. Для обучающей выборки постройте два графика (по графику на каждую модель), на которых изобразите зависимость истинного и предсказанного значения таргета от каждого признака.

...

Что можно заметить на этих графиках?

...

7. Обучите регрессию Хьюбера на данных из задачи 2 и сравните качество модели с простой линейной регрессией, которую вы построили в задаче 2.

...

Вывод:

...

---

Задача 7.

Рассмотрим модель одномерной регрессии $y(x) = \theta x$, где $x \in \mathbb{R}$ — одномерный признак, $y \in \mathbb{R}$ — целевой признак, $\theta \in \mathbb{R}$ — неизвестный параметр. Имеется выборка размера $n$, полученная по правилу $$Y_i = \theta x_i + \varepsilon_i,\ \ \ i=1,...,n,$$ где $\varepsilon_i$ — случайная ошибка измерений.

Предложите точный алгоритм поиска оценки параметра $\theta$ методом наименьших модулей, то есть $$\sum_{i=1}^n \left|Y_i - \theta x_i\right| \to \min_\theta,$$ работающий за время $O(n \log n)$. Приведите его описание и теоретическое обоснование. Реализация в коде не требуется.

Решение:

...