 Insea 2025             Statistiques Bayésiennes

TP3: Introduction à PyMC

                 Author: Hicham Janati

Objectifs:

Découvrir la librairie PyMC
Implémenter les premiers modèles bayésiens et faire le diagnostic de convergence
Interpréter les résultats et comparer avec les statistiques fréquentistes

%pip install pymc arviz numpy pandas matplotlib ipywidgets

On importe les librariries et un crée un générateur aléatoire:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pymc as pm
import arviz as az


seed = 42
rng = np.random.default_rng(seed)

Modèle Poisson-Gamma simple

On considère l’ex 1 du TD1. On observe le nombre de sinistres par année \(N | \lambda \sim \mathcal{P}(\lambda)\), où \(\lambda\) suit une loi a priori \(\text{Gamma}(a, b)\). On suppose que les données historiques mènent au choix a = 4 et b = 2. On utilise la définition de Gamma où b correspond au “rate” et non pas au scale (comme en TD), la moyenne de cette Gamma est a/b. (voir wikipedia).

On suppose que les 5 observations sont données par:

data = np.array([4, 0, 2, 1, 0])
rng = np.random.default_rng(42)
lambda_true = 2
data = rng.poisson(lambda_true, size=5)

On définit le modèle bayésien dans un contexte avec pymc ou on précise la distribution, le nom et les paramètres de chaque variable

a = 4
b = 2
with pm.Model() as model:
    lambda_ = pm.Gamma("lambda", a, b) # non observée
    N_obs = pm.Poisson("N_obs", mu=lambda_, observed=data) # observée
    trace = pm.sample(1000) # on simule une chaine MCMC de la loi a posteriori

On voit que pymc a automatiquement choisi NUTS et a simulé 4 chaînes avec 2000 échantillons chacune dont 1000 jetés (tuning / burn-in). Voyons ce que contient l’objet trace:

trace

C’est un objet InferenceData du package arviz. On peut obtenir les échantillons simulés dans l’attribut posterior:

trace.posterior["lambda"].data

trace.posterior["lambda"].data.shape

On a effectivement généré 4 chaines avec 1000 échantillons chacune. On peut les visualiser:

plt.figure()
plt.plot(trace.posterior["lambda"].data.T)
plt.grid(True)
plt.title("Trace plot")
plt.show()

Ou utiliser la librairie arviz directement qui donne également une estimation de la densité a posteriori:

az.plot_trace(trace)

On fait le diagnostic de convergence:

az.plot_autocorr(trace)
plt.show()

az.summary(trace)

az.summary(trace)

Les tracés sont bien mélangés
Rhat = 1.0 < 1.01 : pas de différence significative entre les 4 chaines
ESS très larges
Autocorrélations diminuent très rapidement

Cette chaîne réussit les diagnostics de convergence.

On peut visualiser la densité a posteriori avec l’intervalle de crédibilité HDI:

az.plot_posterior(trace, hdi_prob=0.94)

Les bornes de cet intervalle sont également présente dans le tableau du summary ci-dessus. On peut calcule un HDI directement:

az.hdi(trace, hdi_prob=0.94).to_pandas()

On peut calculer un ESS relatif (divisé par le nombre d’échantillons):

az.ess(trace, method="bulk", relative=True).to_pandas()

Ainsi, on en déduit que 44% des échantillons “sont efficaces” pour estimer “le centre” (bulk) de la distribution

Question 1: Augmenter le nombre de sinistres observés de 5 à 10 puis 100, comment changent les statistiques du az.summary ?

Question 2: Les métriques de convergence sont-elle très différentes ? Est-ce surprenant ?

Question 3: Comment peut-on interpréter le HDI obtenu ?

Question 4: En utilisant le fait que la loi a priori soit conjuguée, générez des échantillons a posteriori directement (sans pymc) et comparez

La loi a posteriori est \(\text{Gamma}(a + \sum_{i=1}^n N_i, b + n)\)

lambda_post_samples = rng.gamma(a + data.sum(), 1/(b + len(data)), size=(4, 1000))

