Введение
В этом руководстве показано, как выполнить квантильную регрессию с использованием scikit-learn. Мы сгенерируем два синтетических набора данных, чтобы показать, как квантильная регрессия может предсказывать не тривиальные условные квантили. Мы будем использовать класс QuantileRegressor для оценки медианы, а также нижнего и верхнего квантилей, установленных соответственно на 5% и 95%. Мы сравним QuantileRegressor с LinearRegression и оценим их производительность с использованием средней абсолютной ошибки (MAE) и средней квадратичной ошибки (MSE).
Советы по работе с ВМ
После запуска ВМ нажмите в левом верхнем углу, чтобы переключиться на вкладку Notebook и получить доступ к Jupyter Notebook для практики.
Иногда вам может потребоваться подождать несколько секунд, пока Jupyter Notebook не загрузится полностью. Валидация операций не может быть автоматизирована из-за ограничений Jupyter Notebook.
Если вы сталкиваетесь с проблемами при обучении, не стесняйтесь обращаться к Labby. Оставьте отзыв после занятия, и мы оперативно решим проблему для вас.
Генерация набора данных
Мы сгенерируем два синтетических набора данных с одинаковым математическим ожиданием, используя линейную зависимость с одной характеристикой 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_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")
Квантильная регрессия
Мы будем использовать класс QuantileRegressor для оценки медианы, а также нижнего и верхнего квантилей, установленных соответственно на 5% и 95%. Мы будем использовать квантили 5% и 95% для поиска выбросов в обучающей выборке за пределами центрального интервала в 90%.
from sklearn.linear_model import QuantileRegressor
## Эта строка нужна, чтобы избежать несовместимости при использовании старой версии SciPy.
## Вы должны использовать `solver="highs"` с более новой версией SciPy.
solver = "highs" if sp_version >= parse_version("1.6.0") else "interior-point"
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, solver=solver)
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
)
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")
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, solver=solver)
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")
Сравнение QuantileRegressor и LinearRegression
Мы сравним QuantileRegressor и LinearRegression и оценим их производительность с использованием средней абсолютной ошибки (MAE) и средней квадратичной ошибки (MSE).
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import cross_validate
linear_regression = LinearRegression()
quantile_regression = QuantileRegressor(quantile=0.5, alpha=0, solver=solver)
y_pred_lr = linear_regression.fit(X, y_pareto).predict(X)
y_pred_qr = quantile_regression.fit(X, y_pareto).predict(X)
print(f"""Training error (in-sample performance)
{linear_regression.__class__.__name__}:
MAE = {mean_absolute_error(y_pareto, y_pred_lr):.3f}
MSE = {mean_squared_error(y_pareto, y_pred_lr):.3f}
{quantile_regression.__class__.__name__}:
MAE = {mean_absolute_error(y_pareto, y_pred_qr):.3f}
MSE = {mean_squared_error(y_pareto, y_pred_qr):.3f}
""")
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(f"""Test error (cross-validated performance)
{linear_regression.__class__.__name__}:
MAE = {-cv_results_lr["test_neg_mean_absolute_error"].mean():.3f}
MSE = {-cv_results_lr["test_neg_mean_squared_error"].mean():.3f}
{quantile_regression.__class__.__name__}:
MAE = {-cv_results_qr["test_neg_mean_absolute_error"].mean():.3f}
MSE = {-cv_results_qr["test_neg_mean_squared_error"].mean():.3f}
""")
Резюме
В этом уроке мы узнали, как выполнять квантильную регрессию с использованием scikit - learn. Мы сгенерировали два синтетических набора данных, чтобы показать, как квантильная регрессия может предсказывать нетривиальные условные квантили. Мы использовали класс QuantileRegressor для оценки медианы, а также нижнего и верхнего квантилей, установленных соответственно на 5% и 95%. Мы сравнили QuantileRegressor с LinearRegression и оценили их производительность с использованием средней абсолютной ошибки (MAE) и средней квадратичной ошибки (MSE).