Au cours d'une période d'exercice, le montant total des indemnisations est donné par
$$ X = \sum_{i = 1}^N U_i, $$ou $N$ est une variable aléatoire de comptage et $(U_i)_{i\geq0}$ est une suite iid de variables aléatoires positives et indépendantes de $N$. Nous disposons d'un historique de données de $t$ périodes d'exercice, nous n'avons cependant pas accès aux données complètes
$$ n_s, \mathbf{u}_s = (u_{s,1},\ldots, u_{s,n_s}),\text{ }s=1,\ldots,t, $$mais seulement à la somme des prestations
$$ x_s = \sum_{i = 1}^{n_s}u_{s,i},\text{ }s=1,\ldots, t. $$qui sont des réalisations de la variable aléatoire $X$. On fait l'hypothèse d'un modèle paramétrique sur la distribution de la fréquence des sinistres $p(n|\theta_{\text{freq}})$ et sur les montants individuels des sinistres $f(u|\theta_{\text{sev}})$.
L'objectif est simple, proposer des estimateurs de $\theta_{\text{freq}}$ et $\theta_{\text{sev}}$ à partir des données $x_1,\ldots, x_t$.
Le modèle sous jacent est le suivant
$$N\sim \text{Geom}(q),\text{ avec }p(n|q) = q^n(1-q),\text{ }n \geq0\text{ et }U\sim \text{Exp}(\delta)\text{, avec }f_U(u|\delta) = \frac{e^{-u/\delta}}{\delta},\text{ }u>0.$$import numpy as np
import pandas as pd
import scipy as sc
import seaborn as sns
from scipy.optimize import minimize
import matplotlib.pyplot as plt
Ecrire une fonction sim_X pour simuler des réalisation de la variable $X$.
t, q, δ = 10, 0.3, 5
def sim_X(t, q, δ):
n_t = np.random.negative_binomial(1, 1-q, size=t)
u_t = [np.random.exponential(scale = δ, size = n) for n in n_t]
x_t = [np.sum(u) for u in u_t]
return(np.array(x_t))
sim_X(t, q, δ)
array([ 0. , 4.09548659, 12.38933536, 0. , 2.47667854,
0. , 0. , 0. , 0. , 0. ])
La distribution de $X$ est mixte au sens ou
$$ \text{d}\mathbb{P}_X(x) = \mathbb{P}(N = 0)\text{d}\delta_0(x) + \sum_{n = 1}^{+\infty}f_U^{\ast n}(x)\mathbb{P}(N = n)\text{d}\lambda(x)= \mathbb{P}(N = 0)\text{d}\delta_0(x) + f^+_X(x)\text{d}\lambda(x) $$Dans le cadre de notre modèle, nous avons
$$ \mathbb{P}(N = 0) = 1-q\text{, et }f^+_X(x) =\text{ ? } $$On ré-arrange l'échantillon avec
$$ \mathbf{x} = (0,\ldots, 0, x_1^+,\ldots, x_{t-t_0}^+), $$où $t_0$ désigne le nombre de $0$ et $x_1^+,\ldots, x_{t-t_0}^+$ les cumuls d'indemnisations des période d'exercice comportant au moins un sinistre.
Comment définir la vraisemblance de l'échantillon $\mathbf{x}$?
$$ L(\mathbf{x}|\theta) = \text{ ? } $$En générale on cherche à maximiser la logvraisemblance notée $$ l(\mathbf{x}|\theta) = \log L(\mathbf{x}|\theta) = ? $$
L'estimateur du maximum de vraisemblance est obtenu en résolvant le sytème
$$ \begin{cases} \frac{\partial l}{\partial q}=0\\ \frac{\partial l}{\partial \delta}=0 \end{cases} $$Ecrire une fonction mle permettant de trouver l'estimateur du maximum de vraisemblance
X = sim_X(1000, q, δ)
def mle(X):
t, t0 = len(X), np.sum(X ==0)
q_mle = (t-t0)/t
δ_mle = t0 / t * np.sum(X) /(t-t0)
return(q_mle, δ_mle)
Faire la même chose en utilisant la fonction minimize de la librairie scipy.optimize minimize function (voir aussi les exemples)
En inférence Bayésienne, le paramètre est supposé aléatoire $\theta \sim \pi(\theta)$ où $\pi(\cdot)$ est la loi a priori qui résume l'information sur le paramètre avant d'avoir vu les données. L'inférence Bayésienne repose sur la mise à jour de l'information a priori par la vraisemblance des données pour obtenir la loi a posteriori
$$ \pi(\theta|\mathbf{x}) = \frac{L(\mathbf{x}|\theta)\pi(\theta)}{\int_{\Theta}L(\mathbf{x}|\theta)\pi(\theta)\text{d}\theta} $$L'intégrale est souvent difficile à évaluer d'où le recours à un algorithme de simulation pour générer des valeurs de paramètres suivant la loi à posteriori.
Nous utiliserons un algorithme MCMC (Markov Chain Monte Carlo) basé sur des pas de Metropolis-Hasting, voici le pseudo code
Pour $i=1,...,I$
Le résultat est un échantillon $\theta_1,\ldots, \theta_I$ distribué suivant la loi a posteriori. Il est nécessaire de choisir judicieusement le paramètre $\Sigma$ de la perturbation pour obtenir un algorithme efficace.
Supposons des lois a priori uniformes et indépendantes
$$ \pi(\theta) = \pi(q)\cdot\pi(\delta) =\frac{1}{b_q- a_q}\mathbb{I}_{[a_q,b_q]}(q)\times \frac{1}{b_\delta- a_\delta}\mathbb{I}_{[a_\delta,b_\delta]}(\delta) $$Ecrire une fonction logp_geom_exp qui prend en argument les données et retourne une fonction qui évalue la vraisemblance pour un jeu de paramètres
t, q, δ = 100, 0.1, 5
X = sim_X(t, q, δ)
def logp_geom_exp(X):
def logp(parms):
q, δ = parms
if np.logical_and(np.all(parms > 0), q < 1 ) :
t0, t = np.sum(X == 0), len(X)
return(t * np.log(1-q) + (t-t0) * np.log(1-q) + (t - t0) * np.log(q) -
(t - t0) * np.log(δ) - (1 - q) / δ * np.sum(X))
else:
return(-np.inf)
return(logp)
logp = logp_geom_exp(X)
logp(np.array([q, δ]))
-46.814485463324615
Ecrire une fonction sim_post_sample_MH pour simuler depuis la loi à posteriori, le résultat doit être un data frame contenant les valeurs simulés pour chacun des paramètres $q$ et $\delta$.
parms_name = ['q', 'δ']
a_q, b_q, a_δ, b_δ = 0, 1, 0, 100
Σ, n_step = np.diag(np.array([0.01, 1])), 10000
def sim_post_sample_MH(X, parms_name, a_q, b_q, a_δ, b_δ, Σ, n_step):
logp = logp_geom_exp(X)
trace = []
accepted = []
q0, δ0 = np.random.uniform(low = a_q, high = b_q, size = 1), np.random.uniform(low = a_δ, high = b_δ, size = 1)
trace.append(np.array([q0[0], δ0[0]])), accepted.append(True)
for j in range(n_step):
ϵ = np.random.multivariate_normal(mean = np.zeros(2), cov = Σ , size = 1)
parms_old, parms_new = trace[-1], trace[-1] + ϵ[0]
u = np.random.uniform(0, 1, size = 1)
if logp(parms_new) - logp(parms_old) > np.log(u):
trace.append(parms_new), accepted.append(True)
else:
trace.append(parms_old), accepted. append(False)
trace, accepted
trace_df = pd.DataFrame(trace)
trace_df.columns = parms_name
trace_df['accepted'] = np.array(accepted)
return(trace_df)
trace = sim_post_sample_MH(X, parms_name, a_q, b_q, a_δ, b_δ, Σ, n_step)
np.mean(trace)
q 0.072419 δ 10.380515 accepted 0.265873 dtype: float64
fig, axs = plt.subplots(1,2)
for j in range(len(parms_name)):
axs[j].plot(trace[parms_name[j]].values)
axs[j].set_title(parms_name[j])
sns.despine()
fig.tight_layout()
trace['q'].values
array([0.1540923 , 0.14132258, 0.14132258, ..., 0.12513508, 0.1114039 ,
0.1114039 ])
sns.pairplot(trace[parms_name].loc[int(n_step/2):n_step])
sns.despine()
Lorsque la dimension de l'espace des paramètres est plus grand que $1$ (deux dans notre cas), il est possible de simuler la loi a posteriori, soit la loi jointe $q,\delta|x$ composante par composante via l'algorithme suivant
On initialise aléatoirement par exemple via la loi a priori $(q_1,\delta_1)\sim\pi(q,\delta)$
Pour $i = 1,\ldots, I$
Le résultat est une suite de réalisations $(q_1,\delta_1),\ldots, (q_I,\delta_I)$ d'une chaine de Markov dont la loi stationnaire est la loi a posteriori.
Cet méthode d'échantillonage, dite de Gibbs, requiert la connaissance des lois conditionnelles de $q|\delta,x$ et $\delta|q,x$ qui sont potentiellement plus facile à déterminer que la loi a posteriori $\pi(q,\delta|x)$.
Supposons que $\delta\sim \text{Inverse-Gamma}(\alpha,\beta)$ de densité
$$ \pi(\delta) = \frac{\beta^\alpha e^{-\beta/\delta}}{\Gamma(\alpha)\delta^{\alpha+1}},\text{ pour }\delta>0. $$Adaptons le code précédent pour échantilloner la loi a posteriori via un échantilloneur de Gibbs. On suppose que $\delta\sim \text{Inverse-Gamma}(a_\delta,b_\delta)$ et $q\sim \text{Unif}([a_q, b_q])$. La mise à jour du paramètre $q$ s'effectue via des itération Metropolis Hasting. Ecrire la fonction sim_post_sample_Gibbs
parms_name = ['q', 'δ']
a_δ, b_δ, a_q, b_q = 1, 1, 0, 1
σ, n_step = 0.01, 10000
def sim_post_sample_Gibbs(X, parms_name, a_δ, b_δ, a_q, b_q, σ, n_step):
logp = logp_geom_exp(X)
t, t0 = len(X), np.sum(X == 0)
trace = []
accepted = []
q0, δ0 = np.random.uniform(low = a_q, high = b_q, size = 1), 1 / np.random.gamma(shape = a_δ, scale = 1 / b_δ, size = 1)
trace.append(np.array([q0[0], δ0[0]])), accepted.append(True)
for j in range(n_step):
q_old, δ_old = trace[-1]
δ_new = 1 / np.random.gamma(shape = (t - t0) + a_δ + 1 , scale = 1 / ((1- q_old)*np.sum(X) + b_δ), size = 1)
q_new = np.random.normal(q_old, σ)
u = np.random.uniform(0, 1, size = 1)
if logp(np.array([q_new, δ_new[0]])) - logp(np.array([q_old, δ_old])) > u:
trace.append(np.array([q_new, δ_new[0]])), accepted.append(True)
else:
trace.append(np.array([q_old, δ_new[0]])), accepted.append(False)
trace_df = pd.DataFrame(trace)
trace_df.columns = parms_name
trace_df['accepted'] = np.array(accepted)
return(trace_df)
trace = sim_post_sample_Gibbs(X, parms_name, a_δ, b_δ, a_q, b_q, σ, n_step)
fig, axs = plt.subplots(1,2)
for j in range(len(parms_name)):
axs[j].plot(trace[parms_name[j]].values)
axs[j].set_title(parms_name[j])
sns.despine()
fig.tight_layout()
trace['q'].values
array([0.06620961, 0.09634418, 0.09634418, ..., 0.08984419, 0.08984419,
0.08984419])
sns.pairplot(trace[parms_name].loc[2500:n_step])
sns.despine()
Lorsque la vraisemblance n'est pas accessible, (ce serait le cas en considérant d'autres distributions pour la fréquence et le montant des sinistres) il existe une technique pour s'affranchir de l'absence de vraisemblance via des simulations. L'algorithme est le suivant
Tant que $i < I $
où $D(\cdot,\cdot)$ est une mesure de dissimilarité entre les données observées $\mathbf{x}$ et simulées $\mathbf{x}^\ast$ Le résultat est un échantillon distribué suivant la loi a posteriori approchée
$$ \pi_{abc}(\theta|\mathbf{x})=\frac{ \int_{\mathbb{R}^t}L(\mathbf{x}^\ast|\theta)\mathbb{I}_{D(\mathbf{x}, \mathbf{x^\ast}) <\epsilon}\text{d}\mathbf{x}^\ast\pi(\theta)}{\int_{\Theta}\int_{\mathbb{R}^t}L(\mathbf{x}^\ast|\theta)\mathbb{I}_{D(\mathbf{x}, \mathbf{x^\ast}) <\epsilon}\text{d}\mathbf{x}^\ast\pi(\theta)\text{d}\theta} $$La précision est bien entendu liée au choix de $\epsilon$, il s'agit d'un compromis précision/temps de calcul. Le paramètre important est la distance $D$ permettant la comparaison des données observées et simulées. Le choix d'une distance euclidienne $$ D(x,x^\ast) = \sqrt{\sum_{s = 1}^t(x_s - x_s^\ast)^2} $$
conduit a ne jamais sélectionner de paramètre eu égard à une trop grande variance (pour de grands échantillons). Une solution consiste à définir des résumés statistiques
$$ S:x\mapsto S(x)\in\mathbb{R}^{d},\text{ avec } d\leq t $$puis de mesurer la dissimilarité entre données observées et données simulées par $D[S(x),S(x^\ast)]$. L'utilisation de statistiques rajoute une approximation puisque la loi a posteriori approché converge vers la loi des paramètres sachant $S(x)$ ce qui correspond à une perte d'information.
Proposer une distance $D$ pour comparer les données simulées et observées. Ecrire une fonction distance prenant en argument le nom d'une distance et renvoyant une fonction prenant en argument $X$ et $X^\ast$ calculant la distance $D[S(x),S(x^\ast)]$.
def distance(name):
def D(X, X_ast):
if name == "Wasserestein":
return(np.sqrt(np.sum(np.sort(X)-np.sort(X_ast))**2))
elif name == "stat_ex":
t0, t0_ast = np.sum(X == 0), np.sum(X_ast == 0)
return(np.sqrt( (t0 - t0_ast)**2 + (np.sum(X) - np.sum(X_ast))**2))
return(D)
Ecrire une fonction sim_post_sample_ABC produisant des échantillons suivant la loi a posteriori approchée.
def sim_post_sample_ABC(X, parms_name, a_q, b_q, a_δ, b_δ, n_step, ϵ, name_distance):
trace, i = [], 0
D = distance(name_distance)
while i < n_step:
q_ast, δ_ast = np.random.uniform(low = a_q, high = b_q, size = 1), np.random.uniform(low = a_δ, high = b_δ, size = 1)
X_ast = sim_X(t, q_ast, δ_ast)
if D(X, X_ast) < ϵ :
trace.append(np.array([q_ast[0], δ_ast[0]]))
i = i +1
trace_df = pd.DataFrame(trace)
trace_df.columns = parms_name
return(trace_df)
t, q, δ = 100, 0.1, 5
X = sim_X(t, q, δ)
parms_name = ['q', 'δ']
a_q, b_q, a_δ, b_δ = 0, 1, 0, 10
n_step_abc, ϵ_abc, name_distance = 1000, 10, "stat_ex"
Σ, n_step_MH =np.diag(np.array([0.01, 1])), 5000
trace_abc, trace_MH = sim_post_sample_ABC(X, parms_name, a_q, b_q, a_δ, b_δ, n_step_abc, ϵ_abc, name_distance), sim_post_sample_MH(X, parms_name, a_q, b_q, a_δ, b_δ, Σ, n_step_MH)
trace_abc["method"] = np.repeat("abc",len(trace_abc) )
trace_MH["method"] = np.repeat("MH",len(trace_MH))
trace = pd.concat([trace_MH.drop(columns = ['accepted']).iloc[int(n_step_MH - n_step_abc):n_step_MH], trace_abc])
sns.pairplot(trace[np.append(parms_name, "method")], hue = "method")
sns.despine()
