20 Défis Actuariels

Code de démarrage

# Chargement
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

df = pd.read_csv('freMTPLfreq.csv')
print(df.shape, df.dtypes)
print(df['ClaimNb'].value_counts(normalize=True))

# Fréquence moyenne par segment
freq_mean = df.groupby('VehBrand').apply(
    lambda x: x['ClaimNb'].sum() / x['Exposure'].sum()
).sort_values(ascending=False)
print(freq_mean)

# Distribution de l'exposition
sns.histplot(df['Exposure'], bins=50)
plt.title("Distribution de l'exposition")
plt.show()

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pandas as pd
import matplotlib.pyplot as plt

df = pd.read_csv('freMTPLfreq.csv')

# Relativité brute = fréquence observée / fréquence globale
freq_globale = df['ClaimNb'].sum() / df['Exposure'].sum()

for var in ['VehPower', 'VehAge', 'DrivAge', 'BonusMalus']:
    rel = df.groupby(var).apply(
        lambda x: (x['ClaimNb'].sum() / x['Exposure'].sum()) / freq_globale
    )
    rel.plot(kind='bar', title=f'Relativités — {var}', figsize=(10, 4))
    plt.axhline(1, color='red', linestyle='--')
    plt.tight_layout()
    plt.show()

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

df = pd.read_csv('freMTPLfreq.csv')
features = ['VehPower', 'VehAge', 'DrivAge', 'BonusMalus', 'Density']
X = StandardScaler().fit_transform(df[features].fillna(0))

kmeans = KMeans(n_clusters=5, random_state=42, n_init=10)
df['cluster'] = kmeans.fit_predict(X)

# Visualisation PCA 2D
pca = PCA(n_components=2)
X2 = pca.fit_transform(X)
plt.scatter(X2[:, 0], X2[:, 1], c=df['cluster'], cmap='tab10', alpha=0.3, s=1)
plt.title('K-Means — visualisation PCA')
plt.show()

# Fréquence par cluster
print(df.groupby('cluster').apply(lambda x: x['ClaimNb'].sum() / x['Exposure'].sum()))

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pandas as pd
import statsmodels.formula.api as smf
import numpy as np

df = pd.read_csv('freMTPLfreq.csv')
df['ClaimNb'] = df['ClaimNb'].clip(upper=4)
df['Exposure'] = df['Exposure'].clip(upper=1)

# GLM Poisson avec offset
glm = smf.glm(
    formula='ClaimNb ~ VehPower + C(VehBrand) + DrivAge + BonusMalus + np.log(Density)',
    data=df,
    family=sm.families.Poisson(link=sm.families.links.log()),
    offset=np.log(df['Exposure'])
).fit()

print(glm.summary())

# Relativités (exp des coefficients)
print(np.exp(glm.params))

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

freq = pd.read_csv('freMTPLfreq.csv')
sev  = pd.read_csv('freMTPLsev.csv')

# Modèle fréquence (Poisson)
freq['ClaimNb'] = freq['ClaimNb'].clip(upper=4)
glm_freq = smf.glm('ClaimNb ~ BonusMalus + DrivAge + VehAge',
    data=freq, family=sm.families.Poisson(),
    offset=np.log(freq['Exposure'])).fit()

# Modèle sévérité (Gamma) — sur sinistres seulement
glm_sev = smf.glm('ClaimAmount ~ BonusMalus + DrivAge',
    data=sev[sev['ClaimAmount'] > 0],
    family=sm.families.Gamma(link=sm.families.links.log())).fit()

# Prime pure = fréquence prédite × sévérité prédite
freq['pred_freq'] = glm_freq.predict(freq)
freq['pred_sev']  = glm_sev.predict(freq)
freq['pure_premium'] = freq['pred_freq'] * freq['pred_sev']
print(freq['pure_premium'].describe())

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

def gini_coefficient(y_true, y_pred, exposure):
    """Gini normalisé — version actuarielle."""
    df = pd.DataFrame({'obs': y_true, 'pred': y_pred, 'exp': exposure})
    df = df.sort_values('pred')
    df['cum_exp'] = df['exp'].cumsum() / df['exp'].sum()
    df['cum_obs'] = (df['obs']).cumsum() / df['obs'].sum()
    gini = 1 - 2 * np.trapz(df['cum_obs'], df['cum_exp'])
    return gini

# Après avoir ajusté votre GLM :
# gini = gini_coefficient(y_test, y_pred, exposure_test)
# print(f"Gini : {gini:.4f}")

# Courbe de Lorenz
plt.plot([0, 1], [0, 1], 'k--', label='Aléatoire')
# plt.plot(cum_exp, cum_obs, label=f'Modèle (Gini={gini:.3f})')
plt.xlabel('Exposition cumulée'); plt.ylabel('Sinistres cumulés')
plt.title('Courbe de Lorenz'); plt.legend(); plt.show()

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pandas as pd, numpy as np
import xgboost as xgb
import lightgbm as lgb

