Accueil › Glossaire › Gibbs Sampling

Gibbs Sampling (Échantillonnage de Gibbs)

Le Gibbs sampling est un algorithme de Monte Carlo par chaînes de Markov (MCMC) qui échantillonne une distribution jointe multidimensionnelle en mettant à jour chaque variable tour à tour, conditionnellement à toutes les autres. C’est l’un des algorithmes les plus utilisés en inférence bayésienne, notamment pour les modèles hiérarchiques et le topic modeling.

L’idée est d’une élégance remarquable. Vous ne savez pas échantillonner la distribution jointe P(x₁, x₂, …, xₙ) directement, mais vous connaissez chaque distribution conditionnelle P(xᵢ | x₁, …, xᵢ₋₁, xᵢ₊₁, …, xₙ). Alors parcourez les variables une par une, échantillonnez chacune conditionnellement aux valeurs actuelles des autres, et répétez. La séquence d’échantillons ainsi produite forme une chaîne de Markov dont la distribution stationnaire est exactement la distribution jointe recherchée.

Le Gibbs sampling a été nommé d’après le physicien Josiah Willard Gibbs (1839-1903), en référence aux distributions de Gibbs en mécanique statistique. L’algorithme a été introduit en traitement d’images par Geman et Geman (1984), puis popularisé en statistique bayésienne par Gelfand et Smith (1990), déclenchant la « révolution MCMC » des années 1990 qui a transformé la pratique de l’inférence statistique.

L’algorithme pas à pas

Le principe

Soit une distribution jointe P(x₁, x₂, …, xd) sur d variables dont on veut échantillonner. Le Gibbs sampling procède ainsi :

1. Initialiser (x₁⁰, x₂⁰, …, xd⁰) à des valeurs arbitraires.

2. À chaque itération t, pour chaque variable i de 1 à d : échantillonner xᵢᵗ depuis la distribution conditionnelle P(xᵢ | x₁ᵗ, …, xᵢ₋₁ᵗ, xᵢ₊₁ᵗ⁻¹, …, xdᵗ⁻¹), en utilisant les valeurs les plus récentes de toutes les autres variables.

3. Répéter l’étape 2 pendant N itérations.

4. Jeter les premières itérations (burn-in) et utiliser le reste comme échantillons de la distribution jointe.

Le point clé est qu’à chaque mise à jour, on utilise la valeur la plus récente de chaque variable. Si on vient de mettre à jour x₁ et x₂, alors pour mettre à jour x₃, on utilise les nouvelles valeurs de x₁ et x₂, pas les anciennes.

Exemple : gaussienne bivariée

Considérons une distribution normale bivariée avec corrélation ρ :

(X, Y) ~ N(0, Σ) avec Σ = [[1, ρ], [ρ, 1]]

Les distributions conditionnelles sont connues analytiquement : X|Y=y ~ N(ρy, 1-ρ²) et Y|X=x ~ N(ρx, 1-ρ²). Le Gibbs sampler alterne entre ces deux tirages :

import numpy as np

def gibbs_bivariate_normal(rho, n_samples=10000, burn_in=1000):
    """Gibbs sampling pour une gaussienne bivariée corrélée."""
    samples = np.zeros((n_samples + burn_in, 2))
    x, y = 0.0, 0.0  # Initialisation arbitraire

    for t in range(n_samples + burn_in):
        # Échantillonner X | Y = y
        x = np.random.normal(rho * y, np.sqrt(1 - rho**2))
        # Échantillonner Y | X = x (utilise le NOUVEAU x)
        y = np.random.normal(rho * x, np.sqrt(1 - rho**2))
        samples[t] = [x, y]

    return samples[burn_in:]  # Jeter le burn-in

# Exécuter
samples = gibbs_bivariate_normal(rho=0.9, n_samples=5000)
print(f"Corrélation empirique : {np.corrcoef(samples.T)[0,1]:.3f}")
# Devrait être proche de 0.9

Visuellement, les échantillons se déplacent le long des axes coordonnés : d’abord horizontalement (mise à jour de X), puis verticalement (mise à jour de Y). Après suffisamment d’itérations, les échantillons couvrent la distribution jointe. Quand ρ est élevé (forte corrélation), les déplacements le long des axes sont petits par rapport à l’étendue de la distribution, et la chaîne explore lentement. C’est la principale faiblesse du Gibbs sampling.

Quand utiliser le Gibbs sampling

Le prérequis fondamental

Le Gibbs sampling nécessite de pouvoir échantillonner les distributions conditionnelles complètes (full conditionals). Cette condition est satisfaite naturellement dans deux cas importants :

Modèles conjugués : quand le prior et la vraisemblance sont conjugués (ex. : prior gaussien + vraisemblance gaussienne, prior Beta + vraisemblance Bernoulli, prior Dirichlet + vraisemblance catégorielle), les conditionnelles complètes sont des distributions connues dont on sait échantillonner directement.

