Une difficile estimation de la distribution a posteriori

Paramètres multidimensionnels

Dans les applications réelles il y a en général plusieurs paramètres. On a donc un vecteur de paramètres \(\boldsymbol{\theta} = (\theta_1,\dots, \theta_d)\). La formule de Bayes donnant la distribution a posteriori à partir de la loi a priori et de la vraisemblance est toujours valable dans ce cas : elle donne la loi a posteriori conjointe des \(d\) paramètres. Toute l’information est contenue dans cette loi conjointe. Mais son calcul numérique n’est pas toujours facile, notamment dans des modèles complexes, et dans certains modèles même la vraisemblance est difficile à calculer. De plus, pour obtenir la distribution a posteriori conjointe, il faut aussi calculer la constante d’intégration \(f(\boldsymbol{y}) = \textstyle \int_{\Theta^d} f(\boldsymbol{y}|\boldsymbol{\theta})\pi(\boldsymbol{\theta})\,\text{d}\boldsymbol{\theta}\). Une solution analytique n’est disponible que dans des cas très particuliers (notamment lors de l’utilisation de lois conjuguées); et dans la majorité des cas pratiques, l’intégrale doit en fait être calculée numériquement. Si \(\theta\) est de dimension \(d\), il s’agit de calculer une intégrale de multiplicité \(d\), ce qui devient difficile lorsque \(d\) est grand (les problèmes numériques apparaissent dès que \(d> 4\)).

Un problème encore plus difficile surgit lorsque l’on veut tirer des conclusions à partir de cette distribution a posteriori conjointe. En général nous sommes intéressés par les valeurs possibles pour chaque paramètre. C’est à dire que l’on a besoin des distributions marginales, unidimensionnelles, de chaque paramètre. Pour les obtenir il faut là aussi intégrer la distribution conjointe (et ce \(d\) fois s’il y a \(d\) paramètres). Le problème est d’autant plus difficile qu’il faut calculer ces intégrales pour chaque valeur possible du paramètre pour reconstituer toute la densité a posteriori numériquement. Dans les problèmes complexes un calcul suffisamment précis de ces intégrales parait impossible et l’on a recours en général à des algorithmes basés sur des simulations d’échantillonnage, en particulier les algorithmes de Monte-Carlo par chaînes de Markov (dits “Markov chain Monte Carlo” — MCMC).

Statistique bayésienne computationnelle

La statistique bayésienne computationnelle cherche des solutions pour pouvoir estimer la distribution a posteriori, y compris lorsqu’on ne connaît que le numérateur dans le théorème de Bayes (loi a posteriori non-normée). Les principales méthodes utilisées s’appuient sur des algorithmes d’échantillonnage permettant de générer un échantillon distribué selon la loi a posteriori. On peut distinguer deux grands types parmi ces algorithmes : i) d’abord les méthodes d’échantillonnage directes, où l’on génère un échantillon à partir d’une loi simple (par exemple uniforme), que l’on transforme afin que le résultat soit distribué selon la loi a posteriori ; ii) les méthodes de Monte-Carlo par chaînes de Markov (MCMC), où l’on construit une chaîne de Markov sur l’espace des paramètres dont la loi invariante correspond à la loi a posteriori.

Méthode de Monte-Carlo

Monte-Carlo (1955) est le nom crypté d’un projet de John von Neumann et Stanislas Ulam au Los Alamos Scientific Laboratory visant à utiliser des nombres aléatoires pour estimer des quantités difficiles (ou impossible) à calculer analytiquement.

En s’appuyant sur la Loi des Grands Nombres, il s’agit d’obtenir un échantillon dit “de Monte-Carlo” permettant de calculer divers fonctionnelles à partir de la distribution de probabilité suivie par l’échantillon. En effet \(\mathbb{E}[f(X)] = \int_x f(x)p_X(x)dx\). Or grâce à la loi (forte) des grands nombres, on a : \(\frac{1}{N}\sum_i f(x_i) \overset{p.s.}{\underset{N \to +\infty}{\longrightarrow}} \mathbb{E}[f(X)]\), à condition que les \(x_i\) forment un échantillon \(iid\) selon la loi de \(X\). On peut alors estimer un certain nombre d’intégrales, à condition d’être capable d’échantillonner selon \(p_X\).

En réalité, il est même possible d’affaiblir la condition d’indépendance dans la Loi des Grands Nombres grâce au théorème ergodique : il suffit que la suite de \(x_i\) soit ergodique (intuitivement, cela peut signifier que la dépendance décroît au fur et à mesure que les observations s’espacent – et à condition que l’espérance de la loi \(p_X\) soit finie).