df = pd.read_csv('freMTPLfreq.csv')
df['ClaimNb'] = df['ClaimNb'].clip(upper=4)
df['Exposure'] = df['Exposure'].clip(upper=1)

features = ['VehPower','VehAge','DrivAge','BonusMalus','Density']
X = df[features]; y = df['ClaimNb']; w = df['Exposure']

# XGBoost Poisson
xgb_model = xgb.XGBRegressor(
    objective='count:poisson', n_estimators=500,
    learning_rate=0.05, max_depth=4,
    monotone_constraints={'BonusMalus': 1}
)
xgb_model.fit(X, y, sample_weight=w)

# LightGBM Poisson
lgb_model = lgb.LGBMRegressor(
    objective='poisson', n_estimators=500,
    learning_rate=0.05, max_depth=4,
    monotone_constraints=[0, 0, 0, 1, 0]
)
lgb_model.fit(X, y, sample_weight=w)

def poisson_deviance(y_true, y_pred, exposure):
    y_pred = y_pred / exposure
    return 200 * np.mean(y_pred - y_true + y_true * np.log((y_true + 1e-10) / y_pred))

print("XGBoost:", poisson_deviance(y, xgb_model.predict(X), w))
print("LightGBM:", poisson_deviance(y, lgb_model.predict(X), w))

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import shap, xgboost as xgb, pandas as pd

df = pd.read_csv('freMTPLfreq.csv')
features = ['VehPower','VehAge','DrivAge','BonusMalus','Density']
X = df[features]

# Après entraînement du modèle XGBoost :
explainer = shap.TreeExplainer(xgb_model)
shap_values = explainer.shap_values(X)

# Beeswarm — importance globale
shap.summary_plot(shap_values, X)

# Dependence plot — effet BonusMalus
shap.dependence_plot('BonusMalus', shap_values, X, interaction_index='DrivAge')

# Explication d'une police individuelle
shap.waterfall_plot(shap.Explanation(
    values=shap_values[0], base_values=explainer.expected_value, data=X.iloc[0]
))

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import tensorflow as tf
from tensorflow import keras
import numpy as np

# Architecture CANN : GLM output + correction neuronale
def build_cann(n_features):
    inp = keras.Input(shape=(n_features,))

    # Branche GLM (log-linéaire)
    glm_out = keras.layers.Dense(1, activation='linear', name='glm')(inp)

    # Branche neuronale (correction)
    x = keras.layers.Dense(20, activation='tanh')(inp)
    x = keras.layers.Dense(15, activation='tanh')(x)
    nn_out = keras.layers.Dense(1, activation='linear', name='nn')(x)

    # Skip connection : GLM + correction
    combined = keras.layers.Add()([glm_out, nn_out])
    output = keras.layers.Activation('exponential')(combined)

    return keras.Model(inp, output)

model = build_cann(5)
model.compile(optimizer='adam', loss='poisson')
# model.fit(X_train, y_train, sample_weight=exposure_train, epochs=50, batch_size=1024)

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pandas as pd
import statsmodels.api as sm
from statsmodels.discrete.count_model import ZeroInflatedPoisson

df = pd.read_csv('freMTPLfreq.csv')
print(f"Proportion zéros : {(df['ClaimNb']==0).mean():.1%}")

X = sm.add_constant(df[['BonusMalus', 'DrivAge', 'VehAge']])

# Zero-Inflated Poisson
zip_model = ZeroInflatedPoisson(df['ClaimNb'], X, exog_infl=X).fit()
print(zip_model.summary())

# Negative Binomial (surdispersion)
nb_model = sm.NegativeBinomial(df['ClaimNb'], X,
    offset=pd.np.log(df['Exposure'])).fit()
print(nb_model.summary())

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pandas as pd
import plotly.express as px

df = pd.read_csv('freMTPLfreq.csv')

# Fréquence par région
freq_region = df.groupby('Region').apply(
    lambda x: x['ClaimNb'].sum() / x['Exposure'].sum()
).reset_index(name='freq')
freq_region.columns = ['region', 'frequence']

# Carte choroplèthe avec Plotly (GeoJSON régions françaises)
fig = px.choropleth(
    freq_region, locations='region',
    color='frequence',
    color_continuous_scale='RdYlGn_r',
    title='Fréquence sinistres par région française'
)
fig.show()
# Note: pour les contours précis, utiliser un GeoJSON des régions FR

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pandas as pd, numpy as np
import matplotlib.pyplot as plt
from pygam import PoissonGAM, s, f

df = pd.read_csv('freMTPLfreq.csv')
df['ClaimNb'] = df['ClaimNb'].clip(upper=4)

