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

(Фотография была выбрана как "гиперсамолет", каламбур)

В предыдущем посте я поделился некоторыми ресурсами, специально ориентированными на новичков в машинном обучении. — Путь самоучки к AI/ML, часть 1. Я очень рекомендую ML Zoomcamp. Хотя он предполагает предыдущий опыт работы с Python и некоторыми пакетами, связанными с ML, такими как Numpy и Pandas, он в основном знакомит студента с ML с прикладной точки зрения с помощью практических проектов. Учебную программу можно считать завершенной, поскольку она включает в себя проекты, начиная от линейной регрессии и классификации и заканчивая развертыванием моделей и глубоким обучением.

В этом посте мы рассмотрим первый простой пример линейной регрессии, который также является одним из первых проектов курса. Проблема заключается в прогнозировании цен на жилье на основе известного набора данных, который можно найти на Kaggle (Цены на жилье в Калифорнии).

Ниже приведена схема и реализация шагов:

  1. Подготовка данных — получаем данные и выбираем, что будем использовать
  2. EDA — мы выполняем простое исследование и наблюдение за данными
  3. Фреймворк Train-Test — разделяем наши данные
  4. Линейная регрессия — модель
  5. Регуляризация и настройка модели — выберите лучшую модель для наших данных

Подготовка данных

Первым шагом является получение данных и проверка их характеристик. Библиотека Pandas очень удобна для таких проблем с табличными данными. Сначала мы импортируем готовую к использованию библиотеку, а затем читаем файл csv, который мы уже скачали ранее с веб-сайта Kaggle (или любого другого источника данных, который нам подходит).

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

import pandas as pd
import matplotlib.pyplot as plt
# Read the csv file containing the data
data = pd.read_csv('housing.csv')
# Print the shape of the dataframe
print("data shape: ",data.shape)
# Print the first few lines of the dataframe
data.head()

форма данных: (20640, 10)

Из всех столбцов, которые у нас есть в нашем наборе данных, следующий вопрос:

Какое количество мы пытаемся предсказать?

В нашей задаче цель состоит в том, чтобы предсказать цены на жилье, как они появляются в медианной стоимости дома. Тот факт, что прогнозируемое значение является непрерывным числом, делает его проблемой регрессии. Очень простой, но весьма полезный первый подход — линейная регрессия. Несмотря на свою простоту и тот факт, что она существует уже несколько столетий, линейная регрессия часто оказывается рабочей лошадкой во многих отраслях (в той или иной форме) и, безусловно, одним из первых подходов, которые необходимо протестировать при разработке и определении базовой линии.

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

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

Тем не менее, мы строим распределение целевого значения.

plt.hist(data['median_house_value'], bins=50)

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

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

# select the columns
cols = ['latitude','longitude','housing_median_age','total_rooms','total_bedrooms','population','households','median_income','median_house_value']
data = data[cols]

Исследовательский анализ данных (EDA)

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

Мы начнем с проверки того, есть ли у нас какие-либо значения Null (также известные как NaN, n/a и т. д.). Это важный первый шаг, поскольку все алгоритмы, которые мы используем в машинном обучении, требуют определенной доступности данных, чтобы они могли работать должным образом. Алгоритму нужно знать данные, с которыми он должен работать, потому что… ну, это абстрактная конструкция, которая имеет очень четко определенную реализацию в виде компьютерной программы.

Строка ниже выводит количество записей, которые имеют значение Null для каждой из функций (столбцов) нашего набора данных.

data.isnull().sum()

Если обнаружены какие-либо значения Null, мы можем решить, как мы будем обрабатывать отсутствующие значения на следующем шаге. Мы видим, что здесь у нас есть 207 пропущенных значений в столбце «total_bedrooms».

Рекомендуемый шаг — также проверить типы данных функций и некоторую дополнительную информацию, относящуюся к числовым функциям (функциям, которые в основном представлены в виде чисел). В отличие от числовых признаков категориальные признаки представляют собой качественную характеристику данных. В нашем случае исходный признак «ocean_proximity» является категориальным, поскольку он качественно представляет, где дом расположен по отношению к морю. Для исходного набора данных (включая столбец «ocean_proximity») у нас было бы:

data.info()

Мы видим, что «ocean_proximity» имеет тип object, что означает, что он имеет тип данных, который нельзя легко рассматривать как число. Обратите внимание, что строки могут иметь тип объекта, и обычно стоит явно проверить, не ошибочно ли он назначен как объект, хотя он должен быть числом. Обратите также внимание на количество ненулевых значений в столбце «total_bedrooms», которое меньше, чем общее количество записей 20 640.

