MCMC et PyMC

Avancé 60 min 13 sections

Quand le posterior n'a plus de formule, on l'échantillonne : Metropolis-Hastings à la main, puis NUTS et ArviZ.

Télécharger le CSV Télécharger le script Python

Objectifs d'apprentissage

Comprendre pourquoi et quand on a besoin de MCMC
Coder l'algorithme de Metropolis-Hastings en numpy pur
Diagnostiquer une chaîne : burn-in, taux d'acceptation, autocorrélation
Utiliser PyMC (NUTS) et lire les diagnostics ArviZ (divergences, r_hat, ESS)
Reconnaître une divergence, comprendre ce qu'elle signale, et la corriger par reparamétrisation plutôt que par target_accept
Mener un test A/B bayésien de bout en bout avec PyMC

Prérequis

Modules Inférence bayésienne et Régression linéaire bayésienne

Théorie

Pourquoi MCMC: quand le posterior n'a plus de formule

Les deux modules précédents vivaient dans un monde confortable: prior conjugue, posterior avec formule exacte (Beta, gaussienne). Ce confort est l'exception, pas la règle.

Où le confort s'arrête. Des que le modèle devient réaliste, le posterior

$$P(\theta \mid \text{données}) = \frac{P(\text{données} \mid \theta) \, P(\theta)}{\int P(\text{données} \mid \theta') \, P(\theta') \, d\theta'}$$

bute sur son dénominateur: une intégrale sur TOUTES les valeurs possibles des paramètres. Avec 2 paramètres, on peut encore quadriller numériquement. Avec 50 (un modèle hiérarchique modeste), l'intégrale est hors de portée de toute machine: c'est la malédiction de la dimension.

L'idée géniale de MCMC (Markov Chain Monte Carlo):

On n'a pas besoin de CALCULER le posterior. On a besoin d'en tirer des ÉCHANTILLONS: avec 10000 tirages de $\theta$, on estime tout ce qu'on veut (moyenne, intervalles, $P(\theta > 0.1)$...) par simple comptage.

Et le coup de maître: pour échantillonner, le numérateur $P(\text{données} \mid \theta) P(\theta)$ suffit. L'intégrale impossible du dénominateur disparaît dans l'algorithme, car il ne manipule que des RAPPORTS de densités entre deux points (le dénominateur, identique en haut et en bas, se simplifie).

La famille d'algorithmes:

Algorithme	Idée	Usage
Metropolis-Hastings (1953-1970)	marche aléatoire + règle d'acceptation	pédagogie, cas simples
Gibbs (1984)	échantillonner un paramètre à la fois	modèles conditionnellement conjugues
HMC / NUTS (2011)	utiliser le GRADIENT pour proposer loin et bien	le standard actuel (PyMC, Stan)

Dans ce module: on code Metropolis-Hastings à la main pour comprendre la mécanique, puis on délègue à PyMC/NUTS pour le travail sérieux.

Théorie

Schéma: un pas de Metropolis-Hastings

À chaque itération, la chaîne propose un petit déplacement au hasard depuis sa position actuelle, puis décide: soit elle accepte la proposition (elle se déplace), soit elle la rejette (elle reste sur place). Voici ce pas en une image:

flowchart TD CUR["Position actuelle theta_t"] PROP["Proposition
theta* = theta_t + bruit"] RATIO["Ratio
r = post(theta*) / post(theta_t)"] ACC["ACCEPTER
theta_t+1 = theta*"] REJ["REJETER
theta_t+1 = theta_t"] CUR --> PROP --> RATIO RATIO -->|"r >= 1
(on monte)"| ACC RATIO -->|"r < 1: accepter
avec probabilité r"| ACC RATIO -->|"sinon"| REJ ACC -.-> CUR REJ -.-> CUR class CUR ml-node-secondary class PROP ml-node-success class RATIO ml-node-main class ACC ml-node-accent class REJ ml-node-danger

Accepter, rejeter: qu'est-ce que ça veut dire concrètement ?

Accepter = la chaîne SE DÉPLACE: la proposition $\theta^*$ devient la nouvelle position ($\theta_{t+1} = \theta^*$) et elle est ajoutée à la liste des échantillons.
Rejeter = la chaîne RESTE SUR PLACE: la position actuelle est conservée ($\theta_{t+1} = \theta_t$) et elle est enregistrée une nouvelle fois dans les échantillons.

Dans les deux cas, on note une valeur à chaque pas. Rejeter n'est donc pas "perdre un tour": c'est voter une fois de plus pour la position actuelle. C'est ainsi que les zones probables accumulent plus d'échantillons que les zones improbables.

La chaîne monte toujours vers les zones plus probables, mais accepte PARFOIS de descendre (avec probabilité $r$): c'est ce qui lui permet d'explorer toute la distribution au lieu de rester coincée sur le sommet.

Exemple numérique d'un pas:

La chaîne est en $\textcolor{#9B7AC4}{\theta_t = 0.10}$ ou le posterior (non normalise) vaut $\textcolor{#9B7AC4}{0.0040}$. On propose $\textcolor{#3498db}{\theta^* = 0.13}$ où il vaut $\textcolor{#3498db}{0.0052}$.

$$r = \frac{\textcolor{#3498db}{0.0052}}{\textcolor{#9B7AC4}{0.0040}} = \textcolor{#27ae60}{1.3} \ge 1 \quad \Rightarrow \quad \text{on ACCEPTE, } \theta_{t+1} = 0.13$$

Pas suivant: de $\theta_t = 0.13$ on propose $\theta^* = 0.20$ ou le posterior vaut $0.0013$.

$$r = \frac{0.0013}{0.0052} = \textcolor{#e74c3c}{0.25} < 1 \quad \Rightarrow \quad \text{on tire } u \sim \text{Uniforme}(0,1) \text{ et on accepte si } u < 0.25$$

Trois quarts du temps la chaîne refusera cette excursion vers une zone peu probable, un quart du temps elle ira quand même: exactement la dose d'exploration nécessaire pour que l'histogramme des positions visitées reproduise le posterior.

Légende des couleurs:

$\textcolor{#9B7AC4}{Violet}$ : position actuelle de la chaîne
$\textcolor{#3498db}{Bleu}$ : position proposée
$\textcolor{#27ae60}{Vert}$ : ratio favorable (acceptation garantie)
$\textcolor{#e74c3c}{Rouge}$ : ratio défavorable (acceptation probabiliste)

Avancé Exercice manuel: À vous de calculer!

Trois pas de Metropolis-Hastings à la main

La chaîne explore un posterior dont voici quelques valeurs (non normalisées):

$\theta$	0.08	0.10	0.12	0.14	0.16
$\tilde{P}(\theta)$	0.0020	0.0040	0.0060	0.0045	0.0015

La chaîne démarre en $\theta = 0.10$. Pour chaque proposition, calculez le ratio $r = \tilde{P}(\theta^*) / \tilde{P}(\theta_t)$ et dites si elle est acceptée à coup sûr, ou avec quelle probabilité.

Pas 1: depuis $\theta = 0.10$, on propose $\theta^* = 0.12$.

Pas 2: la chaîne est maintenant en $0.12$; on propose $\theta^* = 0.16$. Le tirage uniforme donne $u = 0.31$. La proposition est-elle acceptée ?

Pas 3: d'où repart le pas suivant ?

Question bonus: pourquoi n'a-t-on jamais eu besoin de la constante de normalisation du posterior dans ces calculs ?

Avancé Solution de l'exercice manuel

Solution détaillée

Pas 1: de 0.10 vers 0.12

$$r = \frac{0.0060}{0.0040} = 1.5 \ge 1 \quad \Rightarrow \quad \text{ACCEPTE à coup sûr}$$

La chaîne monte vers une zone plus probable: acceptation automatique. Position: $0.12$.

Pas 2: de 0.12 vers 0.16

$$r = \frac{0.0015}{0.0060} = 0.25 < 1 \quad \Rightarrow \quad \text{accepte seulement si } u < 0.25$$

Or $u = 0.31 > 0.25$: la proposition est REJETÉE.

Pas 3:

En cas de rejet, la chaîne RESTE sur place: le pas suivant repart de $\theta = 0.12$ (et la valeur $0.12$ est comptée une fois de plus dans l'histogramme: les rejets "épaississent" les zones probables, c'est voulu).

Question bonus:

Le posterior normalise s'écrit $P(\theta) = \tilde{P}(\theta) / Z$ ou $Z = \int \tilde{P}(\theta') d\theta'$ est l'intégrale incalculable. Dans le ratio:

$$r = \frac{\tilde{P}(\theta^*) / Z}{\tilde{P}(\theta_t) / Z} = \frac{\tilde{P}(\theta^*)}{\tilde{P}(\theta_t)}$$

$Z$ se simplifie. C'est LA raison d'être de MCMC: explorer une distribution qu'on ne sait pas normaliser.

Code

Metropolis-Hastings code à la main

Ctrl+Entrée

Cliquez sur "Exécuter" pour voir le résultat

Contenu verrouillé

5 / 13

Continuez votre apprentissage

Vous avez exploré 5 sections de ce module. Connectez-vous pour débloquer le reste du cours, incluant les exercices pratiques et les solutions.

Se connecter ou Créer un compte gratuit