Initiation aux statistiques Bayésiennes

L’objectif de ce document est de proposer une initiation aux stats Bayésiennes, dans le contexte des essais cliniques, et auprès d’un public ayant déjà un certain bagage de statistiques fréquentistes.

On souhaite également ne pas rester dans l’abstraction, et proposer des petits exemples applicatifs que chacun puisse répliquer/modifier.

Pour les utilisateurs de R, on commence par charger un certain nombre de modules intéressants.

# le package moderne standard pour manipuler ses données
library(tidyverse)
# pour faire des graphes multiples
library(cowplot)
# pour afficher joliment les tableaux dans un document rmarkdown
library(knitr)
# pour la loi beta-binomiale dont je parle dans le doc
library(extraDistr)

# quelques couleurs manuelles pour intégrer les figures à la présentation associée
bleufonce <- "#3d5468"
bleuclair <- "#5b7c98"
rose <- "#ff5555"

Exemple introductif fréquentiste

Données binaires

Les stats Bayésiennes reposent sur un vocabulaire précis, qu’on va essayer d’aborder pas à pas autour d’un exemple minimaliste.

Cet exemple est le suivant : imaginez que vous étudiez un échantillon de \(n\) patients, classés en “répondeur” ou “non-répondeur”. La donnée observée est donc très simple : il s’agit de \(n_r\), le nombre de patients répondeurs.

On fixe la valeur de \(n\) et le nombre observé \(n_r\) de patients répondeurs à l’issue de l’étude, pour tout le document.

n <- 120
n_r <- 90

Modèle fréquentiste et vraisemblance

Si vous pensez comme un fréquentiste, vous créez probablement dans votre esprit le modèle suivant : il existe une probabilité \(p\) fixée d’être répondeur, et une probabilité \(1-p\) d’être non-répondeur.

Chaque individu est indépendant de tous les autres, et le nombre \(N_r\) de répondeur suit donc une loi binomiale :

\[ N_r \sim \mathcal{B}(n, p) \]

Si vous fixez \(p\), vous êtes capable de faire des simulations sous votre modèle et de représenter la distribution attendue pour vos données. La figure suivante nous illustre ça si \(p = 0.6\).

p <- 0.6
d <- tibble(x = 0:n) %>%
    mutate(density = dbinom(x, size=n, prob=p))
d %>% ggplot(aes(x = x, y = density)) +
    geom_bar(stat="identity", fill=bleuclair) +
    geom_vline(xintercept = n_r, col=rose) +
    xlab("Nombre de répondeurs") +
    ylab("Probabilité si p vaut 0.6") +
    theme_bw()

Estimation en fréquentiste

Si votre problème est un problème d’estimation du paramètre \(p\), vous allez typiquement chercher le paramètre \(\hat{p}\) maximisant la vraisemblance, c’est à dire :

\[ \hat{p} = \text{argmax}_p ~\mathbb{P}( N_r = n_r | p ) \]

Dans ce cas très précis, on pourrait montrer que ça donne exactement l’estimation qu’on aurait eu envie de faire intuitivement : \(\hat{p} = N_r / n\) ce qui donne, dans cet exemple, une estimation ponctuelle à 0.75. Visuellement, on a donc choisi \(p\) tel que la vraisemblance ressemble à ça :

p <- round(n_r/n,2)
d <- tibble(x = 0:n) %>%
    mutate(density = dbinom(x, size=n, prob=p))
d %>% ggplot(aes(x = x, y = density)) +
    geom_bar(stat="identity", fill=bleuclair) +
    geom_vline(xintercept = n_r, col=rose) +
    xlab("Nombre de répondeurs") +
    ylab(paste("Probabilité si p vaut ", p, sep="")) +
    theme_bw()

Dans un problème d’estimation classique, on aimerait en plus obtenir un intervalle de confiance pour \(p\), c’est à dire deux variables aléatoires \(\hat{a}\) et \(\hat{b}\) vérifiant :

\[ \mathbb{P}(p \in (\hat{a},\hat{b})) \geq 95\% \]

Notons que ce type de résultat mathématique n’est pas en soi si facile à comprendre. Ici, ce sont les bornes de l’intervalle qui sont des variables aléatoires. Celles qui seront finalement calculées ne serot plus les estimateurs (aléatoires), mais leurs estimations. La seule garantie donnée par ce résultat, c’est que si on répétait l’expérience un grand nombre de fois, 95% des intervalles qu’on donnerait contiendrait la vraie valeur de \(p\).

