Un cours en R, Stan, et brms
Ladislas Nalborczyk (LPC, LNC, CNRS, Aix-Marseille Univ)
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\:}\]
Principes de l’analyse bayésienne :
The first idea is that Bayesian inference is reallocation of credibility across possibilities. The second foundational idea is that the possibilities, over which we allocate credibility, are parameter values in meaningful mathematical models (Kruschke, 2015).
Inférence bayésienne : On infère (ou plutôt on déduit) la probabilité que le paramètre ait telle ou telle valeur sachant les données (et le prior) via le théorème de Bayes.
\[ \color{purple}{p(\theta \given y)} = \frac{\color{orangered}{p(y \given \theta)}\ \color{steelblue}{p(\theta)}}{\color{green}{p(y)}} \propto \color{orangered}{p(y \given \theta)}\ \color{steelblue}{p(\theta)} \]
Objectif de l’analyse de données bayésienne : Faire évoluer une connaissance a priori sur les paramètres \(\color{steelblue}{p(\theta)}\) en une connaissance a posteriori \(\color{purple}{p(\theta \given y)}\), intégrant l’information contenue dans les nouvelles données via \(\color{orangered}{p(y \given \theta)}\).
Les étapes de l’analyse de données bayésienne :
1. Définir le modèle - Identifier les paramètres du modèle génératif, définir une distribution a priori pour ces paramètres.
2. Mettre à jour le modèle - Calculer la distribution a posteriori des paramètres (ou une bonne approximation).
3. Interpréter la distribution postérieure - Comparaison de modèles, estimation des paramètres, vérification des prédictions du modèle.
Objectif du cours : Illustrer les différentes étapes de cette démarche à l’aide d’un modèle simple (un seul paramètre), le modèle Beta-Binomial.
Pourquoi ce modèle ?
Le modèle Beta-Binomial couvre de nombreux problèmes de la vie courante :
C’est un modèle simple
S’applique à toutes les situations où le processus de génération des données ne peut résulter qu’en deux issues mutuellement exclusives (e.g., un lancer de pièce). À chaque essai, si on admet que \(\Pr(\text{face}) = \theta\), alors \(\Pr(\text{pile}) = 1 - \theta\).
Depuis Bernoulli, on sait calculer la probabilité du résultat d’un lancer de pièce, du moment que l’on connait le biais de la pièce \(\theta\). Admettons que \(Y = 0\) lorsqu’on obtient pile, et que \(Y = 1\) lorsqu’on obtient face. Alors \(Y\) est distribuée selon une loi de Bernoulli :
\[p(y \given \theta) = \Pr(Y = y \given \theta) = \theta^{y} (1 - \theta)^{(1 - y)}\]
En remplaçant \(y\) par \(0\) ou \(1\), on retombe bien sur nos observations précédentes :
\[\Pr(Y = 1 \given \theta) = \theta^{1} (1 - \theta)^{(1 - 1)} = \theta \times 1 = \theta\]
\[\Pr(Y = 0 \given \theta) = \theta^{0} (1 - \theta)^{(1 - 0)} = 1 \times (1 - \theta) = 1 - \theta\]
Si l’on dispose d’une suite de lancers \(\{Y_i\}\) indépendants et identiquement distribués (i.e., chaque lancer a une distribution de Bernoulli de probabilité \(\theta\)), l’ensemble de ces lancers peut être décrit par une distribution binomiale.
Par exemple, imaginons que l’on dispose de la séquence de cinq lancers suivants : Pile, Pile, Pile, Face, Face. On peut recoder cette séquence en \(\{0, 0, 0, 1, 1\}\).
Rappel : La probabilité de chaque \(1\) est \(\theta\) est la probabilité de chaque \(0\) est \(1 - \theta\).
Quelle est la probabilité d’obtenir 2 faces sur 5 lancers ?
Sachant que les essais sont indépendants les uns des autres, la probabilité d’obtenir cette séquence est de \((1 - \theta) \times (1 - \theta) \times (1 - \theta) \times \theta \times \theta\), c’est à dire : \(\theta^{2} (1 - \theta)^{3}\).
On peut généraliser ce résultat pour une séquence de \(n\) lancers et \(y\) “succès” :
\[\theta^{y} (1 - \theta)^{n - y}\]
On a jusqu’ici considéré seulement une seule séquence résultant en 2 succès pour 5 lancers, mais il existe de nombreuses séquences pouvant résulter en 2 succès pour 5 lancers (e.g., \(\{0, 0, 1, 0, 1\}\), \(\{0, 1, 1, 0, 0\}\))…
Le coefficient binomial nous permet de calculer le nombre de combinaisons possibles résultant en \(y\) succès pour \(n\) lancers de la manière suivante (lu “\(y\) parmi \(n\)” ou “nombre de combinaisons de \(y\) parmi \(n\)”)1 :
\[ \left(\begin{array}{l} n \\ y \end{array}\right) = C_{y}^{n} = \frac{n !}{y !(n - y) !} \]
Par exemple pour \(y = 1\) et \(n = 3\), on sait qu’il existe 3 combinaisons possibles : \(\{0, 0, 1\}, \{0, 1, 0\}, \{1, 0, 0\}\). On peut vérifier ça par le calcul, en appliquant la formule ci-dessus.
\[ \left(\begin{array}{l} 3 \\ 1\end{array}\right) = C_{1}^{3} = \frac{3 !}{1 !(3 - 1) !} = \frac{3 \times 2 \times 1}{1 \times 2 \times 1} = \frac{6}{2} = 3 \]
\[ p(y \given \theta) = \Pr(Y = y \given \theta) = \left(\begin{array}{l} n \\ y \end{array}\right) \theta^{y}(1 - \theta)^{n - y} \]
La loi binomiale nous permet de calculer la probabilité d’obtenir \(y\) succès sur \(n\) essais, pour un \(\theta\) donné. Exemple de la distribution binomiale pour une pièce non biaisée (\(\theta = 0.5\)), indiquant la probabilité d’obtenir \(n\) faces sur 10 lancers (en R: dbinom(x = 0:10, size = 10, prob = 0.5)
).
library(tidyverse)
set.seed(666) # for reproducibility
sample(x = c(0, 1), size = 500, prob = c(0.4, 0.6), replace = TRUE) %>% # theta = 0.6
data.frame() %>%
mutate(x = seq_along(.), y = cummean(.) ) %>%
ggplot(aes(x = x, y = y) ) +
geom_line(lwd = 1) +
geom_hline(yintercept = 0.6, lty = 3) +
labs(x = "Nombre de lancers", y = "Proportion de faces") +
ylim(0, 1)
Fonction de vraisemblance (likelihood)
La fonction de vraisemblance s’écrit de la manière suivante :
\[ \color{orangered}{\mathcal{L}(\theta \given y, n) = p(y \given \theta, n) = \left(\begin{array}{l} n \\ y \end{array}\right) \theta^{y}(1 - \theta)^{n - y} \propto \theta^{y}(1 - \theta)^{n - y}} \]
On lance à nouveau une pièce de biais \(\theta\) (où \(\theta\) représente la probabilité d’obtenir Face). On lance cette pièce deux fois et on obtient une Face et un Pile.
On peut calculer la probabilité d’observer une Face sur deux lancers de pièce en fonction de différentes valeurs de \(\theta\) de la manière suivante :
\[ \begin{aligned} \Pr(F, P \given \theta) + \Pr(P, F \given \theta) &= 2 \times \Pr(P \given \theta) \times \Pr(F \given \theta) \\ &= \theta(1 - \theta) + \theta(1 - \theta) \\ &= 2 \theta(1 - \theta) \end{aligned} \]
Cette probabilité est définie pour un jeu de données fixe et une valeur de \(\theta\) variable. On peut représenter cette fonction visuellement.
# Représentation graphique de la fonction de vraisemblance de theta pour y = 1 et n = 2
y <- 1 # nombre de faces
n <- 2 # nombre d'essais
data.frame(theta = seq(from = 0, to = 1, length.out = 1e3) ) %>%
mutate(likelihood = dbinom(x = y, size = n, prob = theta) ) %>%
ggplot(aes(x = theta, y = likelihood) ) +
geom_area(color = "orangered", fill = "orangered", alpha = 0.5) +
labs(x = expression(paste(theta, " - Pr(face)") ), y = "Vraisemblance (likelihood)")
Si on calcule l’aire sous la courbe de cette fonction, on obtient :
\[\int_{0}^{1} 2 \theta(1 - \theta) \mathrm{d} \theta = \frac{1}{3}\]
Quand on varie \(\theta\), la fonction de vraisemblance n’est pas une distribution de probabilité valide (i.e., son intégrale n’est pas égale à 1). On utilise le terme de vraisemblance, pour distinguer ce type de fonction des fonctions de densité de probabilité. On utilise la notation suivante pour mettre l’accent sur le fait que la fonction de vraisemblance est une fonction de \(\theta\), et que les données sont fixes : \(\mathcal{L}(\theta \given \text{data}) = p(\text{data} \given \theta)\).
Nombre de Faces (y) |
||||
---|---|---|---|---|
theta | 0 | 1 | 2 | Total |
0 | 1.00 | 0.00 | 0.00 | 1 |
0.2 | 0.64 | 0.32 | 0.04 | 1 |
0.4 | 0.36 | 0.48 | 0.16 | 1 |
0.6 | 0.16 | 0.48 | 0.36 | 1 |
0.8 | 0.04 | 0.32 | 0.64 | 1 |
1 | 0.00 | 0.00 | 1.00 | 1 |
Total | 2.20 | 1.60 | 2.20 |
Notons que la vraisemblance de \(\theta\) pour une donnée particulière est égale à la probabilité de cette donnée pour cette valeur de \(\theta\). Cependant, la distribution de ces vraisemblances (en colonne) n’est pas une distribution de probabilité. Dans l’analyse bayésienne, les données sont considérées comme fixes et la valeur de \(\theta\) est considérée comme une variable aléatoire.
Comment définir un prior dans le cas du lancer de pièce ?
Aspect sémantique \(~\rightarrow~\) le prior doit pouvoir rendre compte :
Aspect mathématique \(~\rightarrow~\) pour une solution entièrement analytique :
\[ \begin{align} \color{steelblue}{p(\theta \given a, b)} \ &\color{steelblue}{= \mathrm{Beta}(\theta \given a, b)} \\ & \color{steelblue}{= \theta^{a - 1}(1 - \theta)^{b - 1} / B(a, b)} \\ & \color{steelblue}{\propto \theta^{a - 1}(1 - \theta)^{b - 1}} \end{align} \]
où \(a\) et \(b\) sont deux paramètres tels que \(a \geq 0\), \(b \geq 0\), et \(B(a, b)\) est une constante de normalisation.
Le niveau de certitude augmente avec la somme \(\kappa = a + b\).
Supposons que l’on dispose d’une estimation de la valeur la plus probable \(\omega\) du paramètre \(\theta\). On peut reparamétriser la distribution Beta en fonction du mode \(\omega\) et du niveau de certitude \(\kappa\) :
\[ \begin{align} a &= \omega(\kappa - 2) + 1 \\ b &= (1 - \omega)(\kappa - 2) + 1 &&\mbox{pour } \kappa > 2 \end{align} \]
Si \(\omega = 0.65\) et \(\kappa = 25\) alors \(p(\theta) = \mathrm{Beta}(\theta \given 15.95, 9.05)\).
Si \(\omega = 0.65\) et \(\kappa = 10\) alors \(p(\theta) = \mathrm{Beta}(\theta \given 6.2, 3.8)\).
Formellement, si \(\mathcal{F}\) est une classe de distributions d’échantillonnage \(p(y|\theta)\), et \(\mathcal{P}\) est une classe de distributions a priori pour \(\theta\), alors \(\mathcal{P}\) est conjuguée à \(\mathcal{F}\) si et seulement si :
\[ p(\theta \given y) \in \mathcal{P} \text{ for all } p(\cdot \given \theta) \in \mathcal{F} \text{ and } p(\cdot) \in \mathcal{P} \]
(p.35, Gelman et al., 2013). En d’autres termes, un prior est appelé conjugué si, lorsqu’il est converti en une distribution a posteriori en étant multiplié par la fonction de vraisemblance, il conserve la même forme. Dans notre cas, le prior Beta est un prior conjugué pour la vraisemblance binomiale, car le posterior est également une distribution Beta.
Le résultat du produit d’un prior Beta et d’une fonction de vraisemblance Binomiale est proportionnel à une distribution Beta. On dit que la distribution Beta est un prior conjugué de la fonction de vraisemblance Binomiale.
Soit un prior défini par : \(\ \color{steelblue}{p(\theta \given a, b) = \mathrm{Beta}(a, b) = \frac{\theta^{a - 1}(1 - \theta)^{b - 1}}{B(a, b)} \propto \theta^{a - 1}(1 - \theta)^{b - 1}}\)
Soit une fonction de vraisemblance associée à \(y\) “Face” pour \(n\) lancers : \(\ \color{orangered}{p(y \given n, \theta) = \mathrm{Bin}(y \given n, \theta) = \left(\begin{array}{l} n \\ y \end{array}\right) \theta^{y}(1 - \theta)^{n - y} \propto \theta^{y}(1 - \theta)^{n - y}}\)
Alors (en omettant les constantes de normalisation) :
\[ \begin{align} \color{purple}{p(\theta \given y, n)} &\propto \color{orangered}{p(y \given n, \theta)} \ \color{steelblue}{p(\theta)} &&\mbox{Théorème de Bayes} \\ &\propto \color{orangered}{\mathrm{Bin}(y \given n, \theta)} \ \color{steelblue}{\mathrm{Beta}(\theta \given a, b)} \\ &\propto \color{orangered}{\theta^{y}(1 - \theta)^{n - y}} \ \color{steelblue}{\theta^{a - 1}(1 - \theta)^{b - 1}} &&\mbox{Application des formules précédentes} \\ &\propto \color{purple}{\theta}^{\color{orangered}{y} + \color{steelblue}{a - 1}}\color{purple}{(1 - \theta)}^{\color{orangered}{n - y} + \color{steelblue}{b - 1}} &&\mbox{En regroupant les puissances des termes identiques} \\ \end{align} \]
Ici, on a ignoré les constantes qui ne dépendent pas de \(\theta\) (i.e., le nombre de combinaisons dans la fonction de vraisemblance binomiale et la fonction Beta \(B(a, b)\) dans le prior Beta).1 En les prenant en compte, on obtient en effet une distribution a posteriori Beta de la forme suivante :
\[\color{purple}{p(\theta \given y, n) = \mathrm{Beta}(y + a, n - y + b)}\]
On observe \(y = 7\) réponses correctes sur \(n = 10\) questions. On choisit un prior \(\mathrm{Beta}(1, 1)\), c’est à dire un prior uniforme sur \([0, 1]\). Ce prior équivaut à une connaissance a priori de 0 succès et 0 échecs (i.e., prior plat).
La distribution postérieure est donnée par :
\[ \begin{align} \color{purple}{p(\theta \given y, n)} &\propto \color{orangered}{p(y \given n, \theta)} \ \color{steelblue}{p(\theta)} \\ &\propto \color{orangered}{\mathrm{Bin}(7 \given 10, \theta)} \ \color{steelblue}{\mathrm{Beta}(\theta \given 1, 1)} \\ &= \color{purple}{\mathrm{Beta}(y + a, n - y + b)} \\ &= \color{purple}{\mathrm{Beta}(8, 4)} \end{align} \]
La moyenne de la distribution postérieure est donnée par :
\[ \color{purple}{\underbrace{\frac{y + a}{n + a + b}}_{posterior}} = \color{orangered}{\underbrace{\frac{y}{n}}_{data}} \underbrace{\frac{n}{n + a + b}}_{weight} + \color{steelblue}{\underbrace{\frac{a}{a + b}}_{prior}} \underbrace{\frac{a + b}{n + a + b}}_{weight} \]
Cas \(n < a + b, (n = 10, a = 4, b = 16)\).
Cas \(n = a + b, (n = 20, a = 4, b = 16)\).