注意
跳转到末尾 下载完整的示例代码。或者通过JupyterLite或Binder在浏览器中运行此示例。
泊松回归和非正态损失#
此示例说明了在法国机动车第三方责任索赔数据集（来自[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，单位为年）内以恒定速率发生的离散事件的数量。
在这里，我们想通过（缩放的）泊松分布对y = ClaimNb / Exposure（条件为X）进行建模，并使用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 中使用的泊松分布，方差与预测的期望值成正比。
因此，在所考虑的估计器中，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% 以上，而其他三个模型能够近似恢复测试组合的索赔总数。
排序能力评估#
对于某些业务应用，我们感兴趣的是模型对投保人风险进行排序的能力，而不管预测的绝对值是多少。在这种情况下，模型评估将问题视为排序问题而不是回归问题。
为了从这个角度比较这三个模型，可以绘制累积索赔比例与累积风险暴露比例的关系图，测试样本按模型预测排序，从最安全到最危险。
该图称为洛伦兹曲线，可以用基尼系数来总结。
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 0x74b4c2e41a00>
正如预期的那样，虚拟回归器无法正确对样本进行排序，因此在此图上的表现最差。
基于树的模型在对投保人进行风险排序方面明显更好，而两个线性模型的表现相似。
所有三个模型都明显优于偶然性，但也远非做出完美的预测。
最后一点是由于问题的性质而预期的：事故的发生主要由数据集列中未捕捉到的偶然原因决定，实际上可以认为是纯粹随机的。
线性模型假设输入变量之间没有相互作用，这很可能导致欠拟合。插入多项式特征提取器（PolynomialFeatures）确实通过 2 个基尼指数点提高了它们的判别能力。特别是它提高了模型识别前 5% 最危险人群的能力。
主要结论#
模型的性能可以通过它们产生良好校准预测和良好排序的能力来评估。
可以通过绘制按预测风险分组的测试样本组的平均观察值与平均预测值来评估模型的校准。
Ridge回归模型的最小二乘损失（以及隐式使用恒等链接函数）似乎导致该模型校准不良。特别是，它往往低估风险，甚至可能预测无效的负频率。使用带有对数链接的泊松损失可以纠正这些问题，并导致一个良好校准的线性模型。
基尼指数反映了模型对预测进行排序的能力，而不管其绝对值如何，因此只评估其排序能力。
尽管校准有所改进，但两个线性模型的排序能力都相当，并且远低于梯度提升回归树的排序能力。
泊松偏差作为评估指标，反映了模型的校准能力和排序能力。它还对响应变量的期望值和方差之间的理想关系作出了线性假设。为简洁起见，我们没有检查该假设是否成立。
传统的回归指标，如均方误差和平均绝对误差，很难对包含许多零的计数值进行有意义的解释。
plt.show()
脚本总运行时间：（0 分 22.346 秒）
相关示例
