Квантильная регрессия

Генерация набора данных

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

import numpy as np

rng = np.random.RandomState(42)
x = np.linspace(start=0, stop=10, num=100)
X = x[:, np.newaxis]
y_true_mean = 10 + 0.5 * x

Мы создадим две задачи, изменив распределение целевой переменной y, сохранив при этом одинаковое математическое ожидание:

• в первом случае добавляется гетероскедастичный нормальный шум;
• во втором случае добавляется асимметричный шум с распределением Парето.

y_normal = y_true_mean + rng.normal(loc=0, scale=0.5 + 0.5 * x, size=x.shape[0])
a = 5
y_pareto = y_true_mean + 10 * (rng.pareto(a, size=x.shape[0]) - 1 / (a - 1))

Для начала визуализируем наборы данных и распределение остатков y - mean(y).

import matplotlib.pyplot as plt

_, axs = plt.subplots(nrows=2, ncols=2, figsize=(15, 11), sharex="row", sharey="row")

axs[0, 0].plot(x, y_true_mean, label="True mean")
axs[0, 0].scatter(x, y_normal, color="black", alpha=0.5, label="Observations")
axs[1, 0].hist(y_true_mean - y_normal, edgecolor="black")


axs[0, 1].plot(x, y_true_mean, label="True mean")
axs[0, 1].scatter(x, y_pareto, color="black", alpha=0.5, label="Observations")
axs[1, 1].hist(y_true_mean - y_pareto, edgecolor="black")

axs[0, 0].set_title("Dataset with heteroscedastic Normal distributed targets")
axs[0, 1].set_title("Dataset with asymmetric Pareto distributed target")
axs[1, 0].set_title(
    "Residuals distribution for heteroscedastic Normal distributed targets"
)
axs[1, 1].set_title("Residuals distribution for asymmetric Pareto distributed target")
axs[0, 0].legend()
axs[0, 1].legend()
axs[0, 0].set_ylabel("y")
axs[1, 0].set_ylabel("Counts")
axs[0, 1].set_xlabel("x")
axs[0, 0].set_xlabel("x")
axs[1, 0].set_xlabel("Residuals")
_ = axs[1, 1].set_xlabel("Residuals")

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

Для целевой переменной с асимметричным распределением Парето мы видим, что положительные остатки ограничены сверху.

Такие зашумлённые целевые переменные затрудняют оценивание с помощью LinearRegression: для получения стабильных результатов нам потребуется больше данных, а крупные выбросы могут существенно повлиять на подобранные коэффициенты. (Иными словами, при постоянной дисперсии оценки методом наименьших квадратов с увеличением объёма выборки значительно быстрее сходятся к истинным коэффициентам.)

В таких асимметричных условиях медиана или различные квантили дают дополнительную информацию. Кроме того, оценка медианы гораздо устойчивее к выбросам и распределениям с тяжёлыми хвостами. Однако следует учитывать, что экстремальные квантили оцениваются по очень малому числу точек данных. Квантиль 95% оценивается примерно по 5% наибольших значений и поэтому также в некоторой степени чувствителен к выбросам.

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

Обучение QuantileRegressor

В этом разделе мы хотим оценить условную медиану, а также нижний и верхний квантили, равные соответственно 5% и 95%. Таким образом, мы получим три линейные модели — по одной для каждого квантиля.

Мы воспользуемся квантилями 5% и 95%, чтобы найти выбросы в обучающей выборке, расположенные за пределами центрального интервала в 90%.

from sklearn.linear_model import QuantileRegressor

quantiles = [0.05, 0.5, 0.95]
predictions = {}
out_bounds_predictions = np.zeros_like(y_true_mean, dtype=np.bool_)
for quantile in quantiles:
    qr = QuantileRegressor(quantile=quantile, alpha=0)
    y_pred = qr.fit(X, y_normal).predict(X)
    predictions[quantile] = y_pred

    if quantile == min(quantiles):
        out_bounds_predictions = np.logical_or(
            out_bounds_predictions, y_pred >= y_normal
        )
    elif quantile == max(quantiles):
        out_bounds_predictions = np.logical_or(
            out_bounds_predictions, y_pred <= y_normal
        )

Теперь можно построить три линейные модели и выделить точки, которые находятся внутри центрального 90%-го интервала, отдельно от точек за его пределами.

plt.plot(X, y_true_mean, color="black", linestyle="dashed", label="True mean")

for quantile, y_pred in predictions.items():
    plt.plot(X, y_pred, label=f"Quantile: {quantile}")

plt.scatter(
    x[out_bounds_predictions],
    y_normal[out_bounds_predictions],
    color="black",
    marker="+",
    alpha=0.5,
    label="Outside interval",
)
plt.scatter(
    x[~out_bounds_predictions],
    y_normal[~out_bounds_predictions],
    color="black",
    alpha=0.5,
    label="Inside interval",
)

plt.legend()
plt.xlabel("x")
plt.ylabel("y")
_ = plt.title("Quantiles of heteroscedastic Normal distributed target")

Поскольку шум по-прежнему имеет нормальное, в частности симметричное, распределение, истинные условное среднее и условная медиана совпадают. Действительно, мы видим, что оценка медианы почти совпадает с истинным средним. Мы наблюдаем влияние увеличения дисперсии шума на квантили 5% и 95%: наклоны этих квантилей сильно различаются, а интервал между ними расширяется при увеличении x.

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

Повторим тот же эксперимент для целевой переменной с асимметричным распределением Парето.

