TP Statistiques Bayésiennes: modèles hierarchiques et MCMC

library(tidyverse)

Objectifs

Ce TP a pour objectifs :

• Rappeler les équations analytiques du modèle Normal — Normal — Inverse-Gamma conjugué (moyenne inconnue et variance inconnue).
• Implémenter en R les formules analytiques des postérieurs (pour la moyenne et la variance) et visualiser ces distributions.
• Mettre en place un échantillonneur de Gibbs pour le même modèle et illustrer le fonctionnement des algorithmes MCMC.
• Comparer empiriquement les distributions analytiques et les distributions empiriques obtenues.

Modèle et rappels

Nous considérons les observations \(y_1, \dots, y_n\) modélisées par :

\[ y_i \mid \mu, \sigma^2 \overset{\text{iid}}{\sim} \mathcal{N}(\mu, \sigma^2), \]

et la hiérarchie de priors conjuguées suivante : \[ \mu \mid \sigma^2 \sim \mathcal{N}(\mu_0, \sigma^2 / \kappa_0),\qquad \sigma^2 \sim \text{Inv-Gamma}(\alpha_0, \beta_0), \]

où la paramétrisation de l’inverse-gamma est :

\[p(\sigma^2) \propto (\sigma^2)^{-\alpha_0-1} \exp\left(-\dfrac{\beta_0}{\sigma^2}\right).\]

Formules analytiques du postérieur

Avec les notations :

• \(n\) = nombre d’observations,
• \(\bar{y} = n^{-1} \sum_i y_i\),
• \(S = \sum_i (y_i - \bar{y})^2\),

les paramètres postérieurs conjugués sont :

\[ \kappa_n = \kappa_0 + n, \qquad \mu_n = \frac{\kappa_0 \mu_0 + n \bar{y}}{\kappa_0 + n}. \]

Pour la variance (paramètres de l’Inverse-Gamma a posteriori) :

\[ \alpha_n = \alpha_0 + \frac{n}{2}, \qquad \beta_n = \beta_0 + \frac{1}{2} S + \frac{\kappa_0 n}{2(\kappa_0 + n)} (\bar{y} - \mu_0)^2. \]

Les distributions a posteriori conjuguées sont alors :

• \(\sigma^2 \mid y \sim \text{Inv-Gamma}(\alpha_n, \beta_n)\),
• \(\mu \mid \sigma^2, y \sim \mathcal{N}(\mu_n, \sigma^2 / \kappa_n)\),
• Marginalement, \(\mu \mid y\) suit une Student-\(t\).

---

TP (format question / indice / réponse)

Exercice 1 — Simulation de données et calcul analytique des postérieurs

Question

Simulez un jeu de données \(n=30\) avec \(\mu_{true}=1.5\) et \(\sigma^2_{true}=2.0\). Choisissez les hyperparamètres \(\mu_0=0,\ \kappa_0=1,\ \alpha_0=2,\ \beta_0=1\) et calculez \(\bar{y}, S\) et les paramètres postérieurs \(\kappa_n,\mu_n,\alpha_n,\beta_n\). Affichez ces valeurs.

Indice

Calculez d’abord la moyenne empirique \(\bar{y}\) et \(S=\sum (y_i-\bar{y})^2\). Ensuite appliquez les formules du rappel.

Réponse

set.seed(2025)

n <- 30
mu_true <- 1.5
sigma2_true <- 2.0

y <- rnorm(n, mean = mu_true, sd = sqrt(sigma2_true))

# hyperparameters
mu0 <- 0
kappa0 <- 1
alpha0 <- 2
beta0 <- 1

# statistics
ybar <- mean(y)
S <- sum((y - ybar)^2)

# posteriors
kappa_n <- kappa0 + n
mu_n <- (kappa0 * mu0 + n * ybar) / kappa_n
alpha_n <- alpha0 + n/2
beta_n <- beta0 + 0.5 * S + (kappa0 * n) / (2 * kappa_n) * (ybar - mu0)^2

