Problème inverse d’EEG: TP1 - échantillonneur de Gibbs
ESIR 2 - Problèmes Inverses
Université de Rennes1
1
Introduction
L’échantillonneur de Gibbs fait partie de la famille des approches Bayésiennes. Les méthodes ap-
partenant a cette famille sont basées sur un modele probabiliste des données. Plus particulièrement,
le vecteur de mesures du potentiel électrique x, les signaux des dipôles de l’espace sources s, et le
bruit n à un instant donné t sont considérés comme réalisations de processus stochastiques vectoriels.
L’objectif des approches Bayésiennes consiste à estimer la probabilité a posteriori P (s|x) des sources
qui, suivant le théorème de Bayes, est donné par
P (s|x) =
P (x|s)P (s)
P (x)
.
(1)
Pour reconstruire les sources, on cherche à estimer la moyenne de la probabilité a posteriori P (s|x)
pour un vecteur d’observations x donné.
L’a priori des observations P (x) est considéré comme une constante de normalisation. En supposant
que le bruit suit un processus stochastique gaussien centré avec variance σ2
n inconnue, la vraisemblance
P (x|s) suit également un modèle gaussien avec moyenne As et variance σ2
n:
P (x|s) = g(x − As; σ2
nIN ).
(2)
Dans les approches Bayésiennes, les hypothèses (a priori) sur les sources sont incorporées dans le
choix de la loi de probabilité des sources P (s). Dans ce TP, pour obtenir un vecteur de sources
parcimonieux, nous supposons que le vecteur s suit un modèle Bernoulli-Gaussien. Plus précisément,
il s’agit d’un processus vectoriel défini en deux étapes dont les éléments sont distribuées de façon
identique et indépendant. Premièrement, la parcimonie des sources est caractérisée par la loi Bernoulli
P (q) = λL(1 − λ)D−L
(3)
où q = [q1, . . . , qD]t est une séquence binaire qui décrit l’activation des dipôles et L = (cid:80)D
d=1 qd est le
nombre d’éléments non-zéro correspondant au nombre de dipôles actifs. Deuxièmement, les amplitudes
s = [s1, . . . , sD]t des sources sont supposeées suivre une loi gaussienne centrée de façon identique et
indépendant conditionnellement a q (c’est-a-dire uniquement les amplitudes des dipôles actifs pour
lesquels qd = 1 suivent une loi gaussienne, les amplitudes des dipôles non-actifs étant zéros):
où diag(q) est une matrice diagonale dont la diagonal est q et σ2
s est la variance inconnue des sources.
s | q ∼ N (0, σ2
s diag(q)),
(4)
1
Publicité
L’objectif consiste alors a estimer tous les parametres inconnus, résumé dans le vecteur
Θ = {s, q, λ, σ2
n, σ2
s },
à partir des observations x. La loi conjointe a posteriori du vecteur Θ est donné par
P (Θ | x) ∝ g(x − As; σ2
nIN ) P (q; λ) g(s; σ2
ou g(·; R) denote la loi gaussienne centrée de covariance R. Pour le parametre λ, on utilise l’a priori
suivante:
s diag(q)) P (σ2
n)P (σ2
s )P (λ)
(5)
ou Be représente la loi bêta et ou α et β sont des hyperparametres. Pour simplifier le probleme, nous
n et σ2
utilisons dans ce TP des valeurs fixes pour σ2
s .
λ ∼ Be(α, β)
1.1
Échantillonneur de Gibbs
Pour estimer la moyenne de la probabilité a posteriori du vecteur de paramètres Θ (et ainsi des vecteurs
q et s caractérisant les dipôles actifs et leurs amplitudes), nous utilisons l’échantillonneur de Gibbs
pour tirer des réalisations aléatoires des variables inconnues. Selon le principe de Monte Carlo, la
moyenne a posteriori peut ensuite être approximée par:
(cid:98)Θ =
1
I − J
I
(cid:88)
k=J+1
Θ(k),
(6)
ou la somme s’étend sur les derniers I − J échantillons. Dans le cadre des modeles MCMC, les
échantillons sont générés itérativement tel que leur distribution converge asymptotiquement vers la loi
conjointe a posteriori de l’équation (5).
Gibbs algorithm
1(cid:13) current configuration Θ(k)
2(cid:13) choose (cyclically or randomly) i
3(cid:13) draw Θ(k+1)
, z
4(cid:13) repeat from 1(cid:13)
∼ Θi | Θ(k) \ Θ(k)
i
i
Pour une meilleure estimation des amplitudes s, nous n’utilisons pas le vecteur s estimé par l’échantillonneur
Publicité
de Gibbs, mais seulement le vecteur q qui caractérise les dipôles actifs. Les amplitudes des dipôles
actifs sont ensuite estimés par un estimateur MAP (maximum a posteriori).
2 Préparation
• Écrire les lois de probabilité conditionnelle à partir de la loi conjointe de l’équation Eq (5):
1. qd | Θ \ {qd, sd}, x,
2. sd | Θ \ {sd}, x,
3. λ | Θ \ λ, x,
N.B. la loi Beta s’écrit Beta(x; α, β) = 1
B(α,β) xα−1(1 − x)β−1
• A partir du schéma de l’échantillonneur de Gibbs, déduire le pseudocode dans Tableau 1.
2
n, σ2
s , α, β
Sample (qi, si) ---------------------
Table 1: Pseudocode de l’échantillonneur de Gibbs
% --------------------- Step 1:
e ← x − As
for i = 1 to D do
1: % Initialization
2: q ← 0; s ← 0
3: Sample λ using the prior law and choose appropriate values for σ2
4: repeat
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19: until Convergence
ei ← e + Aisi % Ai is the i-th column of A
s (cid:107)Ai(cid:107)2)
nσ2
i ← σ2
σ2
i /σ2
µi ← (σ2
νi ← λ(σi/σs) exp (µ2
Publicité
λi ← νi/(νi + 1 − λ)
Sample qi ∼ Bi(λi)
Sample si ∼ N (µi, σ2
e ← ei − Aisi % Update e
end for
% --------------------- Step 2:
Sample λ ∼ Be(α + L, β + D − L), with L = (cid:80)
i ) if qi = 1, si = 0 otherwise
n + σ2
iei
s /(σ2
n)At
i /(2σ2
d qd
i ))
Sample λ ------------------------
3 Manipulations
3.1 Choix des paramètres α, β, σ2
n, σ2
s
Nous allons d’abord déterminer des valeurs appropriées pour les paramètres σ2
s . Pour cela, nous
nous servons de la matrice de bruit et de la matrice des sources simulées (qui sont connues dans notre
cas; en pratique, on ne disposerait pas de cette information et devrait estimer ces valeurs d’une autre
maniere). Pour cela, calculer les variances du bruit et des sources actives a l’échantillon correspondant
au maximum de la pointe épileptique pour SNR=10.
Le parametre λ est tiré aléatoirement d’une loi bêta dont il faut ajuster les parametres α et β.
Visualiser la loi Beta(x; α, β) (utiliser la fonction betadist) pour (α, β) = (0.5, 0.5), (1, 10), (2, 100),
(10, 1), et (100, 2). Par la suite, nous allons utiliser les valeurs (α, β) = (1.01, 2000). Justifier le choix
de ces valeurs pour le problème de localisation de sources traité dans ce TP.
n et σ2
3.2 Implémentation du Gibbs sampler
• Compléter la fonction Matlab Gibbs sampler à l’aide du pseudocode dans Tableau 1 (attention
aux parenthèses!) en respectant les consignes dans la fonction sur les noms des variables
utilisées. Pour échantillonner les lois bêta et gaussienne, utiliser les fonctions betarnd et
randn. Pour obtenir un échantillon e de la loi Bernoulli, vous pouvez vous servir de la loi
uniforme: e=rand(1)<lambda.
• Dans le scripte TP inverse problems.m, appliquer la fonction Gibbs sampler au vecteur
de données correspondant au maximum de la pointe. Visualiser le vecteur ˆs estimé à l’aide de
la fonction trisurf et comparer les sources estimées avec les sources originales.
• Analyser l’influence du RSB (allant de 0.1 à 10) sur le résultat du Gibbs sampler.
3