Problème inverse d’EEG: TP1 - Échantillonneur de Gibbs

Ce document présente un travail pratique (TP) consacré à l’échantillonneur de Gibbs appliqué au problème inverse en électroencéphalographie (EEG). Il s’adresse aux étudiants de niveau master en traitement du signal ou en probabilités appliquées, souhaitant comprendre et implémenter une méthode bayésienne pour la localisation de sources cérébrales à partir de mesures EEG.

Introduction au problème inverse d’EEG et approche bayésienne

Le problème inverse d’EEG consiste à estimer les sources électriques cérébrales à partir des mesures des potentiels électriques enregistrés à la surface du crâne. L’approche bayésienne modélise les données et les sources comme des processus stochastiques vectoriels. Le vecteur de mesures x, les signaux des dipôles s, et le bruit n à un instant t sont considérés comme des réalisations aléatoires.

La probabilité a posteriori des sources, notée P(s|x), est donnée par le théorème de Bayes :

P(s|x) = P(x|s) P(s) / P(x)

où P(x) est une constante de normalisation. Le bruit est supposé gaussien centré avec une variance inconnue σn², ce qui implique que la vraisemblance suit un modèle gaussien :

P(x|s) = g(x − As; σn² I_N)

où A est la matrice de transfert (modèle direct), et g(·; R) désigne une loi gaussienne centrée de covariance R.

Pour favoriser une solution parcimonieuse, le vecteur des sources s est modélisé par un processus Bernoulli-Gaussien. Ce modèle est défini en deux étapes :

• Une variable binaire q = [q1, ..., qD]^t décrit l’activation des dipôles, avec une loi de Bernoulli :

P(q) = λ^L (1 − λ)^(D−L)

où L = ∑_d q_d est le nombre de dipôles actifs.

• Les amplitudes des sources conditionnelles à q suivent une loi gaussienne centrée :

s | q ∼ N(0, σs² diag(q))

Les amplitudes des dipôles non actifs (qd = 0) sont nulles.

L’objectif est d’estimer le vecteur de paramètres inconnus Θ = {s, q, λ, σn², σs²} à partir des observations x. La loi conjointe a posteriori s’écrit :

P(Θ | x) ∝ g(x − As; σn² I_N) P(q; λ) g(s; σs² diag(q)) P(σn²) P(σs²) P(λ)

où λ suit une loi bêta Be(α, β) avec des hyperparamètres α et β. Dans ce TP, σs² est fixé pour simplifier.

Échantillonneur de Gibbs

L’échantillonneur de Gibbs est une méthode de Monte Carlo par chaînes de Markov (MCMC) permettant de générer des échantillons de la distribution conjointe a posteriori de Θ. La moyenne a posteriori est approximée par la moyenne des échantillons générés après une phase de convergence :

Θ̂ = (1 / (I − J)) ∑_{k=J+1}^I Θ^(k)

où J est le nombre d’itérations de « burn-in » et I le nombre total d’itérations.

Le principe de l’algorithme de Gibbs est :

1. Initialiser la configuration Θ^(0).
2. Choisir cycliquement ou aléatoirement une variable Θ_i à mettre à jour.
3. Tirer un nouvel échantillon Θ_i^(k+1) selon la loi conditionnelle P(Θ_i | Θ \ Θ_i).
4. Répéter.

Pour une meilleure estimation des amplitudes s, on utilise uniquement le vecteur q estimé par Gibbs pour identifier les dipôles actifs, puis on calcule les amplitudes via un estimateur MAP (maximum a posteriori).

Pseudocode de l’échantillonneur de Gibbs

% Initialisation
q ← 0; s ← 0
Sample λ selon la loi a priori (bêta)
Fixer σn² et σs²