ax = az.plot_posterior(trace)
ax.hist(lambda_post_samples, bins=100, density=True, alpha=0.7)
plt.show()

trace_iid = az.convert_to_inference_data(dict(iid=lambda_post_samples))
az.summary(trace_iid)

lambda_post_samples.mean()

az.summary(trace)

Modèle Poisson-Gamma à plusieurs conducteurs

On considère désormais les données de plusieurs individus avec un \(\lambda_i\) différent mais un a priori commun. Chaque conducteur \(i\) a ses données \(N_i^1, \dots, N_i^m\). Le modèle pymc s’adapte facilement en changeant le shape des paramètres:

data = np.array([[4, 0, 0],
                [0, 2, 1],
                [2, 2, 0],
                [1, 3, 0],
                [0, 0, 1]])

# donnéees de 3 conducteurs

n_drivers = data.shape[1]

a = 4
b = 2
with pm.Model() as model:
    lambda_ = pm.Gamma("lambda", a, b, shape=n_drivers) # non observée, on précise le nombre de lambda
    # le shape des données par défaut est n_observation x n_features, pymc associe chaque lambda_i a une colonne de data
    N_obs = pm.Poisson("N_obs", mu=lambda_, observed=data, ) # observée
    trace = pm.sample(1000) # on simule une chaine MCMC de la loi a posteriori
    pm.compute_log_likelihood(trace)

L’objet trace contient désormais plusieurs variables lambda:

print(az.summary(trace))

Question 4: Comparez l’estimation fréquentiste avec l’estimation bayésienne avec (a, b) = (4, 2) puis (a, b) = (10, 1).

data.mean(0)

Question 5: Déterminez un intervalle de confiance fréquentiste de niveau 95% pour chaque \(\lambda_i\).

\([\bar{N} \pm \frac{Q_{97.5}\hat{\sigma}}{\sqrt{n}}]\) avec \(Q_{97.5}\) le quantile de la loi de student à n-1 degrés de liberté.

from scipy.stats import t

n = 5
q975 = t.ppf(0.975, df=n-1)
sigma = data.std(0)

data.mean(0) - q975 * sigma / n ** 0.5, data.mean(0) + q975 * sigma / n ** 0.5,

az.summary(trace)

Question 6: On reprend à présent une loi a priori Uniforme([0, 5]). Déterminez des HDI de niveau 95% pour chaque \(\lambda_i\). Comment se comparent-ils aux intervalles de confiance fréquentistes ?

data = np.array([[4, 0, 0],
                [0, 2, 1],
                [2, 2, 0],
                [1, 3, 0],
                [0, 0, 1]])

# donnéees de 3 conducteurs

n_drivers = data.shape[1]

a = 4
b = 2
with pm.Model() as model:
    lambda_ = pm.Uniform("lambda", 0, 5, shape=n_drivers) # non observée, on précise le nombre de lambda
    # le shape des données par défaut est n_observation x n_features, pymc associe chaque lambda_i a une colonne de data
    N_obs = pm.Poisson("N_obs", mu=lambda_, observed=data, ) # observée
    trace = pm.sample(1000) # on simule une chaine MCMC de la loi a posteriori

data = np.array([[4, 0, 0],
                [0, 2, 1],
                [2, 2, 0],
                [1, 3, 0],
                [0, 0, 1]])

# donnéees de 3 conducteurs

n_drivers = data.shape[1]

a = 4
b = 2
with pm.Model() as model_unif:
    # lambda_ = pm.Gamma("lambda", a, b, shape=n_drivers) # non observée, on précise le nombre de lambda
    lambda_ = pm.Uniform("lambda", 0, 100, shape=n_drivers)
    # le shape des données par défaut est n_observation x n_features, pymc associe chaque lambda_i a une colonne de data
    N_obs = pm.Poisson("N_obs", mu=lambda_, observed=data, ) # observée
    trace_unif = pm.sample(1000) # on simule une chaine MCMC de la loi a posteriori

az.summary(trace_unif)

data.mean(0) - q975 * sigma / n ** 0.5, data.mean(0) + q975 * sigma / n ** 0.5,