# print some statistics
data.describe()

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

Настройте платформу для проверки и проверки поезда.

Золотое правило структуры train-validation-test очень простое, и на практике его можно сформулировать так:

Как только вы получите данные, отделите тестовый набор и забудьте об этом!

Это действительно избавит вас от многих разочарований и ошибочных подходов к данным, как только вы это сделаете. Обоснование этого довольно простое. Вы же не хотите утечки информации между поездом и тестовыми наборами. Единственная цель тестового набора — смоделировать реальный сценарий, в котором ваша модель не будет иметь доступа к будущим данным, с которыми она может столкнуться. Следовательно, вы не хотите, чтобы обучение вашей модели каким-либо образом зависело от набора тестов.

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

Тем не менее, золотое правило — это довольно простой и общий подход. Как только вы сделаете честное разделение, где вы убедитесь, что данные тщательно перетасованы, а ваши обучающие примеры максимально приближаются к состоянию i.i.d (независимым и одинаково распределенным), вы забываете о тестовых данных до самого конца, где вы сделаете одну и единственную последнюю проверку модели.

Мы можем создать функцию, которая будет использоваться для разделения данных и возврата набора поездов, набора проверки и набора тестов. Приведенная ниже функция представляет собой простую и не самую эффективную реализацию функции разделения. Типичным коэффициентом разделения будет 60% данных обучения, 20% данных проверки и 20% данных испытаний.

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

def split_data(data, train_ratio = 0.6, test_ratio =0.2, val_ratio = 0.2):
    # we create a random shuffle of the dataset indexes
    idx = np.arange(data.shape[0])
    np.random.shuffle(idx)
    # we find the indexes that we'll split
    train = int(train_ratio*len(idx))
    test = int(test_ratio*len(idx))
    val = int(val_ratio*len(idx))
    
    assert train_ratio + val_ratio + test_ratio == 1
    y = data['median_house_value'].copy()
    X = data.drop(['median_house_value'], axis=1).copy()
    
    # we create our train, validation and test datasets
    X_train = X.iloc[:train]
    X_val = X.iloc[train:train+val]
    X_test = X.iloc[-test:]
    y_train = np.log1p(y.iloc[:train])
    y_val = np.log1p(y.iloc[train:train+val])
    y_test = np.log1p(y.iloc[-test:])
    return X_train, X_val, X_test, y_train, y_val, y_test

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

mean_bedrooms = X_train['total_bedrooms'].mean()
X_train_mean = X_train.fillna(mean_bedrooms)

Линейная регрессия в векторной форме, нормальное уравнение

Следующим шагом является построение модели. В нашем случае модель представляет собой простую линейную регрессию. Мы строим модель линейной регрессии в векторной форме, также известную как нормальное уравнение, которая является решением задачи наименьших квадратов (см. Обычные наименьшие квадраты — Википедия). При условии, что столбцы матрицы X линейно независимы и количество признаков не создает переопределенной ситуации, так что обращение матриц возможно,

предполагая решение вида:

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

def train_linear_regression(X, y):
    ones = np.ones(X.shape[0])
    X = np.column_stack([ones, X])
    XTX = X.T @ X
    XTX_inv = np.linalg.inv(XTX)
    w = XTX_inv @ (X.T) @ (y)
    return w[0], w[1:]

Сначала мы добавляем столбец единиц к данным, которые соответствуют постоянному члену в линейном уравнении. Мы вычисляем матрицу Грама X^T*X, а затем ее обратную, используя функции линейной алгебры, предоставляемые Numpy. В качестве последнего шага мы вычисляем вектор коэффициентов w в соответствии с приведенным выше уравнением.

Хорошо известно, что этот метод не самый эффективный, особенно для многомерных пространств, потому что обращение матриц имеет временную сложность O(n³) по количеству признаков. Чтобы преодолеть эту высокую вычислительную сложность, обычно используются численные итерационные методы, такие как стохастический градиентный спуск и т. д.

Базовая модель

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

Другой моделью, которая может выступать в качестве базовой линии, будет модель, рассчитанная по нормальному уравнению. Как мы увидим позже, мы можем создавать модели, содержащие гиперпараметры, которые можно тонко настроить для достижения наилучшей производительности. Затем эти модели сравнивают с исходными и выбирают лучшую.

Метрики ошибок и производительность