C’est de plus (en général) un problème difficile ! Même pour ce cas très simple, ça nécessite de reposer sur des théorèmes asymptotiques, et de produire des “beaux” résultats de maths.

Pour tous les modèles relativement simples / classiques, tout ça a déjà été fait pour nous, et il ne nous reste qu’à faire confiance aux outils implémentés pour obtenir l’estimation de notre intervalle :

ci <- prop.test(x=n_r, n=n)$conf.int

Ici l’intervalle de confiance à 95% englobe les valeurs [0.66, 0.82].

Modèle introductif version Bayésien

Vraisemblance identique

Si vous pensez comme un Bayésien, vous commencez par vous dire la même chose que précédemment : si vous connaissiez \(p\), la probabilité d’être répondeur, alors vous sauriez que le nombre de répondeur est :

\[ N_r | p ~~ \sim ~~ \mathcal{B}(n, p) \]

La loi de \(N_r\) (les données) sachant \(p\) (le paramètre) est appelée vraisemblance et est notée \(\mathbb{P}(N_r | p)\).

Choix d’un prior

Seulement voilà : le raisonnement Bayésien suppose qu’il n’y a pas une unique valeur de \(p\), mais qu’il y a en fait une distribution de probabilité sur le paramètre lui-même ! C’est ce qu’on appelle prior : \(\mathbb{P}(p)\).

Puisque \(p\) vit dans (0,1), il faut l’équiper d’une loi de probabilité sur (0,1), par exemple une loi uniforme.

tibble(x = seq(-0.1, 1.1, 0.001)) %>%
    mutate(y = dunif(x)) %>%
    ggplot(aes(x=x, y=y)) +
    geom_line(col = bleuclair) +
    xlab("Valeur de p") +
    ylab("densité") +
    theme_bw()

Marginale des observations

Dans un monde Bayésien, ce n’est que maintenant qu’on a choisi une vraisemblance et un prior qu’on peut simuler une réalisation de \(N_r\) sous notre modèle. On simule une réalisation de \(p\) tiré dans le prior, puis une réalisation de loi binomiale ayant cette valeur de \(p\) :

p_simu = runif(1, min=0, max=1)
n_r_simu = rbinom(1, size=n, prob=p_simu)

Ca revient à tirer une valeur dans la loi de \(N_r\), qu’on appelle loi marginale des observations. Cette loi peut s’exprimer comme :

\[ \mathbb{P}(N_r) = \int_0^1 \mathbb{P}(N_r | p) \mathbb{P}(p \in dp) dp \]

Il est très facile de simuler dans cette loi, mais elle n’est pas toujours facile à caractériser plus finement. Dans notre scénario, on peut toujours représenter l’allure de la distribution obtenue par simulation :

1. on simule une valeur de \(p\) dans le prior de \(p\),
2. on simule une valeur de \(N_r\) dans la distribution de \(N_r\) sachant que \(p\) est fixé à la valeur simulée à la première étape.

sim_n_r <- function(n, sample_size){
    p_simu = runif(sample_size, min=0, max=1)
    return (rbinom(sample_size, size=n, prob=p_simu))
}
tibble(x = sim_n_r(n, 1e5)) %>% 
    group_by(x) %>%
    summarize(n = n()/1e5) %>%
    ggplot(aes(x = x, y = n)) +
    geom_bar(stat="identity", fill=bleuclair) +
    theme_bw() +
    xlab("Nombre de répondeurs") +
    ylab("Probabilité")

Posterior

C’est LA quantité d’intérêt pour le statisticien ! Dans la réalité, on a observé UNE valeur \(n_r\). On cherche donc :

\[ \mathbb{P}(p | N_r = n_r) = \frac{ \mathbb{P}(N_r = n_r | p) \mathbb{P}(p) }{ \mathbb{P}(N_r = n_r) } \]

On obtient le posterior en combinant l’idée qu’on avait a priori sur \(p\) (le prior) avec ce qu’on apprend de notre expérience (la vraisemblance).

J’insiste également une nouvelle fois sur ce changement de paradigme majeur : on cherche une distribution de probabilité sur \(p\), et on arrête de penser qu’il y a une valeur unique de \(p\).