list(ybar = round(ybar,4), S = round(S,4), kappa_n = kappa_n, mu_n = round(mu_n,4),
     alpha_n = alpha_n, beta_n = round(beta_n,4))

$ybar
[1] 1.7924

$S
[1] 55.2818

$kappa_n
[1] 31

$mu_n
[1] 1.7346

$alpha_n
[1] 17

$beta_n
[1] 30.1955

Exercice 2 — Visualisation des postérieurs analytiques

Question

Tracez les graphiques suivants :

• la densité a posteriori de \(\sigma^2\) (Inv-Gamma avec \(\alpha_n,\beta_n\)),
• la densité conditionnelle \(\mu\mid\sigma^2\) pour une valeur représentative de \(\sigma^2\) (par ex. la moyenne postérieure)

Indice

Utilisez la densité de l’Inverse-Gamma et simulez \((\sigma^2,\mu)\) en tirant \(\sigma^2\sim\text{Inv-Gamma}(\alpha_n,\beta_n)\) puis \(\mu\mid\sigma^2\sim\mathcal{N}(\mu_n,\sigma^2/\kappa_n)\).

Réponse

# grid for sigma2
sigma2_grid <- seq(0.01, 6, length.out = 400)
# Inverse-Gamma density function
inv_gamma_density <- function(x, shape, scale) {
  beta <- scale
  alpha <- shape
  (beta^alpha) / gamma(alpha) * x^(-alpha - 1) * exp(-beta / x)
}

# posterior density of sigma2
p_sigma2 <- data.frame(sigma2 = sigma2_grid,
                       density = inv_gamma_density(sigma2_grid, alpha_n, beta_n))

# conditional mu | sigma2 at sigma2 = posterior mean of sigma2 (approx: beta_n/(alpha_n-1))
sigma2_rep <- beta_n / (alpha_n - 1)
mu_sd_cond <- sqrt(sigma2_rep / kappa_n)

p_mu_cond <- data.frame(mu = seq(mu_n - 3*mu_sd_cond, mu_n + 3*mu_sd_cond, length.out = 300))
p_mu_cond$density <- dnorm(p_mu_cond$mu, mean = mu_n, sd = mu_sd_cond)

# plots
p1 <- ggplot(p_sigma2, aes(x=sigma2, y=density)) + geom_line(color='darkred') +
  labs(title='Posterior of sigma^2 (Inv-Gamma)', x=expression(sigma^2), y='densite') + theme_minimal()

p2 <- ggplot(p_mu_cond, aes(x=mu, y=density)) + geom_line(color='darkblue') +
  labs(title=expression('Conditional posterior ' ~ mu ~ '| ' ~ sigma^2), x=expression(mu), y='densite') + theme_minimal()

p1

p2

Implémentation d’un Gibbs « à la main »

Aide pour l’implémentation d’un échantillonneur de Gibbs

Pour implémenter un échantillonneur de Gibbs, il faut connaître les distributions conditionnelles suivantes:

• \(p(\sigma^2 \mid \mu, y)\) (paramètres de l’Inverse-Gamma a posteriori) ;
• \(p(\mu \mid \sigma^2, y)\) (moyenne et variance conditionnelles).

On rappelle que ces distributions sont :

• \(p(\sigma^2 \mid \mu, y) \sim \text{Inv-Gamma}\big(\alpha_0 + n/2,\; \beta_0 + \tfrac{1}{2}\sum_i (y_i-\mu)^2\big)\).
• \(p(\mu \mid \sigma^2, y) \sim \mathcal{N}\big(\mu_n,\; \sigma^2 / \kappa_n\big)\), avec \(\kappa_n = \kappa_0 + n\) et \(\mu_n = (\kappa_0\mu_0 + n\bar y)/(\kappa_0 + n)\).

Ces formules serviront directement pour tirer les mises à jour dans un Gibbs exact.

Exercice 3 — Implémenter une mise à jour

Question

Implémentez en R une fonction qui effectue une mise à jour unique (un pas de Gibbs) :

Indice

