MCMC et PyMC

Avance 60 min 12 sections

Quand le posterior n'a pas de formule, on l'echantillonne. Codez Metropolis-Hastings a la main en numpy, puis passez a l'outil professionnel PyMC et aux diagnostics ArviZ, jusqu'au test A/B complet.

Telecharger CSV Telecharger 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 chaine: burn-in, taux d'acceptation, autocorrelation
Utiliser PyMC (NUTS) et lire les diagnostics ArviZ (r_hat, ESS)
Mener un test A/B bayesien de bout en bout avec PyMC

Prerequis

Modules Inference Bayesienne et Regression Lineaire Bayesienne

Theorie

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

Les deux modules precedents vivaient dans un monde confortable: prior conjugue,

posterior avec formule exacte (Beta, gaussienne). Ce confort est l'exception, pas la regle.

Ou le confort s'arrete. Des que le modele devient realiste, le posterior

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

bute sur son denominateur: une integrale sur TOUTES les valeurs possibles des parametres.

Avec 2 parametres, on peut encore quadriller numeriquement. Avec 50 (un modele hierarchique

modeste), l'integrale est hors de portee de toute machine: c'est la malediction de la dimension.

L'idee geniale de MCMC (Markov Chain Monte Carlo):

On n'a pas besoin de CALCULER le posterior. On a besoin d'en tirer des ECHANTILLONS:

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 maitre: pour echantillonner, le numerateur $P(\text{donnees} \mid \theta) P(\theta)$

suffit. L'integrale impossible du denominateur disparait dans l'algorithme, car il ne

manipule que des RAPPORTS de densites entre deux points (le denominateur, identique en haut

et en bas, se simplifie).

La famille d'algorithmes:

Algorithme	Idee	Usage
Metropolis-Hastings (1953-1970)	marche aleatoire + regle d'acceptation	pedagogie, cas simples
Gibbs (1984)	echantillonner un parametre a la fois	modeles 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 a la main pour comprendre la mecanique,

puis on delegue a PyMC/NUTS pour le travail serieux.

Theorie

Schema: un pas de Metropolis-Hastings

L'algorithme 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 probabilite r"| ACC RATIO -->|"sinon"| REJ ACC -.-> CUR REJ -.-> CUR style CUR fill:#E5D7F5,color:#1A1A1A style PROP fill:#d4edda,color:#1A1A1A style RATIO fill:#9B7AC4,color:#fff style ACC fill:#F7E64D,color:#1A1A1A style REJ fill:#f8d7da,color:#1A1A1A

La chaine monte toujours vers les zones plus probables, mais accepte PARFOIS de descendre

(avec probabilite $r$): c'est ce qui lui permet d'explorer toute la distribution au lieu

de rester coincee sur le sommet.

Exemple numerique d'un pas:

La chaine 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}$

ou 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 chaine refusera cette excursion vers une zone peu probable,

un quart du temps elle ira quand meme: exactement la dose d'exploration necessaire pour

que l'histogramme des positions visitees reproduise le posterior.

Legende des couleurs:

$\textcolor{#9B7AC4}{Violet}$ : position actuelle de la chaine
$\textcolor{#3498db}{Bleu}$ : position proposee
$\textcolor{#27ae60}{Vert}$ : ratio favorable (acceptation garantie)
$\textcolor{#e74c3c}{Rouge}$ : ratio defavorable (acceptation probabiliste)

Avance Exercice manuel: A vous de calculer!

## Trois pas de Metropolis-Hastings a la main

La chaine explore un posterior dont voici quelques valeurs (non normalisees):

__MATH_e4ba04ae__	0.08	0.10	0.12	0.14	0.16
__MATH_61a7cd8e__	0.0020	0.0040	0.0060	0.0045	0.0015

La chaine demarre en $\theta = 0.10$. Pour chaque proposition, calculez le ratio

$r = \tilde{P}(\theta^*) / \tilde{P}(\theta_t)$ et dites si elle est acceptee a coup sur,

ou avec quelle probabilite.

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

Pas 2: la chaine est maintenant en $0.12$; on propose $\theta^* = 0.16$.

Le tirage uniforme donne $u = 0.31$. La proposition est-elle acceptee ?

Pas 3: d'ou repart le pas suivant ?

Question bonus: pourquoi n'a-t-on jamais eu besoin de la constante de

normalisation du posterior dans ces calculs ?

Avance Solution de l'exercice manuel

## Solution detaillee

Pas 1: de 0.10 vers 0.12

$$r = \frac{0.0060}{0.0040} = 1.5 \ge 1 \quad \Rightarrow \quad \text{ACCEPTE a coup sur}$$

La chaine 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 REJETEE.

Pas 3:

En cas de rejet, la chaine RESTE sur place: le pas suivant repart de $\theta = 0.12$

(et la valeur $0.12$ est comptee une fois de plus dans l'histogramme: les rejets

"epaississent" les zones probables, c'est voulu).

Question bonus:

Le posterior normalise s'ecrit $P(\theta) = \tilde{P}(\theta) / Z$ ou

$Z = \int \tilde{P}(\theta') d\theta'$ est l'integrale 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'etre de MCMC: explorer une distribution

qu'on ne sait pas normaliser.

Code

Metropolis-Hastings code a la main

Cliquez sur "Executer" pour voir le resultat

Contenu verrouille section restante

5 / 12

Continuez votre apprentissage

Vous avez explore 5 sections de ce module. Connectez-vous pour debloquer le reste du cours, incluant les exercices pratiques et les solutions.

Se connecter ou Creer un compte gratuit