<jupyter_start><jupyter_text>Poisson regression for predicting insurance claim rateshttps://scikit-learn.org/stable/auto_examples/linear_model/plot_poisson_regression_non_normal_loss.html?highlight=poisson%20regression<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
!pip install -U scikit-learn
import sklearn
print(sklearn.__version__)
from sklearn.linear_model import PoissonRegressor<jupyter_output>0.24.2<jupyter_text>Data<jupyter_code>from sklearn.datasets import fetch_openml


df = fetch_openml(data_id=41214, as_frame=True).frame
df
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, ax = plt.subplots()
ax.set_title("Frequency (number of claims per year)")
_ = df["Frequency"].hist(bins=30, log=True, ax=ax)
print(df["Frequency"][:-20])

from sklearn.model_selection import train_test_split

df_train, df_test = train_test_split(df, test_size=0.33, random_state=0)
fig, ax = plt.subplots()
n_bins = 20
_ = df_test["Frequency"].hist(bins=np.linspace(-1, 30, n_bins), ax=ax)
ax.set_yscale('log')<jupyter_output><empty_output><jupyter_text>Feature engineering<jupyter_code>from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, OneHotEncoder
from sklearn.preprocessing import StandardScaler, KBinsDiscretizer
from sklearn.compose import ColumnTransformer


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),
            ["VehAge", "DrivAge"]),
        ("log_scaled_numeric", log_scale_transformer,
            ["Density"]),
        ("onehot_categorical", OneHotEncoder(),
            ["VehBrand", "VehPower", "VehGas", "Region", "Area"]),
    ],
    remainder="drop",
)<jupyter_output><empty_output><jupyter_text>Evaluation metrics<jupyter_code>from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_poisson_deviance


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(f"WARNING: Estimator yields invalid, non-positive predictions "
              f" for {n_masked} samples out of {n_samples}. These predictions "
              f"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]))
def plot_predictions(models):
  nmodels = len(models)
  height = 5
  width = 5*(nmodels+1)
  fig, axes = plt.subplots(nrows=1, ncols=nmodels+1, figsize=(width,height), sharey=True)
  n_bins = 20
  df = df_test.copy()
  ax = axes[0]
  df["Frequency"].hist(bins=np.linspace(-1, 30, n_bins), ax=ax)
  ax.set_title("Data")
  ax.set_yscale('log')
  ax.set_xlabel("y (observed Frequency)")
  #ax.set_ylim([1e1, 5e5])
  ax.set_ylabel("#samples")

  for idx, model in enumerate(models):
    ax = axes[idx+1]
    y_pred = model.predict(df)
    pd.Series(y_pred).hist(bins=np.linspace(-1, 4, n_bins),ax=ax)
    ax.set(
        title=model[-1].__class__.__name__,
        yscale='log',
        xlabel="y_pred (predicted expected Frequency)"
    )

  plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Dummy modelJust predicts overall mean.<jupyter_code>from sklearn.dummy import DummyRegressor
from sklearn.pipeline import Pipeline

dummy = Pipeline([
    ("preprocessor", linear_model_preprocessor),
    ("regressor", DummyRegressor(strategy='mean')),
]).fit(df_train, df_train["Frequency"])
y_pred = dummy.predict(df_test)
print(y_pred[0:5])<jupyter_output>[0.26520328 0.26520328 0.26520328 0.26520328 0.26520328]<jupyter_text>We need to weight the examples by exposure, for reasons explained here:https://github.com/scikit-learn/scikit-learn/issues/18059<jupyter_code># weighted version
from sklearn.dummy import DummyRegressor
from sklearn.pipeline import Pipeline

dummy = Pipeline([
    ("preprocessor", linear_model_preprocessor),
    ("regressor", DummyRegressor(strategy='mean')),
]).fit(df_train, df_train["Frequency"],
       regressor__sample_weight=df_train["Exposure"])
y_pred = dummy.predict(df_test)
print(y_pred[0:5])
print("Constant mean frequency evaluation:")
score_estimator(dummy, df_test)<jupyter_output>Constant mean frequency evaluation:
MSE: 0.564
MAE: 0.189
mean Poisson deviance: 0.625<jupyter_text>Linear regressionLinear regression is an okay baseline, but not the most appropriate model for count data... We see the L2 regularization to a small value, since the training set is large.<jupyter_code>from sklearn.linear_model import Ridge


ridge_glm = Pipeline([
    ("preprocessor", linear_model_preprocessor),
    ("regressor", Ridge(alpha=1e-6)),
]).fit(df_train, df_train["Frequency"],
       regressor__sample_weight=df_train["Exposure"])
print("Ridge evaluation:")
score_estimator(ridge_glm, df_test)<jupyter_output>Ridge evaluation:
MSE: 0.560
MAE: 0.177
WARNING: Estimator yields invalid, non-positive predictions  for 1315 samples out of 223745. These predictions are ignored when computing the Poisson deviance.
mean Poisson deviance: 0.601<jupyter_text>Poisson linear regressionWe set the L2 regularizer to the same value as in ridge regression, but divided by the number of training samples, since poisson regression penalizes the average log likelihood.<jupyter_code>from sklearn.linear_model import PoissonRegressor

n_samples = df_train.shape[0]

poisson_glm = Pipeline([
    ("preprocessor", linear_model_preprocessor),
    ("regressor", PoissonRegressor(alpha=1e-12, max_iter=300))
])
poisson_glm.fit(df_train, df_train["Frequency"],
                regressor__sample_weight=df_train["Exposure"])


print("PoissonRegressor evaluation:")
score_estimator(poisson_glm, df_test)<jupyter_output>PoissonRegressor evaluation:
MSE: 0.560
MAE: 0.186
mean Poisson deviance: 0.594<jupyter_text>Comparison<jupyter_code>plot_predictions([ridge_glm, poisson_glm])
plot_predictions([dummy, ridge_glm, poisson_glm])
plt.savefig('poisson_regr_insurance_pred.pdf')
for model in [dummy, ridge_glm, poisson_glm]:
  print(model[-1].__class__.__name__)
  score_estimator(model, df_test)<jupyter_output>DummyRegressor
MSE: 0.564
MAE: 0.189
mean Poisson deviance: 0.625
Ridge
MSE: 0.560
MAE: 0.177
WARNING: Estimator yields invalid, non-positive predictions  for 1315 samples out of 223745. These predictions are ignored when computing the Poisson deviance.
mean Poisson deviance: 0.601
PoissonRegressor
MSE: 0.560
MAE: 0.186
mean Poisson deviance: 0.594<jupyter_text>Calibration plot<jupyter_code>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


def plot_calibration(models):
  print(f"Actual number of claims: {df_test['ClaimNb'].sum()}")
  nmodels = len(models)
  height = 5
  width = 5*(nmodels+1)
  fig, ax = plt.subplots(nrows=1, ncols=nmodels, figsize=(width, height))
  plt.subplots_adjust(wspace=0.3)
  for axi, model in zip(ax.ravel(), models):
      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]}: "
            f"{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()

plot_calibration([dummy, ridge_glm, poisson_glm])
plt.savefig('poisson_regr_insurance.pdf')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11792041823990705, 'fraction_numerical': 0.022848291218898246, 'mean_word_length': 4.198792148968294, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 36, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '190851', 'license': 'MIT', 'n_tokens_mistral': 3754, 'n_tokens_neox': 3448, 'n_words': 763, 'path': 'notebooks/poisson_regression_insurance.ipynb', 'repo_name': 'arpitvaghela/probml-notebooks', 'size': '10329'}