Probl`eme inverse d’EEG: TP1 - ´echantillonneur de Gibbs
ESIR 2 - Probl`emes Inverses
Universit´e de Rennes1
1
Introduction
L’´echantillonneur de Gibbs fait partie de la famille des approches Bay´esiennes. Les m´ethodes ap-
partenant a cette famille sont bas´ees sur un modele probabiliste des donn´ees. Plus particuli`erement,
le vecteur de mesures du potentiel ´electrique x, les signaux des dipˆoles de l’espace sources s, et le
bruit n `a un instant donn´e t sont consid´er´es comme r´ealisations de processus stochastiques vectoriels.
L’objectif des approches Bay´esiennes consiste `a estimer la probabilit´e a posteriori P (s|x) des sources
qui, suivant le th´eor`eme de Bayes, est donn´e par
P (s|x) =
P (x|s)P (s)
P (x)
.
(1)
Pour reconstruire les sources, on cherche `a estimer la moyenne de la probabilit´e a posteriori P (s|x)
pour un vecteur d’observations x donn´e.
L’a priori des observations P (x) est consid´er´e comme une constante de normalisation. En supposant
que le bruit suit un processus stochastique gaussien centr´e avec variance σ2
n inconnue, la vraisemblance
P (x|s) suit ´egalement un mod`ele gaussien avec moyenne As et variance σ2
n:
P (x|s) = g(x − As; σ2
nIN ).
(2)
Dans les approches Bay´esiennes, les hypoth`eses (a priori) sur les sources sont incorpor´ees dans le
choix de la loi de probabilit´e des sources P (s). Dans ce TP, pour obtenir un vecteur de sources
parcimonieux, nous supposons que le vecteur s suit un mod`ele Bernoulli-Gaussien. Plus pr´ecis´ement,
il s’agit d’un processus vectoriel d´efini en deux ´etapes dont les ´el´ements sont distribu´ees de fa¸con
identique et ind´ependant. Premi`erement, la parcimonie des sources est caract´eris´ee par la loi Bernoulli
P (q) = λL(1 − λ)D−L
(3)
o`u q = [q1, . . . , qD]t est une s´equence binaire qui d´ecrit l’activation des dipˆoles et L = (cid:80)D
d=1 qd est le
nombre d’´el´ements non-z´ero correspondant au nombre de dipˆoles actifs. Deuxi`emement, les amplitudes
Publicité
s = [s1, . . . , sD]t des sources sont suppose´ees suivre une loi gaussienne centr´ee de fa¸con identique et
ind´ependant conditionnellement a q (c’est-a-dire uniquement les amplitudes des dipˆoles actifs pour
lesquels qd = 1 suivent une loi gaussienne, les amplitudes des dipˆoles non-actifs ´etant z´eros):
o`u 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
L’objectif consiste alors a estimer tous les parametres inconnus, r´esum´e dans le vecteur
Θ = {s, q, λ, σ2
n, σ2
s },
`a partir des observations x. La loi conjointe a posteriori du vecteur Θ est donn´e par
P (Θ | x) ∝ g(x − As; σ2
nIN ) P (q; λ) g(s; σ2
ou g(·; R) denote la loi gaussienne centr´ee 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´esente la loi bˆeta 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
´Echantillonneur de Gibbs
Pour estimer la moyenne de la probabilit´e a posteriori du vecteur de param`etres Θ (et ainsi des vecteurs
q et s caract´erisant les dipˆoles actifs et leurs amplitudes), nous utilisons l’´echantillonneur de Gibbs
pour tirer des r´ealisations al´eatoires des variables inconnues. Selon le principe de Monte Carlo, la
moyenne a posteriori peut ensuite ˆetre approxim´ee par:
(cid:98)Θ =
1
I − J
Publicité
I
(cid:88)
k=J+1
Θ(k),
(6)
ou la somme s’´etend sur les derniers I − J ´echantillons. Dans le cadre des modeles MCMC, les
´echantillons sont g´en´er´es it´erativement tel que leur distribution converge asymptotiquement vers la loi
conjointe a posteriori de l’´equation (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´e par l’´echantillonneur
de Gibbs, mais seulement le vecteur q qui caract´erise les dipˆoles actifs. Les amplitudes des dipˆoles
actifs sont ensuite estim´es par un estimateur MAP (maximum a posteriori).
2 Pr´eparation
• ´Ecrire les lois de probabilit´e conditionnelle `a partir de la loi conjointe de l’´equation Eq (5):
1. qd | Θ \ {qd, sd}, x,
2. sd | Θ \ {sd}, x,
3. λ | Θ \ λ, x,
N.B. la loi Beta s’´ecrit Beta(x; α, β) = 1
B(α,β) xα−1(1 − x)β−1
• A partir du sch´ema de l’´echantillonneur de Gibbs, d´eduire le pseudocode dans Tableau 1.
2
n, σ2
s , α, β
Sample (qi, si) ---------------------
Table 1: Pseudocode de l’´echantillonneur de Gibbs
% --------------------- Step 1:
e ← x − As
for i = 1 to D do
Publicité
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
λ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
Publicité
iei
s /(σ2
n)At
i /(2σ2
d qd
i ))
Sample λ ------------------------
3 Manipulations
3.1 Choix des param`etres α, β, σ2
n, σ2
s
Nous allons d’abord d´eterminer des valeurs appropri´ees pour les param`etres σ2
s . Pour cela, nous
nous servons de la matrice de bruit et de la matrice des sources simul´ees (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’´echantillon correspondant
au maximum de la pointe ´epileptique pour SNR=10.
Le parametre λ est tir´e al´eatoirement d’une loi bˆeta 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`eme de localisation de sources trait´e dans ce TP.
n et σ2
3.2 Impl´ementation du Gibbs sampler
• Compl´eter la fonction Matlab Gibbs sampler `a l’aide du pseudocode dans Tableau 1 (attention
aux parenth`eses!) en respectant les consignes dans la fonction sur les noms des variables
utilis´ees. Pour ´echantillonner les lois bˆeta et gaussienne, utiliser les fonctions betarnd et
randn. Pour obtenir un ´echantillon 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´ees correspondant au maximum de la pointe. Visualiser le vecteur ˆs estim´e `a l’aide de
la fonction trisurf et comparer les sources estim´ees avec les sources originales.
• Analyser l’influence du RSB (allant de 0.1 `a 10) sur le r´esultat du Gibbs sampler.
3