Но как мы оцениваем наши модели? Как мы устанавливаем наш базовый уровень и как мы их сравниваем?

Для этого нам нужна метрика производительности. В случае линейной регрессии общей метрикой является среднеквадратическая ошибка (RMSE). Это просто (начиная слева) корень из среднего квадрата разницы между прогнозируемым целевым значением и фактическим целевым значением. Другими словами, мы сначала вычисляем различия прогнозируемых значений и фактических значений для каждой из наших выборок данных, затем возводим эти различия в квадрат, берем их среднее значение (для невзвешенного среднего это просто суммирование и деление на количество выборок) и, наконец, мы извлекаем из него квадратный корень. Таким образом, RMSE измеряется в тех же единицах, что и целевая переменная y.

def RMSE(y_true, y_pred):
    return (np.sum((y_true - y_pred)**2) / y_true.shape[0]) ** 0.5

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

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

# baseline model
w0, w = train_linear_regression(X_train_mean, y_train)

Наша модель в основном представляет собой коэффициенты w0 и w, где w0 — константа, а w — вектор из 8 чисел, столько же признаков (столбцов), которые используются для прогнозирования.

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

y_pred = w0 + X_val_mean @ w

Где X_val_mean — это набор данных проверки, в котором отсутствующие значения заполняются средним значением соответствующего признака в X_train.

X_val_mean = X_val.fillna(mean_bedrooms)

Теперь пришло время рассчитать RMSE, для этого мы просто делаем:

RMSE(y_val, y_pred)

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

Разработка функций

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

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

Регуляризация

Часто бывает так, что модель подходит к обучающим данным. То есть модель обучается, запоминая данные. Это, конечно, нежелательно, поскольку мы хотим, чтобы наша модель могла хорошо обобщать невидимые данные. Другими словами, мы хотим, чтобы наша модель работала относительно хорошо даже за счет снижения производительности в диапазоне различных значений функций. Это распространенная проблема, когда модель слишком сложна или у нас очень мало примеров данных для обучения по сравнению с количеством функций. Интуитивно говоря, наличие большого количества функций заставляет модель находить конкретную комбинацию wкомпонентов, которая минимизирует потери. Система недоопределена, что означает, что в ней больше неизвестных (компонентов вектора w), чем уравнений (обучающих примеров, удовлетворяющих линейным уравнениям). Недоопределенная система может иметь несколько решений (конкретные комбинации значений компонентов w).

На помощь приходят методы регуляризации! Идея, стоящая за этим, состоит в том, чтобы каким-то образом наказать вектор неизвестных, добавив дополнительный член к функции потерь, которую алгоритм пытается минимизировать. Обычно этот штраф представляет собой норму L² или L¹, что просто означает, что мы либо пытаемся минимизировать евклидову норму вектора w, либо сумму абсолютных значений его компонентов. Важно помнить, что штраф — это, по сути, ограничение на возможные w-векторы, которые мы можем иметь, что связано с необходимостью поддерживать общий убыток на оптимальном уровне.

В подходе с нормальным уравнением может происходить переобучение из-за линейно зависимых признаков, что может привести к тому, что матрица Грама будет необратимой. Регуляризация L² проявляется как дополнительный член, включающий единичную матрицу. Не слишком углубляясь в математику, это делает обращение матрицы Грама всегда возможным, а это означает, что мы всегда можем найти решение.

def train_linear_regression_reg(X, y, r=0.0):
    ones = np.ones(X.shape[0])
    X = np.column_stack([ones, X])
    XTX = X.T.dot(X)
    reg = r * np.eye(XTX.shape[0])
    XTX = XTX + reg
    XTX_inv = np.linalg.inv(XTX)
    w = XTX_inv.dot(X.T).dot(y)
    
    return w[0], w[1:]

Настройка модели

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

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

for r in [0, 0.000001, 0.0001, 0.001, 0.01, 0.1, 1, 5, 10]:
    w0, w = train_linear_regression_reg(X_train_zero, y_train, r)
    y_pred = w0 + X_val_zero.dot(w)
    print(f'{r:10.6f} {RMSE(y_pred, y_val):8.4f}')

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

Использование модели

Наконец-то мы готовы использовать модель! Мы просто вставляем значения признаков в линейное уравнение.

y = w0 + X @ w

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

Заключение

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

Надеюсь, вы нашли это полезным.

Следите за новыми статьями по ML!

Удачного машинного обучения!

LinkedIn: Мариос Кокмотос | LinkedIn