Cours n°01 : Introduction à l’inférence bayésienne Cours n°02 : Modèle Beta-Binomial Cours n°03 : Introduction à brms, modèle de régression linéaire Cours n°04 : Modèle de régression linéaire (suite) Cours n°05 : Markov Chain Monte Carlo Cours n°06 : Modèle linéaire généralisé Cours n°07 : Comparaison de modèles Cours n°08 : Modèles multi-niveaux Cours n°09 : Modèles multi-niveaux généralisés Cours n°10 : Data Hackathon
\[\newcommand\given[1][]{\:#1\vert\:}\]
Rappels de notation
La notation \(p(y \given \theta)\) peut faire référence à deux choses selon le contexte : la fonction de vraisemblance et le modèle d’observation. De plus, on trouve de nombreuses notations ambigues en statistique. Essayons de clarifier ci-dessous.
\(\Pr(Y = y \given \Theta = \theta)\) désigne une probabilité (e.g., dbinom(x = 2, size = 10, prob = 0.5)).
\(p(Y = y \given \Theta = \theta)\) désigne une densité de probabilité (e.g., dbeta(x = 0.4, shape1 = 2, shape2 = 3)).
\(p(Y = y \given \Theta)\) désigne une fonction de vraisemblance (likelihood) discrète ou continue, \(y\) est connu/fixé, \(\Theta\) est une variable aléatoire, la somme (ou l’intégrale) de cette distribution n’est pas égale à 1 (e.g., dbinom(x = 2, size = 10, prob = seq(0, 1, 0.1) )).
\(p(Y \given \Theta = \theta)\) désigne une fonction de masse (ou densité) de probabilité (dont la somme ou l’intégrale est égale à 1), qu’on appelle aussi “modèle d’observation” (observation model) ou “distribution d’échantillonnage” (sampling distribution), \(Y\) est une variable aléatoire, \(\theta\) est connu/fixé (e.g., dbinom(x = 0:10, size = 10, prob = 0.5))
Le but de l’analyse bayésienne (i.e., ce qu’on obtient à la fin d’une telle analyse) est la distribution postérieure \(p(\theta \given y)\). On peut la résumer pour faciliter la communication des résultats, mais toute l’information souhaitée est contenue dans la distribution toute entière (pas seulement sa moyenne, son mode, ou autre).
######################################################################### On définit un modèle avec : ## Une fonction de vraisemblance Gaussienne : y ~ Normal(mu, sigma) ## Un prior Gaussien pour la moyenne : mu ~ Normal(100, 10) ## Et un prior Exponentiel pour l'écart-type : sigma ~ Exponential(0.1) ########################################################################## on simule 10.000 observations issues d'une distribution Gaussienne sans incertitude (épistémique)rnorm(n =1e4, mean =100, sd =10) |>hist(breaks ="FD")# on tire 10.000 échantillons issus du prior Gaussien pour mu (i.e., p(mu))# ce prior représente ce qu'on sait de mu avant de voir les données...mu_prior <-rnorm(n =1e4, mean =100, sd =10)# 10.000 observations issues d'une distribution Gaussienne avec incertitude sur murnorm(n =1e4, mean = mu_prior, sd =10) |>hist(breaks ="FD")# on tire 10.000 échantillons issus du prior Exponentiel pour sigma (i.e., p(sigma))# ce prior représente ce qu'on sait de sigma avant de voir les données...sigma_prior <-rexp(n =1e4, rate =0.1)# 10.000 observations issues d'une distribution Gaussienne avec incertitude sur mu ET sigma# ce que le modèle suppose à propos de y sur la base de nos priors pour mu et sigma...rnorm(n =1e4, mean = mu_prior, sd = sigma_prior) |>hist(breaks ="FD")
Petit problème : La constante de normalisation (en vert) s’obtient en calculant la somme (pour des variables discrètes) ou l’intégrale (pour des variables continues) de la densité conjointe \(p(\text{data}, \theta)\) sur toutes les valeurs possibles de \(\theta\). Cela se complique lorsque le modèle comprend plusieurs paramètres et/ou que la forme de la distribution postérieure est complexe…
Le problème avec la distribution postérieure
Rappels Cours n°02
Trois méthodes pour résoudre (contourner) ce problème :
La distribution a priori est un prior conjugué de la fonction de vraisemblance (e.g., modèle Beta-Binomial). Dans ce cas, il existe une solution analytique (i.e., qu’on peut calculer de manière exacte) pour la distribution postérieure.
Autrement, pour des modèles simples, on peut utiliser la méthode par grille. On calcule la valeur exacte de la probabilité postérieure en un nombre fini de points dans l’espace des paramètres.
Pour les modèles plus complexes, explorer tout l’espace des paramètres n’est pas tractable. On va plutôt échantillonner intelligemment un grand nombre de points dans l’espace des paramètres.
Objectifs du cours
\(\longrightarrow~\) Présenter le principe de base de l’échantillonnage : Markov Chain Monte Carlo
\(\longrightarrow~\) Présenter deux algorithmes (Metropolis-Hastings et HMC)
\(\longrightarrow~\) Montrer les forces mais aussi les faiblesses de ces méthodes
\(\longrightarrow~\) Donner des outils de contrôle sur ces méthodes
\(\longrightarrow~\) Appliquer ces méthodes à un cas simple
Markov Chain Monte Carlo
Markov chain Monte Carlo \(\longrightarrow~\) Échantillonnage aléatoire \(\longrightarrow~\) Le résultat est un ensemble de valeurs du paramètre
Markov chain Monte Carlo \(\longrightarrow~\) Les valeurs sont générées sous forme de séquences (liaison de dépendance) \(\longrightarrow~\) Indice temporel pour identifier la place dans la chaîne \(\longrightarrow~\) Le résultat est de la forme : \(\theta^1, \theta^2, \theta^3, \dots, \theta^t\)
Markov chain Monte Carlo \(\longrightarrow~\) La valeur de paramètre générée ne dépend que de la valeur du paramètre précédent \(\Pr(\theta^{t+1} \given \theta^{t}, \theta^{t-1}, \ldots, \theta^{1}) = \Pr(\theta^{t + 1} \given \theta^{t})\)
Méthodes Monte Carlo
Le terme de méthode de Monte-Carlo désigne une famille d’algorithmes visant à calculer (ou approcher) une valeur numérique en utilisant des procédés aléatoires, c’est-à-dire des techniques probabilistes. Cette méthode a été formalisée en 1947 par Nicholas Metropolis, et publiée pour la première fois en 1949 dans un article co-écrit avec Stanislaw Ulam.
Méthodes Monte Carlo : Estimation de \(\pi\)
Soit un point \(M\) de coordonnées \((x, y)\), où \(0 < x < 1\) et \(0 < y < 1\). On tire aléatoirement les valeurs de \(x\) et \(y\) entre \(0\) et \(1\) suivant une loi uniforme. Le point \(M\) appartient au disque de centre \((0, 0)\) de rayon \(r = 1\) si et seulement si \(\sqrt{x^{2} + y^{2}} \leqslant 1\). On sait que le quart de disque est de surface \(\sigma = \pi r^{2} / 4 = \pi / 4\) et que le carré qui le contient est de surface \(s = r^{2} = 1\). Si la loi de probabilité du tirage de point est uniforme, la probabilité que le point \(M\) appartienne au disque est donc de \(\sigma / s = \pi / 4\). En faisant le rapport du nombre de points dans le disque au nombre de tirages \(\frac{N_{\text{inner}}}{N_{\text{total}}}\), on obtient alors une approximation de \(\pi / 4\).
Méthodes Monte Carlo : Estimation de \(\pi\)
trials <-1e5# nombre d'échantillonsradius <-1# rayon du cerclex <-runif(n = trials, min =0, max = radius) # tirages pour xy <-runif(n = trials, min =0, max = radius) # tirages pour ydistance <-sqrt(x^2+ y^2) # distance à l'origineinside <- distance < radius # à l'intérieur (ou pas) du quart de cercle ?pi_estimate <-4*sum(inside) / trials # estimation de pi