Utilisez rgamma() / 1/rgamma() pour tirer l’Inverse-Gamma, et rnorm() pour tirer la normale conditionnelle.

Réponse

# inputs: mu_curr, sigma2_curr, y, hyperparametres
# outputs: mu_new, sigma2_new

gibbs_step <- function(mu_curr, sigma2_curr, y, mu0, kappa0, alpha0, beta0) {
n <- length(y)

# sigma2 | mu, y  ~ Inv-Gamma(alpha_post, beta_post)

alpha_post <- alpha0 + n/2
S_mu <- sum((y - mu_curr)^2)
beta_post <- beta0 + 0.5 * S_mu
sigma2_new <- 1 / rgamma(1, shape = alpha_post, rate = beta_post)

# mu | sigma2, y ~ Normal(mu_n, sigma2 / kappa_n)

kappa_n <- kappa0 + n
mu_n <- (kappa0 * mu0 + n * mean(y)) / kappa_n
sd_mu <- sqrt(sigma2_new / kappa_n)
mu_new <- rnorm(1, mean = mu_n, sd = sd_mu)

return(list(mu = mu_new, sigma2 = sigma2_new))
}

# quick test

set.seed(1)
res_test <- gibbs_step(mu_curr = mean(y), sigma2_curr = var(y), y = y,
mu0 = mu0, kappa0 = kappa0, alpha0 = alpha0, beta0 = beta0)
res_test

$mu
[1] 2.075588

$sigma2
[1] 2.037998

Exercice 4 — Implémenter la boucle Gibbs complète

Question

À partir de la fonction précédente, codez la boucle Gibbs complète qui produit niter itérations, applique un burn-in et retourne des échantillons conservés pour mu et sigma2. Affichez des diagnostics (traceplots, densités) et comparez sommairement aux valeurs analytiques (moyennes postérieures).

Indice

Stockez les chaînes dans des vecteurs numériques et utilisez plot() / ggplot2 pour les diagnostics. Pensez à fixer la graine pour reproductibilité.

Réponse

set.seed(42)

niter <- 6000
nburn <- 1000
mu_chain <- numeric(niter)
sigma2_chain <- numeric(niter)


# initialisation

for(t in 2:niter) {
  step <- gibbs_step(mu_chain[t-1], sigma2_chain[t-1], y, mu0, kappa0, alpha0, beta0)

  mu_chain[t] <- step$mu
  mu_chain[1] <- mean(y)

  sigma2_chain[1] <- var(y)
  sigma2_chain[t] <- step$sigma2
}

# extract post-burn-in

mu_samps_gibbs <- mu_chain[(nburn+1):niter]
sigma2_samps_gibbs <- sigma2_chain[(nburn+1):niter]

# diagnostics

plot(mu_samps_gibbs, type='l', main='Traceplot mu (Gibbs)', ylab='mu')

plot(sigma2_samps_gibbs, type='l', main='Traceplot sigma2 (Gibbs)', ylab='sigma2')

dens_mu = tibble(x = density(mu_samps_gibbs)$x, y = density(mu_samps_gibbs)$y)
dens_sigma = tibble(x = density(sigma2_samps_gibbs)$x, y = density(sigma2_samps_gibbs)$y)

# comparison with analytical posteriors

p1 + geom_line(data = dens_sigma, aes(x=x, y=y), color='red', linetype='dashed')

p2 + geom_line(data = dens_mu, aes(x=x, y=y), color='blue', linetype='dashed')

# summary statistics

mean_mu_gibbs <- mean(mu_samps_gibbs)
mean_sigma2_gibbs <- mean(sigma2_samps_gibbs)

list(mean_mu_gibbs = mean_mu_gibbs,
mean_mu_analytical = mu_n,
mean_sigma2_gibbs = mean_sigma2_gibbs,
mean_sigma2_analytical = beta_n / (alpha_n - 1))

$mean_mu_gibbs
[1] 1.729404

$mean_mu_analytical
[1] 1.734625

$mean_sigma2_gibbs
[1] 1.855577

$mean_sigma2_analytical
[1] 1.887219