Plusieurs possibilités peuvent se présenter pour le calcul du posterior. Mais toutes utilisent déjà une règle commune : on va faire abstraction du dénominateur (la marginale de \(N_r\)) qui était pénible à calculer et qui ne sert ici que de facteur de normalisation. On commence donc toujours par ré-écrire l’équation précédente comme :

\[ \mathbb{P}(p | N_r = n_r) \propto \mathbb{P}(N_r = n_r | p) \mathbb{P}(p) \]

Puis, dans ce cas précis, tout est très simple car le prior de \(p\) est uniforme, ce qui fait qu’on peut ré-écrire notre équation comme :

\[\begin{align*} \mathbb{P}(p | N_r = n_r) &\propto \binom{90}{120} p^{90} (1-p)^{120 - 90} \\ &\propto p^{90} (1-p)^{120 - 90} \end{align*}\]

Où on reconnaît une loi Bêta de paramètres \(\alpha = 90+1\) et \(\beta = 120-90+1\). Le posterior a l’allure suivante :

alpha <- n_r + 1
beta <- n - n_r + 1
d_posterior <- tibble(x = seq(0, 1, 0.001)) %>%
    mutate(y = dbeta(x, shape1=alpha, shape2=beta))
plot_posterior <- d_posterior %>%
    ggplot(aes(x=x, y=y)) +
    geom_line(col = bleuclair) +
    xlab("Valeur de p") +
    ylab("densité du posterior") +
    theme_bw()
plot_posterior

Si on souhaite un estimateur ponctuel, il est d’usage d’utiliser le maximum a posteriori, c’est à dire la valeur de \(p\) qui maximise le posterior. Ici, le calcul nous donne exactement la même solution que dans le cas fréquentiste.

plot_posterior +
    geom_vline(xintercept = 90/120, col=rose)

Intervalle de crédibilité

Si on accepte qu’il n’existe pas une unique valeur du paramètre \(p\), l’intervalle de crédibilité du Bayésien est beaucoup plus intuitif que l’intervalle de confiance du fréquentiste.

Il s’agit en effet d’un intervalle \((a,b)\) tel que :

\[ \mathbb{P}(p \in (a,b) | N_r = 90) \geq 1-\alpha \]

Typiquement, si on souhaitait avoir un intervalle de crédibilité à 95%, on prendrait les valeurs, sur le support du posterior, contenues entre les quantiles de niveau 0.025 et 0.975.

low <- qbeta(0.025, shape1=alpha, shape2=beta)
up <- qbeta(0.975, shape1=alpha, shape2=beta)
plot_posterior +
    geom_vline(xintercept = c(low, up), col=rose)

Ici, l’intervalle de crédibilité à 95% est [0.67, 0.82].

Test d’hypothèse

Enfin, les tests d’hypothèse sont réalisés de façon bien plus intuitive que dans le cas fréquentiste.

1. on calcule la probabilité a posteriori de chaque hypothèse,
2. on base le choix sur ces deux probabilités.

Par exemple, si on souhaite décider entre :

• H0 : \(p \leq 0.7\)
• H1 : \(p > 0.7\),

Alors il nous suffit de s’accorder sur un threshold de probabilité a posteriori à partir duquel on choisit H1, et de calculer nos probabilités. Visuellement, on se base donc sur :

plot_posterior +
    geom_area(data = d_posterior %>% filter(x <= 0.7), fill=rose, alpha=0.5)

Les MCMCs, la cuisine du Bayésien

Tout ce qu’on a effectué ci-dessus était réalisé dans un cas simple idéalisé pour lequel on arrivait à faire des calculs (ou plutôt, raisonner quasi sans calculs) très facilement. Ce n’est malheureusement pas la majorité des cas.

Un prior plus complexe

Dans beaucoup de situations, le prior n’est pas aussi simple, et le posterior est donc plus difficile à caractériser. Modifions notre exemple précédent : plutôt que d’avoir \(p\) uniforme sur [0,1], des données préliminaires issues de la littérature, ou des avis d’experts du domaine, suggèrent que \(p\) serait centré sur la valeur 0.6, avec un écart-type de 0.3.