X = df[['BonusMalus', 'DrivAge', 'VehAge', 'Density', 'VehPower']].values
y = df['ClaimNb'].values

gam = PoissonGAM(
    s(0) + s(1) + s(2) + s(3) + f(4)  # splines + facteur
).fit(X, y)

# Visualiser les effets partiels
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
titles = ['BonusMalus', 'DrivAge', 'VehAge']
for i, (ax, title) in enumerate(zip(axes, titles)):
    XX = gam.generate_X_grid(term=i)
    ax.plot(XX[:, i], gam.partial_dependence(term=i, X=XX))
    ax.set_title(title)
plt.tight_layout(); plt.show()

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pandas as pd
from sklearn.ensemble import IsolationForest
from sklearn.neighbors import LocalOutlierFactor
from sklearn.preprocessing import StandardScaler

df = pd.read_csv('freMTPLfreq.csv')
features = ['VehPower','VehAge','DrivAge','BonusMalus','Density','Exposure']
X = StandardScaler().fit_transform(df[features].fillna(0))

# Isolation Forest
iso = IsolationForest(contamination=0.01, random_state=42)
df['anomaly_iso'] = iso.fit_predict(X)

# Local Outlier Factor
lof = LocalOutlierFactor(n_neighbors=20, contamination=0.01)
df['anomaly_lof'] = lof.fit_predict(X)

# Polices identifiées comme anomalies par les deux
anomalies = df[(df['anomaly_iso'] == -1) & (df['anomaly_lof'] == -1)]
print(f"{len(anomalies)} polices atypiques ({len(anomalies)/len(df):.2%})")
print(anomalies[features].describe())

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pandas as pd
from lifelines import KaplanMeierFitter, CoxPHFitter
import matplotlib.pyplot as plt

df = pd.read_csv('freMTPLfreq.csv')
# duration = Exposure (temps d'observation)
# event = 1 si au moins 1 sinistre, 0 sinon (censuré)
df['event'] = (df['ClaimNb'] > 0).astype(int)
df['duration'] = df['Exposure']

# Kaplan-Meier par groupe BonusMalus
kmf = KaplanMeierFitter()
for group in ['low', 'high']:
    mask = df['BonusMalus'] > 100 if group == 'high' else df['BonusMalus'] <= 100
    kmf.fit(df[mask]['duration'], df[mask]['event'], label=f'BonusMalus {group}')
    kmf.plot_survival_function()
plt.title('Courbe de survie — temps avant premier sinistre')
plt.show()

# Modèle de Cox
cph = CoxPHFitter()
cph.fit(df[['duration','event','BonusMalus','DrivAge','VehAge']],
        duration_col='duration', event_col='event')
cph.print_summary()

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pandas as pd, numpy as np
from fairlearn.reductions import ExponentiatedGradient, BoundedGroupLoss
from fairlearn.metrics import MetricFrame
import lightgbm as lgb

df = pd.read_csv('freMTPLfreq.csv')
features = ['VehPower','VehAge','DrivAge','BonusMalus','Density']
X = df[features]; y = df['ClaimNb']
sensitive = df['Region']  # variable sensible

# Modèle de base
base_model = lgb.LGBMRegressor(objective='poisson')

# Réduction d'inéquité via ExponentiatedGradient
mitigator = ExponentiatedGradient(
    base_model,
    constraints=BoundedGroupLoss(loss=None, upper_bound=0.1)
)
mitigator.fit(X, y, sensitive_features=sensitive)

# Analyse par groupe
preds = mitigator.predict(X)
mf = MetricFrame(
    metrics={'mean_pred': np.mean},
    y_true=y, y_pred=preds,
    sensitive_features=sensitive
)
print(mf.by_group)

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import optuna, lightgbm as lgb, numpy as np, pandas as pd
from sklearn.model_selection import KFold

df = pd.read_csv('freMTPLfreq.csv')
features = ['VehPower','VehAge','DrivAge','BonusMalus','Density']
X = df[features].values; y = df['ClaimNb'].values; w = df['Exposure'].values

def poisson_deviance(y_true, y_pred, w):
    mu = y_pred / w
    return 2 * np.mean(w * (mu - y_true + y_true * np.log((y_true + 1e-9) / mu)))

def objective(trial):
    params = {
        'objective': 'poisson', 'n_estimators': 500, 'verbose': -1,
        'learning_rate': trial.suggest_float('lr', 0.01, 0.3, log=True),
        'max_depth': trial.suggest_int('depth', 3, 8),
        'num_leaves': trial.suggest_int('leaves', 15, 127),
        'min_child_samples': trial.suggest_int('min_child', 20, 200),
    }
    scores = []
    for tr, va in KFold(3).split(X):
        m = lgb.LGBMRegressor(**params).fit(X[tr], y[tr], sample_weight=w[tr])
        scores.append(poisson_deviance(y[va], m.predict(X[va]), w[va]))
    return np.mean(scores)

