泊松回归和非正态损失#
本示例说明了对来自 法国机动车第三方责任索赔数据集 的对数线性泊松回归的使用，该数据集来自 [1]，并将其与使用通常的最小二乘误差拟合的线性模型和使用泊松损失（和对数链接）拟合的非线性 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, 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 中使用的泊松分布，方差与预测的期望值成正比。
因此，在所考虑的估计器中，与对目标变量的分布做出错误假设的 Ridge 模型相比，PoissonRegressor 和 HistGradientBoostingRegressor 先验地更适合于对非负数据的长尾分布进行建模。
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 0x7f4dfc3be520>
正如预期的那样，虚拟回归器无法正确地对样本进行排名，因此在此图上的表现最差。
基于树的模型在按风险对投保人进行排名方面明显更好，而两个线性模型的表现相似。
所有三个模型都明显优于随机模型，但也远未做出完美的预测。
最后一点是预料之中的，因为问题的性质：事故的发生主要由数据集中各列未捕获的偶然原因主导，并且实际上可以被视为纯粹的随机事件。
线性模型假设输入变量之间没有交互作用，这可能会导致欠拟合。插入多项式特征提取器 (PolynomialFeatures) 确实可以通过 2 个基尼指数点来提高它们的判别能力。特别是，它提高了模型识别前 5% 最危险配置文件的能力。
主要结论#
模型的性能可以通过其产生良好校准的预测和良好排名的能力来评估。
可以通过绘制按预测风险分类的测试样本组的平均观察值与平均预测值的关系图来评估模型的校准。
岭回归模型的最小二乘损失（以及隐式使用恒等链接函数）似乎导致该模型校准不良。特别是，它往往会低估风险，甚至会预测无效的负频率。
使用带有对数链接的泊松损失可以纠正这些问题，并导致一个校准良好的线性模型。
基尼系数反映了模型对预测进行排名的能力，而不管其绝对值如何，因此只能评估其排名能力。
尽管校准有所改进，但两个线性模型的排名能力相当，远低于梯度提升回归树的排名能力。
作为评估指标计算的泊松偏差反映了模型的校准和排名能力。它还对响应变量的期望值和方差之间的理想关系做出了线性假设。为了简洁起见，我们没有检查这个假设是否成立。
传统的回归指标，如均方误差和平均绝对误差，很难对具有许多零的计数值进行有意义的解释。
plt.show()
脚本总运行时间：（0 分钟 27.378 秒）
相关示例