Sans trop réfléchir, on se propose donc de voir ce que ça donnerait avec le prior suivant : une loi normale de moyenne 0.6 et écart-type 0.3. Seulement voilà, la cuisine commence : la loi normale est définie sur \(\mathbb{R}\). Qu’à cela ne tienne : considérons qu’on tronque cette distribution pour qu’elle ait un support sur [0,1].

mu <- 0.6
sd <- 0.3
weight_inf_0 <- pnorm(0, mean=mu, sd=sd)
weight_sup_1 <- 1 - pnorm(1, mean=mu, sd=sd)
tibble(x = seq(0, 1, 0.001)) %>%
    mutate(y = dnorm(x, mean=mu, sd=sd)/(1 - weight_inf_0 - weight_sup_1)) %>%
    ggplot(aes(x=x, y=y)) +
    geom_line(col = bleuclair) +
    xlab("Valeur de p") +
    ylab("densité du prior") +
    theme_bw()

La marginale des observations

Comme précédemment, il n’est pas évident de caractériser la loi de \(N_r\), mais on peut l’approcher par simulations :

sim_n_r <- function(n, sample_size){
    res <- c()
    for (i in 1:sample_size){
        p_simu = rnorm(1, mean=mu, sd=sd)
        while (p_simu < 0 | p_simu > 1){
            p_simu = rnorm(1, mean=mu, sd=sd)
        }
        n_r_simu <- rbinom(1, size=n, prob=p_simu)
        res <- c(res, n_r_simu)
    }
    return (res)
}
tibble(x = sim_n_r(n, 1e5)) %>% 
    group_by(x) %>%
    summarize(n = n()/1e5) %>%
    ggplot(aes(x = x, y = n)) +
    geom_bar(stat="identity", fill=bleuclair) +
    theme_bw() +
    xlab("Nombre de répondeurs") +
    ylab("Probabilité")

Le posterior

Comme précédemment, on peut écrire la règle de proportionnalité de base suivante :

\[ \mathbb{P}(p | N_r = n_r) \propto \mathbb{P}(N_r = n_r | p) \mathbb{P}(p) \]

Puis, en remplaçant avec ce qu’on connaît de la vraisemblance et du prior :

\[\begin{align*} \mathbb{P}(p | N_r = n_r) &= 0 ~~~\text{ si } p \notin [0,1]\\ &\propto p^{90} (1-p)^{120 - 90} e^{-\frac{1}{2} \left( \frac{p - 0.6}{0.3} \right)^2} ~~~\text{ sinon.} \end{align*}\]

Cette fois, impossible de reconnaître quoi que ce soit ! Dans ce cas de figure, qui est d’ailleurs la norme dans les applications de stats Bayésiennes, on doit s’appuyer sur des outils de simulations numériques pour échantillonner des valeurs dans le posterior sans arriver à le caractériser plus finement que ça.

Et l’outil central pour réaliser ces simulations s’appelle le MCMC, pour Markov Chain Monte Carlo.

Principe d’un MCMC

L’idée générale des techniques MCMC consiste à construire une chaîne de Markov, une suite de variables aléatoires bien particulière, qui converge vers notre distribution d’intérêt : le posterior qu’on note ici \(\nu\).

La façon la plus classique/facile de construire une telle chaîne consiste à utiliser l’algorithme de Metropolis-Hastings. Celui-ci repose sur le choix adhoc d’une loi de proposition de movement \(q\), avec \(\forall x,y, q(x, y)\) qui donne la probabilité de proposer l’état \(y\) sachant qu’on est dans l’état \(x\).

A deux-trois conditions techniques près, l’algorithme consiste à passer de l’état \(x_i\) au pas \(i\), à \(x_{i+1}\) au pas \(i+1\) en :

1. proposant un état \(y\) en tirant dans la loi de \(q(x_i,.)\),
2. calculant la quantité suivante : \[ r = \frac{\nu(y)q(y,x_i)}{\nu(x_i)q(x_i,y)} \]
3. Si \(r \geq 1\), \(x_{i+1} = y\). Sinon \(x_{i+1} = y\) avec probabilité \(r\) et \(x_{i+1} = x_i\) sinon.

Il ne reste ensuite qu’à :

• supprimer les premières itérations de la chaîne (pour réduire la dépendance aux conditions initiales),
• échantillonner une valeur sur 10, ou sur 100, ou … (pour réduire l’auto-corrélation des valeurs).