study = optuna.create_study(direction='minimize')
study.optimize(objective, n_trials=50)
print("Meilleurs params:", study.best_params)

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pandas as pd, numpy as np

df = pd.read_csv('freMTPLfreq.csv')

# Statistiques par région
reg = df.groupby('Region').agg(
    n=('ClaimNb', 'sum'),
    exposure=('Exposure', 'sum')
).assign(freq=lambda x: x['n'] / x['exposure'])

# Fréquence globale
mu = reg['n'].sum() / reg['exposure'].sum()

# Variance entre régions (between) et intra (within)
grand_mean = mu
between_var = np.var(reg['freq'], ddof=1)
within_var = (reg['n'] / reg['exposure']**2).mean()  # approximation

# Paramètre k de Bühlmann-Straub
k = within_var / between_var

# Facteur de crédibilité Z = m / (m + k)
reg['Z'] = reg['exposure'] / (reg['exposure'] + k)
reg['freq_credibility'] = reg['Z'] * reg['freq'] + (1 - reg['Z']) * mu

print(reg[['freq', 'freq_credibility', 'Z']].sort_values('Z'))

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import numpy as np, pandas as pd
from scipy import stats
from scipy.special import logsumexp

df = pd.read_csv('freMTPLfreq.csv')
y = df['ClaimNb'].values; w = df['Exposure'].values

# Poisson Mixture à 2 composantes — algorithme EM
K = 2
lambdas = np.array([0.05, 0.30])  # init
pis = np.array([0.8, 0.2])

for iteration in range(100):
    # E-step : responsabilités
    log_resp = np.array([
        np.log(pis[k]) + stats.poisson.logpmf(y, lambdas[k] * w)
        for k in range(K)
    ])
    log_norm = logsumexp(log_resp, axis=0)
    resp = np.exp(log_resp - log_norm)

    # M-step
    Nk = resp.sum(axis=1)
    pis = Nk / Nk.sum()
    lambdas = np.array([
        (resp[k] * y).sum() / (resp[k] * w).sum()
        for k in range(K)
    ])

print(f"Composante 1 : λ={lambdas[0]:.4f}, poids={pis[0]:.3f}")
print(f"Composante 2 : λ={lambdas[1]:.4f}, poids={pis[1]:.3f}")

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

# 1. Tracking MLflow
import mlflow, mlflow.sklearn
import lightgbm as lgb

with mlflow.start_run():
    model = lgb.LGBMRegressor(objective='poisson', n_estimators=500)
    model.fit(X_train, y_train, sample_weight=w_train)

    mlflow.log_param("n_estimators", 500)
    mlflow.log_metric("poisson_deviance", score)
    mlflow.sklearn.log_model(model, "model")

# 2. API FastAPI pour scoring temps réel
from fastapi import FastAPI
from pydantic import BaseModel
import numpy as np

app = FastAPI()

class Police(BaseModel):
    VehPower: int; VehAge: int; DrivAge: int
    BonusMalus: int; Density: float; Exposure: float

@app.post("/score")
def score(police: Police):
    X = np.array([[police.VehPower, police.VehAge, police.DrivAge,
                   police.BonusMalus, police.Density]])
    freq = model.predict(X)[0]
    return {"frequence": round(float(freq), 6),
            "prime_annuelle": round(float(freq * police.Exposure * 1500), 2)}

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.

Code de démarrage

import pymc as pm, numpy as np, pandas as pd, arviz as az

df = pd.read_csv('freMTPLfreq.csv').sample(5000, random_state=42)
y = df['ClaimNb'].values
exposure = df['Exposure'].values
bonus = (df['BonusMalus'] - df['BonusMalus'].mean()) / df['BonusMalus'].std()
age   = (df['DrivAge']   - df['DrivAge'].mean())   / df['DrivAge'].std()

with pm.Model() as model:
    # Priors
    alpha = pm.Normal('alpha', mu=-3, sigma=1)
    b_bonus = pm.Normal('b_bonus', mu=0, sigma=0.5)
    b_age   = pm.Normal('b_age',   mu=0, sigma=0.5)

    # Espérance
    mu = pm.math.exp(alpha + b_bonus * bonus + b_age * age) * exposure

    # Vraisemblance Poisson
    obs = pm.Poisson('obs', mu=mu, observed=y)

    # MCMC
    trace = pm.sample(1000, tune=500, target_accept=0.9, return_inferencedata=True)

az.plot_posterior(trace, var_names=['b_bonus', 'b_age'])
az.summary(trace)

💡 Ce code est un point de départ — à toi de l'adapter, compléter et améliorer.