Modèles graphiques : dans les modèles graphiques probabilistes, la conditionnelle complète de chaque variable ne dépend que de sa couverture de Markov (Markov blanket) : ses parents, ses enfants, et les co-parents de ses enfants. Cela simplifie considérablement le calcul des conditionnelles.

Quand les conditionnelles ne sont pas de forme connue, on peut remplacer l’échantillonnage direct par un pas de Metropolis-Hastings pour les variables concernées (« Metropolis-within-Gibbs »). C’est l’approche utilisée par les logiciels comme BUGS et JAGS.

Variantes du Gibbs sampling

Systematic scan vs Random scan

Le systematic scan parcourt les variables dans un ordre fixe (1, 2, …, d, 1, 2, …) à chaque itération. C’est l’implémentation la plus courante. Le random scan choisit la variable à mettre à jour uniformément au hasard à chaque pas. Les deux convergent vers la même distribution stationnaire, mais le random scan a des propriétés théoriques légèrement plus simples.

Blocked Gibbs Sampling

Au lieu de mettre à jour une variable à la fois, le blocked Gibbs sampling regroupe plusieurs variables en blocs et met à jour chaque bloc conjointement. Par exemple, si x₁ et x₂ sont fortement corrélés, les mettre à jour ensemble (en échantillonnant P(x₁, x₂ | x₃, …, xd)) réduit l’autocorrélation de la chaîne.

Le cas extrême du blocking est un seul bloc contenant toutes les variables, ce qui revient à échantillonner directement la distribution jointe. Le compromis est entre l’efficacité de la chaîne (blocs plus grands = meilleure exploration) et la difficulté d’échantillonner la conditionnelle de chaque bloc (blocs plus grands = conditionnelles plus complexes).

Collapsed Gibbs Sampling

Le collapsed Gibbs sampling marginalise (intègre) certaines variables avant de construire la chaîne. Si on peut marginaliser z analytiquement, on échantillonne P(x | données) au lieu de P(x, z | données), ce qui réduit la dimensionnalité de la chaîne et améliore le mixing.

C’est la technique standard pour le Latent Dirichlet Allocation (LDA). En marginalisant les distributions de Dirichlet (les distributions de thèmes par document et de mots par thème), on obtient un Gibbs sampler où les seules variables à échantillonner sont les assignations de thèmes pour chaque mot, avec des conditionnelles qui prennent la forme simple d’une distribution Dirichlet-multinomiale.

Applications majeures

Topic Modeling (LDA)

Le Latent Dirichlet Allocation est l’application emblématique du Gibbs sampling en NLP. LDA modélise chaque document comme un mélange de thèmes, et chaque thème comme une distribution sur les mots. Le collapsed Gibbs sampler pour LDA est devenu le standard de l’inférence dans ce modèle.

L’algorithme est remarquablement simple : pour chaque mot dans le corpus, on échantillonne son assignation de thème conditionnellement aux assignations de tous les autres mots. La conditionnelle est proportionnelle au produit de deux termes : la fréquence du thème dans le document courant et la fréquence du mot dans le thème candidat.

Modèles de mélange

Les modèles de mélange gaussien (GMM) s’échantillonnent naturellement par Gibbs. On introduit des variables latentes d’allocation Z qui assignent chaque observation à une composante du mélange. Le Gibbs sampler alterne entre échantillonner les Z (quelle composante pour chaque point ?) conditionnellement aux paramètres, et échantillonner les paramètres (moyennes, variances, poids) conditionnellement aux Z. Avec des priors conjugués (Normal-Inverse-Wishart pour les paramètres, Dirichlet pour les poids), chaque conditionnelle a une forme analytique.

Modèles hiérarchiques bayésiens

Les modèles hiérarchiques (par exemple, des effets aléatoires par groupe avec des hyperparamètres communs) sont le terrain naturel du Gibbs sampling. La structure hiérarchique produit des conditionnelles complètes bien définies à chaque niveau. Les logiciels BUGS (Bayesian inference Using Gibbs Sampling) et JAGS (Just Another Gibbs Sampler) ont été spécifiquement conçus pour ces modèles.

Imputation de données manquantes (MICE)

La méthode MICE (Multiple Imputation by Chained Equations) utilise un Gibbs sampler pour imputer les données manquantes dans des datasets multivariés. Chaque variable manquante est échantillonnée conditionnellement aux valeurs observées et aux imputations des autres variables. La méthode gère naturellement les types de données mixtes (continues, binaires, catégorielles) en spécifiant un modèle de régression conditionnel pour chaque variable.

Traitement d’images