quantiles = [0.05, 0.5, 0.95]
predictions = {}
out_bounds_predictions = np.zeros_like(y_true_mean, dtype=np.bool_)
for quantile in quantiles:
    qr = QuantileRegressor(quantile=quantile, alpha=0)
    y_pred = qr.fit(X, y_pareto).predict(X)
    predictions[quantile] = y_pred

    if quantile == min(quantiles):
        out_bounds_predictions = np.logical_or(
            out_bounds_predictions, y_pred >= y_pareto
        )
    elif quantile == max(quantiles):
        out_bounds_predictions = np.logical_or(
            out_bounds_predictions, y_pred <= y_pareto
        )

plt.plot(X, y_true_mean, color="black", linestyle="dashed", label="True mean")

for quantile, y_pred in predictions.items():
    plt.plot(X, y_pred, label=f"Quantile: {quantile}")

plt.scatter(
    x[out_bounds_predictions],
    y_pareto[out_bounds_predictions],
    color="black",
    marker="+",
    alpha=0.5,
    label="Outside interval",
)
plt.scatter(
    x[~out_bounds_predictions],
    y_pareto[~out_bounds_predictions],
    color="black",
    alpha=0.5,
    label="Inside interval",
)

plt.legend()
plt.xlabel("x")
plt.ylabel("y")
_ = plt.title("Quantiles of asymmetric Pareto distributed target")

Из-за асимметрии распределения шума истинное среднее и оценка условной медианы различаются. Мы также видим, что у каждой квантильной модели разные параметры, позволяющие лучше подобрать нужный квантиль. В идеале в этом случае все квантили были бы параллельны; это было бы заметнее при большем количестве точек данных или менее экстремальных квантилях, например 10% и 90%.

Сравнение QuantileRegressor и LinearRegression

В этом разделе мы подробнее рассмотрим различия между функциями потерь, которые минимизируют QuantileRegressor и LinearRegression.

Действительно, LinearRegression использует метод наименьших квадратов, минимизируя среднеквадратичную ошибку (MSE) между целевыми значениями в обучающей выборке и прогнозами. В отличие от него, QuantileRegressor с quantile=0.5 минимизирует среднюю абсолютную ошибку (MAE).

Почему это важно? Функции потерь определяют, что именно стремится прогнозировать модель; см. руководство пользователя по выбору функции оценки. Кратко говоря, модель, минимизирующая MSE, прогнозирует среднее (математическое ожидание), а модель, минимизирующая MAE, прогнозирует медиану.

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

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, mean_squared_error

linear_regression = LinearRegression()
quantile_regression = QuantileRegressor(quantile=0.5, alpha=0)

y_pred_lr = linear_regression.fit(X, y_pareto).predict(X)
y_pred_qr = quantile_regression.fit(X, y_pareto).predict(X)

print(
    "Training error (in-sample performance)\n"
    f"{'model':<20}   MAE   MSE\n"
    f"{linear_regression.__class__.__name__:<20} "
    f"{mean_absolute_error(y_pareto, y_pred_lr):5.3f} "
    f"{mean_squared_error(y_pareto, y_pred_lr):5.3f}\n"
    f"{quantile_regression.__class__.__name__:<20} "
    f"{mean_absolute_error(y_pareto, y_pred_qr):5.3f} "
    f"{mean_squared_error(y_pareto, y_pred_qr):5.3f}"
)

Training error (in-sample performance)
model                  MAE   MSE
LinearRegression     1.805 6.486
QuantileRegressor    1.670 7.025

На обучающей выборке MAE у QuantileRegressor ниже, чем у LinearRegression. И наоборот, MSE у LinearRegression ниже, чем у QuantileRegressor. Эти результаты подтверждают, что MAE — функция потерь, минимизируемая QuantileRegressor, а MSE минимизируется LinearRegression.

Аналогичную оценку можно провести, рассмотрев ошибку на тестовой выборке, полученную с помощью кросс-валидации.

from sklearn.model_selection import cross_validate

cv_results_lr = cross_validate(
    linear_regression,
    X,
    y_pareto,
    cv=3,
    scoring=["neg_mean_absolute_error", "neg_mean_squared_error"],
)
cv_results_qr = cross_validate(
    quantile_regression,
    X,
    y_pareto,
    cv=3,
    scoring=["neg_mean_absolute_error", "neg_mean_squared_error"],
)
print(
    "Test error (cross-validated performance)\n"
    f"{'model':<20}   MAE   MSE\n"
    f"{linear_regression.__class__.__name__:<20} "
    f"{-cv_results_lr['test_neg_mean_absolute_error'].mean():5.3f} "
    f"{-cv_results_lr['test_neg_mean_squared_error'].mean():5.3f}\n"
    f"{quantile_regression.__class__.__name__:<20} "
    f"{-cv_results_qr['test_neg_mean_absolute_error'].mean():5.3f} "
    f"{-cv_results_qr['test_neg_mean_squared_error'].mean():5.3f}"
)

Test error (cross-validated performance)
model                  MAE   MSE
LinearRegression     1.732 6.690
QuantileRegressor    1.679 7.129

При оценке на данных, не участвовавших в обучении, мы приходим к аналогичным выводам.

Общее время выполнения скрипта: (0 минут 0.506 секунд)

Download Jupyter notebook: plot_quantile_regression.ipynb

Download Python source code: plot_quantile_regression.py

Download zipped: plot_quantile_regression.zip

Связанные примеры

Интервалы прогнозирования для градиентного бустинга в задаче регрессии

Основные нововведения scikit-learn 1.1

Признаки с лагами для прогнозирования временных рядов

Демонстрация предположений метода k-means