repeat
  e ← x − As
  for i = 1 to D do
    ei ← e + A_i s_i   % A_i : i-ème colonne de A
    σi² ← (σn² σs²) / (σs² ||A_i||² + σn²)
    µi ← (σi² / σn²) A_i^T e_i
    νi ← λ (σi / σs) exp(µi² / (2 σi²))
    λi ← νi / (νi + 1 − λ)
    Sample q_i ∼ Bernoulli(λi)
    Sample s_i ∼ N(µi, σi²) si q_i = 1, sinon s_i = 0
    e ← e − A_i s_i
  end for

  Sample λ ∼ Beta(α + L, β + D − L), avec L = ∑ q_i
until convergence

Préparation et choix des paramètres

Avant l’implémentation, il est nécessaire de déterminer des valeurs adaptées pour les paramètres σn², σs², α et β :

• Calculer les variances du bruit et des sources actives à partir des données simulées au maximum de la pointe épileptique pour un rapport signal sur bruit (RSB) de 10.
• Visualiser la loi bêta Beta(x; α, β) pour différents couples (α, β) : (0.5, 0.5), (1, 10), (2, 100), (10, 1), (100, 2) afin de comprendre leur influence.
• Choisir (α, β) = (1.01, 2000) pour favoriser une forte parcimonie, c’est-à-dire un faible taux d’activation des dipôles, adapté au problème de localisation des sources EEG.

Implémentation du Gibbs sampler

• Compléter la fonction Matlab de l’échantillonneur de Gibbs en suivant le pseudocode ci-dessus. Attention à la syntaxe et aux noms des variables.
• Utiliser les fonctions Matlab betarnd pour échantillonner la loi bêta, randn pour la loi gaussienne, et pour la loi de Bernoulli, générer un échantillon via rand(1) < lambda.
• Appliquer l’échantillonneur au vecteur de données correspondant au maximum de la pointe épileptique.
• Visualiser les sources estimées ˆs avec la fonction trisurf et comparer avec les sources originales simulées.
• Étudier l’influence du rapport signal sur bruit (RSB) variant de 0.1 à 10 sur la qualité de la reconstruction.

Glossaire des termes clés

• Bernoulli-Gaussien : Modèle probabiliste combinant une variable binaire (Bernoulli) pour l’activation et une variable gaussienne pour l’amplitude conditionnelle.
• Échantillonneur de Gibbs : Algorithme MCMC qui génère des échantillons d’une distribution conjointe en échantillonnant successivement chaque variable conditionnellement aux autres.
• Estimateur MAP (Maximum a posteriori) : Estimation qui maximise la probabilité a posteriori d’un paramètre.
• Hyperparamètres : Paramètres des lois a priori, ici α et β pour la loi bêta de λ.
• Loi bêta (Be) : Distribution continue sur [0,1] utilisée comme a priori pour des probabilités, paramétrée par α et β.
• Modèle direct (A) : Matrice de transfert reliant les sources s aux mesures x.
• Parcimonie : Propriété d’un vecteur où la majorité des éléments sont nuls ou proches de zéro.
• Probabilité a posteriori : Probabilité d’un paramètre conditionnellement aux données observées.
• Processus stochastique : Modèle probabiliste décrivant une évolution aléatoire dans le temps ou l’espace.
• Rapport signal sur bruit (RSB) : Rapport entre la puissance du signal et celle du bruit, indicateur de qualité des données.
• Vraisemblance : Fonction de la probabilité des données en fonction des paramètres du modèle.

Points clés à retenir

• L’approche bayésienne permet d’intégrer des connaissances a priori pour résoudre le problème inverse d’EEG.
• Le modèle Bernoulli-Gaussien favorise une solution parcimonieuse, adaptée à la localisation des sources actives.
• L’échantillonneur de Gibbs est une méthode efficace pour approximer la distribution a posteriori complexe.
• Le choix des hyperparamètres α et β influence fortement la parcimonie et la qualité de l’estimation.
• La variance du bruit et des sources doit être estimée ou fixée pour un bon fonctionnement de l’algorithme.
• La phase de burn-in est nécessaire pour garantir la convergence de la chaîne de Markov.
• Les amplitudes des sources actives sont estimées via un estimateur MAP après identification des dipôles actifs.
• La visualisation des résultats permet de comparer les sources estimées aux sources réelles simulées.