注意
跳转至末尾以下载完整示例代码，或通过 JupyterLite 或 Binder 在浏览器中运行此示例
泊松回归与非正态损失#
本示例说明了在 法国机动车第三方责任索赔数据集 上使用对数线性泊松回归，并将其与使用常用最小二乘误差拟合的线性模型以及使用泊松损失（和对数连接）拟合的非线性 GBRT 模型进行比较。
一些定义
保单是保险公司与个人之间的合同：此处的投保人即车辆驾驶员。
索赔是投保人向保险公司提出的，要求赔偿保险承保范围内损失的请求。
风险敞口是给定保单的保险覆盖持续时间，以年为单位。
索赔频率是索赔数量除以风险敞口，通常以每年索赔数量来衡量。
在此数据集中，每个样本对应一份保险保单。可用特征包括驾驶员年龄、车辆年龄、车辆功率等。
我们的目标是根据投保人群的历史数据，预测新投保人因交通事故造成的索赔预期频率。
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
法国机动车第三方责任索赔数据集#
让我们从 OpenML 加载机动车索赔数据集：https://www.openml.org/d/41214
from sklearn.datasets import fetch_openml
df = fetch_openml(data_id=41214, as_frame=True).frame
df
索赔数量 (ClaimNb) 是一个正整数，可以建模为泊松分布。它被假定为在给定时间间隔内（Exposure，以年为单位）以恒定速率发生的离散事件的数量。
在这里，我们希望通过一个（缩放的）泊松分布，在 X 的条件下，对频率 y = ClaimNb / Exposure 进行建模，并使用 Exposure 作为 sample_weight。
df["Frequency"] = df["ClaimNb"] / df["Exposure"]
print(
"Average Frequency = {}".format(np.average(df["Frequency"], weights=df["Exposure"]))
)
print(
"Fraction of exposure with zero claims = {0:.1%}".format(
df.loc[df["ClaimNb"] == 0, "Exposure"].sum() / df["Exposure"].sum()
)
)
fig, (ax0, ax1, ax2) = plt.subplots(ncols=3, figsize=(16, 4))
ax0.set_title("Number of claims")
_ = df["ClaimNb"].hist(bins=30, log=True, ax=ax0)
ax1.set_title("Exposure in years")
_ = df["Exposure"].hist(bins=30, log=True, ax=ax1)
ax2.set_title("Frequency (number of claims per year)")
_ = df["Frequency"].hist(bins=30, log=True, ax=ax2)
Average Frequency = 0.10070308464041304
Fraction of exposure with zero claims = 93.9%
剩余的列可用于预测索赔事件的频率。这些列非常异构，混合了不同尺度、可能分布非常不均匀的类别变量和数值变量。
为了使用这些预测器拟合线性模型，因此需要进行标准特征转换，如下所示：
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import (
FunctionTransformer,
KBinsDiscretizer,
OneHotEncoder,
StandardScaler,
)
log_scale_transformer = make_pipeline(
FunctionTransformer(np.log, validate=False), StandardScaler()
)
linear_model_preprocessor = ColumnTransformer(
[
("passthrough_numeric", "passthrough", ["BonusMalus"]),
(
"binned_numeric",
KBinsDiscretizer(
n_bins=10, quantile_method="averaged_inverted_cdf", random_state=0
),
["VehAge", "DrivAge"],
),
("log_scaled_numeric", log_scale_transformer, ["Density"]),
(
"onehot_categorical",
OneHotEncoder(),
["VehBrand", "VehPower", "VehGas", "Region", "Area"],
),
],
remainder="drop",
)
常数预测基线#
值得注意的是，超过 93% 的投保人没有索赔。如果我们将此问题转换为二分类任务，它将严重不平衡，即使是只预测均值的简单模型也能达到 93% 的准确率。
为了评估所用指标的合理性，我们将一个“哑”估计器作为基线，该估计器始终预测训练样本的平均频率。
from sklearn.dummy import DummyRegressor
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
df_train, df_test = train_test_split(df, test_size=0.33, random_state=0)
dummy = Pipeline(
[
("preprocessor", linear_model_preprocessor),
("regressor", DummyRegressor(strategy="mean")),
]
).fit(df_train, df_train["Frequency"], regressor__sample_weight=df_train["Exposure"])
让我们使用 3 种不同的回归指标计算此常数预测基线的性能
from sklearn.metrics import (
mean_absolute_error,
mean_poisson_deviance,
mean_squared_error,
)
def score_estimator(estimator, df_test):
"""Score an estimator on the test set."""
y_pred = estimator.predict(df_test)
print(
"MSE: %.3f"
% mean_squared_error(
df_test["Frequency"], y_pred, sample_weight=df_test["Exposure"]
)
)
print(
"MAE: %.3f"
% mean_absolute_error(
df_test["Frequency"], y_pred, sample_weight=df_test["Exposure"]
)
)
# Ignore non-positive predictions, as they are invalid for
# the Poisson deviance.
mask = y_pred > 0
if (~mask).any():
n_masked, n_samples = (~mask).sum(), mask.shape[0]
print(
"WARNING: Estimator yields invalid, non-positive predictions "
f" for {n_masked} samples out of {n_samples}. These predictions "
"are ignored when computing the Poisson deviance."
)
print(
"mean Poisson deviance: %.3f"
% mean_poisson_deviance(
df_test["Frequency"][mask],
y_pred[mask],
sample_weight=df_test["Exposure"][mask],
)
)
print("Constant mean frequency evaluation:")
score_estimator(dummy, df_test)
Constant mean frequency evaluation:
MSE: 0.564
MAE: 0.189
mean Poisson deviance: 0.625
（广义）线性模型#
我们首先使用（L2 惩罚的）最小二乘线性回归模型（更常见的是岭回归）来建模目标变量。我们使用较低的惩罚参数 alpha，因为我们预计这样的线性模型在此类大型数据集上会出现欠拟合。
泊松偏差无法根据模型预测的非正值进行计算。对于确实返回少量非正预测值的模型（例如 Ridge），我们忽略相应的样本，这意味着获得的泊松偏差是近似值。另一种方法是使用 TransformedTargetRegressor 元估计器将 y_pred 映射到严格正值域。
print("Ridge evaluation:")
score_estimator(ridge_glm, df_test)
Ridge evaluation:
MSE: 0.560
MAE: 0.186
WARNING: Estimator yields invalid, non-positive predictions for 595 samples out of 223745. These predictions are ignored when computing the Poisson deviance.
mean Poisson deviance: 0.597
接下来，我们对目标变量拟合泊松回归器。我们将正则化强度 alpha 设置为大约样本数量的 1e-6 倍（即 1e-12），以模仿岭回归器，其 L2 惩罚项随样本数量的变化而不同。
由于泊松回归器在内部对预期目标值的对数而不是直接预期值进行建模（对数连接函数与恒等连接函数），因此 X 和 y 之间的关系不再是严格线性的。因此，泊松回归器被称为广义线性模型 (GLM)，而不是像岭回归那样的普通线性模型。
from sklearn.linear_model import PoissonRegressor
n_samples = df_train.shape[0]
poisson_glm = Pipeline(
[
("preprocessor", linear_model_preprocessor),
("regressor", PoissonRegressor(alpha=1e-12, solver="newton-cholesky")),
]
)
poisson_glm.fit(
df_train, df_train["Frequency"], regressor__sample_weight=df_train["Exposure"]
)
print("PoissonRegressor evaluation:")
score_estimator(poisson_glm, df_test)
PoissonRegressor evaluation:
MSE: 0.560
MAE: 0.186
mean Poisson deviance: 0.594
用于泊松回归的梯度提升回归树#
最后，我们将考虑一个非线性模型，即梯度提升回归树。基于树的模型不需要对类别数据进行独热编码：相反，我们可以使用 OrdinalEncoder 对每个类别标签进行任意整数编码。通过这种编码，树会将类别特征视为有序特征，这可能并非总是理想的行为。然而，对于足够深的树，这种影响是有限的，它们能够恢复特征的类别性质。OrdinalEncoder 相对于 OneHotEncoder 的主要优点是它会使训练更快。
梯度提升还提供了使用泊松损失（带隐式对数连接函数）而不是默认最小二乘损失来拟合树的可能性。这里我们仅使用泊松损失拟合树，以使本示例简洁。
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.preprocessing import OrdinalEncoder
tree_preprocessor = ColumnTransformer(
[
(
"categorical",
OrdinalEncoder(),
["VehBrand", "VehPower", "VehGas", "Region", "Area"],
),
("numeric", "passthrough", ["VehAge", "DrivAge", "BonusMalus", "Density"]),
],
remainder="drop",
)
poisson_gbrt = Pipeline(
[
("preprocessor", tree_preprocessor),
(
"regressor",
HistGradientBoostingRegressor(loss="poisson", max_leaf_nodes=128),
),
]
)
poisson_gbrt.fit(
df_train, df_train["Frequency"], regressor__sample_weight=df_train["Exposure"]
)
print("Poisson Gradient Boosted Trees evaluation:")
score_estimator(poisson_gbrt, df_test)
Poisson Gradient Boosted Trees evaluation:
MSE: 0.566
MAE: 0.184
mean Poisson deviance: 0.575
与上述泊松 GLM 类似，梯度提升树模型最小化泊松偏差。然而，由于其更高的预测能力，它达到了更低的泊松偏差值。
使用单一训练/测试集拆分评估模型容易出现随机波动。如果计算资源允许，应验证交叉验证的性能指标是否会得出相似的结论。
这些模型之间的定性差异也可以通过比较观测目标值和预测值的直方图来可视化
fig, axes = plt.subplots(nrows=2, ncols=4, figsize=(16, 6), sharey=True)
fig.subplots_adjust(bottom=0.2)
n_bins = 20
for row_idx, label, df in zip(range(2), ["train", "test"], [df_train, df_test]):
df["Frequency"].hist(bins=np.linspace(-1, 30, n_bins), ax=axes[row_idx, 0])
axes[row_idx, 0].set_title("Data")
axes[row_idx, 0].set_yscale("log")
axes[row_idx, 0].set_xlabel("y (observed Frequency)")
axes[row_idx, 0].set_ylim([1e1, 5e5])
axes[row_idx, 0].set_ylabel(label + " samples")
for idx, model in enumerate([ridge_glm, poisson_glm, poisson_gbrt]):
y_pred = model.predict(df)
pd.Series(y_pred).hist(
bins=np.linspace(-1, 4, n_bins), ax=axes[row_idx, idx + 1]
)
axes[row_idx, idx + 1].set(
title=model[-1].__class__.__name__,
yscale="log",
xlabel="y_pred (predicted expected Frequency)",
)
plt.tight_layout()
实验数据中 y 呈现长尾分布。在所有模型中，我们预测的是随机变量的预期频率，因此极端值必然少于该随机变量的观测实现值。这解释了模型预测直方图的众数不一定对应于最小值。此外，Ridge 中使用的正态分布具有恒定方差，而 PoissonRegressor 和 HistGradientBoostingRegressor 中使用的泊松分布，其方差与预测的期望值成正比。
因此，在所考虑的估计器中，PoissonRegressor 和 HistGradientBoostingRegressor 先验地更适合对非负数据的长尾分布进行建模，而 Ridge 模型对目标变量的分布做出了错误的假设。
`HistGradientBoostingRegressor` 估计器具有最大的灵活性，能够预测更高的期望值。
请注意，我们可以为 HistGradientBoostingRegressor 模型使用最小二乘损失。这会像 Ridge 模型一样错误地假设响应变量呈正态分布，并且可能还会导致略微负的预测值。然而，梯度提升树仍然会表现相对良好，特别是由于树的灵活性与大量训练样本相结合，其性能优于 PoissonRegressor。
预测校准评估#
为了确保估计器对不同类型的投保人产生合理的预测，我们可以根据每个模型返回的 y_pred 对测试样本进行分箱。然后对于每个箱，我们比较平均预测 y_pred 与平均观测目标值。
from sklearn.utils import gen_even_slices
def _mean_frequency_by_risk_group(y_true, y_pred, sample_weight=None, n_bins=100):
"""Compare predictions and observations for bins ordered by y_pred.
We order the samples by ``y_pred`` and split it in bins.
In each bin the observed mean is compared with the predicted mean.
Parameters
----------
y_true: array-like of shape (n_samples,)
Ground truth (correct) target values.
y_pred: array-like of shape (n_samples,)
Estimated target values.
sample_weight : array-like of shape (n_samples,)
Sample weights.
n_bins: int
Number of bins to use.
Returns
-------
bin_centers: ndarray of shape (n_bins,)
bin centers
y_true_bin: ndarray of shape (n_bins,)
average y_pred for each bin
y_pred_bin: ndarray of shape (n_bins,)
average y_pred for each bin
"""
idx_sort = np.argsort(y_pred)
bin_centers = np.arange(0, 1, 1 / n_bins) + 0.5 / n_bins
y_pred_bin = np.zeros(n_bins)
y_true_bin = np.zeros(n_bins)
for n, sl in enumerate(gen_even_slices(len(y_true), n_bins)):
weights = sample_weight[idx_sort][sl]
y_pred_bin[n] = np.average(y_pred[idx_sort][sl], weights=weights)
y_true_bin[n] = np.average(y_true[idx_sort][sl], weights=weights)
return bin_centers, y_true_bin, y_pred_bin
print(f"Actual number of claims: {df_test['ClaimNb'].sum()}")
fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(12, 8))
plt.subplots_adjust(wspace=0.3)
for axi, model in zip(ax.ravel(), [ridge_glm, poisson_glm, poisson_gbrt, dummy]):
y_pred = model.predict(df_test)
y_true = df_test["Frequency"].values
exposure = df_test["Exposure"].values
q, y_true_seg, y_pred_seg = _mean_frequency_by_risk_group(
y_true, y_pred, sample_weight=exposure, n_bins=10
)
# Name of the model after the estimator used in the last step of the
# pipeline.
print(f"Predicted number of claims by {model[-1]}: {np.sum(y_pred * exposure):.1f}")
axi.plot(q, y_pred_seg, marker="x", linestyle="--", label="predictions")
axi.plot(q, y_true_seg, marker="o", linestyle="--", label="observations")
axi.set_xlim(0, 1.0)
axi.set_ylim(0, 0.5)
axi.set(
title=model[-1],
xlabel="Fraction of samples sorted by y_pred",
ylabel="Mean Frequency (y_pred)",
)
axi.legend()
plt.tight_layout()
Actual number of claims: 11935
Predicted number of claims by Ridge(alpha=1e-06): 11933.4
Predicted number of claims by PoissonRegressor(alpha=1e-12, solver='newton-cholesky'): 11932.0
Predicted number of claims by HistGradientBoostingRegressor(loss='poisson', max_leaf_nodes=128): 12196.1
Predicted number of claims by DummyRegressor(): 11931.2
哑回归模型预测一个常数频率。该模型不会将相同的并列排名归因于所有样本，但它在全球范围内仍然是良好校准的（用于估计整个群体的平均频率）。
`Ridge` 回归模型可以预测与数据不符的非常低的预期频率。因此，它可能会严重低估某些投保人的风险。
`PoissonRegressor` 和 `HistGradientBoostingRegressor` 在预测目标和观测目标之间显示出更好的一致性，特别是对于低预测目标值。
所有预测的总和也证实了 `Ridge` 模型的校准问题：它低估了测试集中总索赔数量的 3% 以上，而其他三个模型则可以近似恢复测试组合的总索赔数量。
排名能力评估#
对于某些商业应用，我们感兴趣的是模型将最危险的投保人与最安全的投保人进行排名的能力，而与预测的绝对值无关。在这种情况下，模型评估将把问题视为排名问题而不是回归问题。
为了从这个角度比较这 3 个模型，可以绘制索赔的累积比例与风险敞口的累积比例图，其中测试样本按模型预测从最安全到最危险的顺序排列。
此图称为洛伦兹曲线，可用基尼指数进行概括
from sklearn.metrics import auc
def lorenz_curve(y_true, y_pred, exposure):
y_true, y_pred = np.asarray(y_true), np.asarray(y_pred)
exposure = np.asarray(exposure)
# order samples by increasing predicted risk:
ranking = np.argsort(y_pred)
ranked_frequencies = y_true[ranking]
ranked_exposure = exposure[ranking]
cumulated_claims = np.cumsum(ranked_frequencies * ranked_exposure)
cumulated_claims /= cumulated_claims[-1]
cumulated_exposure = np.cumsum(ranked_exposure)
cumulated_exposure /= cumulated_exposure[-1]
return cumulated_exposure, cumulated_claims
fig, ax = plt.subplots(figsize=(8, 8))
for model in [dummy, ridge_glm, poisson_glm, poisson_gbrt]:
y_pred = model.predict(df_test)
cum_exposure, cum_claims = lorenz_curve(
df_test["Frequency"], y_pred, df_test["Exposure"]
)
gini = 1 - 2 * auc(cum_exposure, cum_claims)
label = "{} (Gini: {:.2f})".format(model[-1], gini)
ax.plot(cum_exposure, cum_claims, linestyle="-", label=label)
# Oracle model: y_pred == y_test
cum_exposure, cum_claims = lorenz_curve(
df_test["Frequency"], df_test["Frequency"], df_test["Exposure"]
)
gini = 1 - 2 * auc(cum_exposure, cum_claims)
label = "Oracle (Gini: {:.2f})".format(gini)
ax.plot(cum_exposure, cum_claims, linestyle="-.", color="gray", label=label)
# Random Baseline
ax.plot([0, 1], [0, 1], linestyle="--", color="black", label="Random baseline")
ax.set(
title="Lorenz curves by model",
xlabel="Cumulative proportion of exposure (from safest to riskiest)",
ylabel="Cumulative proportion of claims",
)
ax.legend(loc="upper left")
<matplotlib.legend.Legend object at 0x7facebf6f250>
正如预期的那样，哑回归器无法正确对样本进行排名，因此在此图上表现最差。
基于树的模型在按风险对投保人进行排名方面表现显著更好，而两个线性模型的表现相似。
所有三个模型都明显优于随机，但离完美预测仍有很大距离。
最后一点是意料之中的，因为问题的性质：事故的发生主要受数据集中未捕获的偶然原因主导，并且确实可以被视为纯粹随机的。
线性模型假设输入变量之间没有交互作用，这可能导致欠拟合。插入多项式特征提取器（PolynomialFeatures）确实将它们的判别能力提高了 2 个基尼指数点。特别是，它提高了模型识别前 5% 最危险配置文件的能力。
主要收获#
模型的性能可以通过其产生良好校准的预测和良好排名的能力来评估。
模型校准可以通过绘制按预测风险分箱的测试样本组的平均观测值与平均预测值来评估。
岭回归模型的最小二乘损失（以及恒等连接函数的隐式使用）似乎导致该模型校准不良。特别是，它倾向于低估风险，甚至可能预测无效的负频率。
使用带对数连接的泊松损失可以纠正这些问题，并产生一个良好校准的线性模型。
基尼指数反映了模型对预测进行排名的能力，而与它们的绝对值无关，因此只评估它们的排名能力。
尽管校准有所改进，但两个线性模型的排名能力相当，远低于梯度提升回归树的排名能力。
作为评估指标计算的泊松偏差反映了模型的校准和排名能力。它还对响应变量的期望值和方差之间的理想关系做出了线性假设。为了简洁起见，我们没有检查这个假设是否成立。
传统的回归指标，如均方误差和平均绝对误差，在包含许多零的计数值上很难有意义地解释。
plt.show()
脚本总运行时间： (0 分钟 19.324 秒)
相关示例