L’application originale de Geman et Geman (1984) était la restauration d’images. Les champs de Markov aléatoires (MRF) modélisent les relations de voisinage entre pixels. Le Gibbs sampling permet d’échantillonner la configuration de pixels qui maximise la postérieure, en mettant à jour chaque pixel conditionnellement à ses voisins. La combinaison du Gibbs sampling avec le simulated annealing (« refroidissement simulé ») dans ce contexte a donné naissance à des techniques efficaces pour la segmentation et la déconvolution d’images.

NLP et modèles de séquences

Au-delà de LDA, le Gibbs sampling est utilisé pour les modèles de Markov cachés (HMM) en traitement automatique du langage : étiquetage morphosyntaxique, reconnaissance d’entités nommées, segmentation de phrases. Les variables latentes discrètes (états cachés) se prêtent naturellement au Gibbs, car les conditionnelles sont des distributions catégorielles simples dans un modèle HMM avec priors conjugués. Les modèles bayésiens non paramétriques (processus de Dirichlet, Chinese Restaurant Process) utilisent également le Gibbs sampling pour inférer automatiquement le nombre de clusters ou de thèmes à partir des données, sans le fixer à l’avance.

Génétique et bioinformatique

L’inférence de structures de population, la localisation de QTL (Quantitative Trait Loci), et l’analyse de réseaux de régulation génique utilisent massivement le Gibbs sampling. Le logiciel STRUCTURE, l’un des plus cités en génétique des populations, repose sur un Gibbs sampler pour inférer l’ascendance de chaque individu à partir de données génomiques. La nature discrète des génotypes et l’existence de priors conjugués rendent le Gibbs sampling particulièrement adapté à ces applications.

Implémentation complète : modèle Normal-Gamma

Voici un Gibbs sampler complet pour inférer la moyenne μ et la précision λ (inverse de la variance) de données gaussiennes avec des priors conjugués :

import numpy as np
from scipy import stats

def gibbs_normal_gamma(data, n_iter=5000, burn_in=1000):
    """
    Gibbs sampling pour le modèle :
        x_i | mu, lambda ~ N(mu, 1/lambda)
        mu | lambda ~ N(mu_0, 1/(kappa_0 * lambda))
        lambda ~ Gamma(alpha_0, beta_0)
    """
    n = len(data)
    x_bar = np.mean(data)
    s2 = np.var(data)

    # Hyperparamètres du prior
    mu_0, kappa_0 = 0.0, 1.0    # Prior sur mu
    alpha_0, beta_0 = 1.0, 1.0  # Prior sur lambda

    # Initialisation
    mu = x_bar
    lam = 1.0 / s2

    samples_mu = np.zeros(n_iter)
    samples_lam = np.zeros(n_iter)

    for t in range(n_iter):
        # 1. Échantillonner mu | lambda, data
        kappa_n = kappa_0 + n * lam
        mu_n = (kappa_0 * mu_0 + n * lam * x_bar) / kappa_n
        mu = np.random.normal(mu_n, 1.0 / np.sqrt(kappa_n))

        # 2. Échantillonner lambda | mu, data
        alpha_n = alpha_0 + n / 2.0
        beta_n = beta_0 + 0.5 * np.sum((data - mu) ** 2)
        lam = np.random.gamma(alpha_n, 1.0 / beta_n)

        samples_mu[t] = mu
        samples_lam[t] = lam

    # Jeter le burn-in
    return samples_mu[burn_in:], samples_lam[burn_in:]

# Données simulées
np.random.seed(42)
data = np.random.normal(loc=5.0, scale=2.0, size=100)

mu_samples, lam_samples = gibbs_normal_gamma(data, n_iter=10000)
print(f"Moyenne postérieure de mu : {mu_samples.mean():.3f}")
print(f"IC 95% mu : [{np.percentile(mu_samples, 2.5):.3f}, {np.percentile(mu_samples, 97.5):.3f}]")
print(f"Écart-type postérieur : {(1/np.sqrt(lam_samples)).mean():.3f}")

Gibbs vs Metropolis-Hastings vs HMC

Outils logiciels

BUGS (Bayesian inference Using Gibbs Sampling) : le logiciel historique (années 1990) qui a rendu le Gibbs sampling accessible aux statisticiens. Versions WinBUGS et OpenBUGS.

JAGS (Just Another Gibbs Sampler) : réécriture moderne de BUGS, cross-plateforme, interfaçable avec R (rjags) et Python (pyjags). Utilise le Gibbs sampling quand les conditionnelles sont conjuguées, et Metropolis-Hastings ou slice sampling sinon.

Stan : successeur spirituel de BUGS/JAGS qui utilise HMC/NUTS au lieu du Gibbs sampling. Stan est généralement préféré aujourd’hui pour les modèles continus, mais ne supporte pas les variables discrètes latentes (où Gibbs/JAGS restent nécessaires).

PyMC : framework Python qui supporte Gibbs, NUTS et d’autres samplers. Sélectionne automatiquement l’algorithme le plus approprié pour chaque variable.