Pour se prouver que c’est très simple à designer, on peut en construire une ici. Par simplicité, prenons une proposition de mouvement \(q\) très simple : la loi uniforme sur [0,1], quelle que soit la valeur de départ.

# notre distribution cible (posterior) est proportionnelle à nu
nu <- function(p){
    res <- p^n_r * (1-p)^(n-n_r) * exp(-0.5 * ((p - mu)/sd)^2 )
    return (res)
}

# on initialise une chaîne vide qu'on remplit un pas après l'autre
# attention, ce morceau de code peut prendre un certain temps à tourner
chain <- c()
nsteps <- 1e5
p <- 0.2
for (i in 1:nsteps){
    chain <- c(chain, p)
    y <- runif(1, min=0, max=1)
    r <- nu(y) / nu(p)
    if (r >= 1){
        p <- y
    }else{
        u = runif(1, min=0, max=1)
        if (u <= r){
            p <- y
        }
    }
}

Tout l’art de la cuisine de MCMC (sauce MH) consiste à choisir empiriquement une proposition de mouvement \(q\) qui “fonctionne” plutôt bien.

Il faut aussi choisir un peu au pif combien de valeurs on supprime au début de la chaîne (burn-in) et avec quelle fréquence on échantillonne. Pour ça, on peut regarder l’allure de la chaîne :

dmcmc <- tibble(p = chain,
                step = seq(1, nsteps)) 
dmcmc %>%
    ggplot(aes(x=step, y=p)) +
    geom_line(col = bleufonce) +
    theme_bw()

Sur cet exemple tout simple, le burn-in est facile à repérer visuellement. Si on enlève les \(10^4\) premières valeurs et qu’on conserve une valeur sur 10, on conserve assez de valeurs pour estimer notre posterior à partir de l’histogramme empirique suivant :

dmcmc %>%
    filter(step > 1e4 & step %% 10 == 0) %>%
    ggplot(aes(x=p)) +
    geom_histogram(aes(y=..density..), bins=50, col=bleuclair, fill=bleuclair, alpha=0.5) +
    geom_density(col=rose, alpha=0.5) +
    ylab("densité de probabilité") +
    theme_bw()

On peut finalement obtenir un intervalle de crédibilité à 95% en récupérant les quantiles empiriques de notre posterior :

ic95_mcmc <- dmcmc %>%
    filter(step > 1e4 & step %% 10 == 0) %>%
    pull(p) %>%
    quantile(probs = c(0.025, 0.975))

Ici, l’intervalle de crédibilité à 95% est [0.66, 0.82].

Le prior conjugué, un outil plus smart

Principe mathématique

On l’a vu, le MCMC permet de faire toute la cuisine qu’on souhaite, en mode Bayésien. Ce n’est pas encore visible avec cet exemple simple qui ne s’intéresse qu’à une seule variable \(p\), mais ça peut devenir très vite lourd d’un point de vue computationnel quand on a un espace de paramètres en plus grande dimension à explorer. En plus de ça, c’est très largement “inélégant” d’un point de vue maths, et on se retrouve souvent à utiliser tout et n’importe quoi comme prior, même des distributions absolument pas faites, à la base, pour avoir le support qu’on leur donne.

Une alternative plus élégante, mais qui nécessite de réfléchir un peu plus en amont, consiste à s’appuyer sur un modèle plus approprié, avec, par exemple, un prior conjugué. Sans entrer dans les détails, pour l’essentiel des distributions de vraisemblance que vous pourriez connaître (qui font partie d’une grande classe de distributions appelée famille exponentielle), il existe ce qu’on appelle des priors/posteriors conjugués, qui ont le bon goût d’avoir la même forme de distribution, avec des paramètres différents.

Prior conjugué sur notre exemple

Sur notre exemple, la vraisemblance est une loi binomiale. Le prior conjugué correspondant est une loi Bêta, avec la propriété suivante :

\[\begin{align*} &p \sim \mathcal{Beta}(\alpha, \beta) \text{ et } N_r | p \sim \mathcal{B}(n, p) \\ \Longrightarrow ~&p | N_r \sim \mathcal{Beta}(\alpha + N_r, \beta + n - N_r) \end{align*}\]

Ce genre de propriété est de plus très facile à démontrer, puisque :

