Регрессия по квантилям

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

Для иллюстрации поведения регрессии по квантилям мы сгенерируем два синтетических набора данных. Истинные генеративные случайные процессы для обоих наборов данных будут состоять из одного и того же математического ожидания с линейной зависимостью от одного признака 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.

Построение регрессии по квантилям

В этом разделе мы хотим оценить условную медиану, а также низкий и высокий квантили, установленные соответственно на 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%.