1. la distribution Beta est telle que \(\mathbb{P}(p) \propto p^{\alpha-1} (1-p)^{\beta-1}\).
2. la distribution Binomiale de la vraisemblance nous donne \(\mathbb{P}(N_r | p) \propto p^{N_r} (1-p)^{n-N_r}\).
3. et donc le posterior s’obtient en écrivant : \[\begin{align*} \mathbb{P}(p | N_r) &\propto \mathbb{P}(p) \mathbb{P}(N_r | p) \\ &\propto p^{alpha + N_r - 1} (1-p)^{\beta + n - N_r - 1} \end{align*}\]

Où on reconnaît, donc, une loi Bêta de paramètres \((\alpha + N_r, \beta + (n - N_r))\). C’est extrêmement pratique, puisque nos hyperparamètres ont maintenant une interprétation intuitive en terme de nombre de répondeurs/non-répondeurs qu’on aurait déjà vu avant de faire notre nouvelle expérience.

Par exemple, on pourrait regarder la biblio, se dire qu’une précédente étude a déjà regardé notre problème d’intérêt. Celle-ci avait enregistré \(\alpha = 50\) répondeurs et \(\beta = 60\) non-répondeurs. On va donc partir sur un prior avec loi Beta de paramètres \((\alpha=50, \beta=60)\).

alpha <- 50
beta <- 60
tibble(x = seq(0, 1, 0.001)) %>%
    mutate(y = dbeta(x, shape1=alpha, shape2=beta)) %>%
    ggplot(aes(x=x, y=y)) +
    geom_line(col = bleuclair) +
    xlab("Valeur de p") +
    ylab("densité du prior") +
    theme_bw()

Loi de la marginale

Dans tous ces cas avec un prior conjugué, la loi de la marginale des observations est aussi connue avec des formules fermées, bien qu’en général pas triviale. Ici, la loi de la marginale est appelée Bêta-binomiale, et son allure est affichée ci-dessous :

tibble(k = seq(0, n, 1)) %>%
    mutate(y = dbbinom(k, size=n, alpha=alpha, beta=beta)) %>%
    ggplot(aes(x=k, y=y)) +
    geom_line(col = bleuclair) +
    xlab("Nombre de répondeurs") +
    ylab("probabilité") +
    theme_bw()

Loi du posterior

Puisqu’on utilise un prior conjugué, le posterior nous est donné sur un plateau d’argent : il s’agit d’une loi Bêta de paramètres \(\alpha + N_r\), \(\beta + n - N_r\).

tibble(x = rep(seq(0, 1, length.out=500), 2),
       distribution = rep(c("prior", "posterior"), each=500)) %>%
    mutate(y = dbeta(x, 
                     shape1=rep(c(alpha, alpha + n_r), each=500), 
                     shape2=rep(c(beta, beta + n - n_r), each=500)) ) %>%
    ggplot(aes(x=x, y=y, col=distribution)) +
    geom_line() +
    xlab("Valeur de p") +
    ylab("densité") +
    theme_bw()

Sur cet exemple, on peut une nouvelle fois calculer un intervalle de crédibilité à 95%,

ic95_conj <- qbeta(c(0.025, 0.975), shape1=alpha+n_r, shape2=beta+n-n_r)

Ici, l’intervalle de crédibilité à 95% est [0.54, 0.67].

On remarque ici une vraie différence de posterior, et donc d’intervalle de crédibilité, par rapport aux exemples précédents. Ceci est entièrement dû à l’influence du prior, qui est pris en compte dans notre posterior et le “shifte” sur la gauche par rapport aux exemples précédents.

Conclusion

Je tente le tableau suivant pour résumer les avantages et inconvénients du fréquentiste et du Bayésien pour les essais cliniques, basés sur ce qui est décrit dans ce notebook :

Les autres points importants sont les suivants :

• Pour les problèmes d’estimation, la notion d’intervalle de crédibilité est plus facile à comprendre que l’intervalle de confiance.
• Pour les tests d’hypothèses, il est possible de faire des choix, en Bayésien, basés directement sur la probabilité des deux hypothèses sachant les données, ce qui est aussi notablement plus facile à comprendre.
• La notion de non-informative prior est assez largement mensongère, à remplacer par celle de flat prior.
• Les priors conjugués sont très intéressants dans les cas simples pour avoir un modèle simple et élégant.
• Dans tous les autres cas, des MCMC permettent de déterminer la distribution à posteriori.