BSS applied on EEG signal

SISEA - Problèmes Inverses

February 14, 2018

1 Prerequisites

We aime to reconstruct the unknown source vector s given its observation :

x = As + (cid:15),

(1)

where A is the known M × N mixing matrix (M (cid:28) N ), x the observation column vector (M × 1),

and (cid:15) the white Gaussian noise of unknown variances σ2

(cid:15) .

To impose sparsity in the solution, it is further supposed that s follows a Bernoulli-Gaussian (BG)

model, an independent, identically distributed (iid) process deﬁned in two stages. Firstly, the sparse

nature is governed by the Bernoulli law:

with q = [q1, . . . , qN ]t a binary sequence and L = (cid:80)N

Secondly, amplitudes s = [s1, . . . , sN ]t are assumed iid zero-mean Gaussian conditionally to q:

n=1 qn, the number of non-zero entries of q.

P (q) = λL(1 − λ)N −L

(2)

s | q ∼ N (0, σ2

s diag(q)),

(3)

where diag(q) denotes a diagonal matrix whose diagonal is q.

(cid:15) , σ2

The problem becomes the estimation of Θ = {s, q, λ, σ2

distribution writes :

s } given x, for which the joint posterior

P (Θ | x) ∝ g(x − As; σ2

(cid:15) IM ) P (q; λ) g(s; σ2

s diag(q)) P (σ2

(cid:15) )P (σ2

s )P (λ)

(4)

where g(·; R) denotes the centered Gaussian density of covariance R. Conjugate prior laws are adopted

for λ, σ2

(cid:15) and σ2

s :

λ ∼ Be(1, 1)?, σ2

(cid:15) ∼ IG(1, 1)? σ2

(cid:15) ∼ IG(1, 1e−4)?

where IG and Be represents respectively the inverse Gamma and Beta distributions.

2 Gibbs Sampler

According to the Monte Carlo principle, a posterior mean estimator of the unknown random variables

can be approximated by:

(cid:98)Θ =

I − J

(cid:88)

k=J+1

Θ(k),

(5)

where the sum extends over the last I −J samples. In the MCMC framework, the samples are generated

iteratively, so that asymptotically converges in distribution to the joint posterior probability in eq. (4).

Metropolis-Hastings algorithm

1(cid:13) current conﬁguration Θ(k)

2(cid:13) draw Θ(cid:48) with a proposition law Q(Θ(cid:48) | Θ(k))

3(cid:13) accept Θ(k+1) = Θ(cid:48) with probability

(cid:26)

min

P (Θ(cid:48))Q(Θ(k) | Θ(cid:48))

P (Θ(k))Q(Θ(cid:48) | Θ(k))

(cid:27)

Gibbs algorithm

1(cid:13) current conﬁguration Θ(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)

4(cid:13) repeat from 1(cid:13)

Exercise 2.1 Show that the Gibbs algorithm

1. is a special case of the Metropolis-Hastings algorithm for which Θ(cid:48) is always accepted;

2. satisﬁes the equilibrium condition and thus Θ(k) converges to P (Θ | z) in distribution.

Exercise 2.2 write the following conditional probabilities from the joint probability in Eq (4):

1. qn | Θ \ {qn, sn}, x,

2. sn | Θ \ {sn}, x,

3. σ2

(cid:15) | Θ \ σ2

(cid:15) , x,

4. σ2

s | Θ \ σ2

s , x,

5. λ | Θ \ λ, x,

N.B. the Beta distribution writes Beta(x; α, β) =

IG(x; α, β) = βα

Γ(α) x−(α+1) exp(−β/x).

Conditional laws:

We ﬁrst detail the joint probability in Eq (4) by :

B(α,β) xα−1(1 − x)β−1, and the inverse Gamma

P (Θ | x) ∝ (cid:0)2πσ2

(cid:15)

(cid:1)−M/2 exp

(cid:26)

−

2σ2

(cid:15)

(cid:27)

(cid:107)x − As(cid:107)2

(cid:80) qn(1 − λ)N −(cid:80) qn

· λ

· (cid:0)(2π)N det(σ2

s diag(q))(cid:1)−1/2

(cid:26)

exp

−

st(σ2

s diag(q))−1s

(cid:27)

· P (σ2

s ) · P (σ2

(cid:15) ).

s )N (cid:81) qn due to its diagonal structure, while a zero-variance normal

Notice that det(σ2

distribution can be considered a Dirac distribution. Notice also that for each binary variable qn, we

should have the normalization constraint

s diag(q)) = (σ2

P (qn = 0 | Θ \ {qn, sn}, x) + P (qn = 1 | Θ \ {qn, sn}, x) = 1.

(6)

Its conditional probability can thus be derived by calculating respectively the two terms up to a

normalization factor :

• for P (qn = 0 | Θ \ {qn, sn}, x), according to the B-G model, sn = 0, thus

P (qn = 0, sn | Θ \ {qn, sn}, x) ∝ exp

−

(cid:26)

(cid:107)x − A−ns−n − Ansn(cid:107)2

(cid:27)

· (1 − λ) · δ0(sn);

2σ2

(cid:15)

for which the A−n and s−n represents the matrix A without the n-th column and the vector s

without its n-th element respectively; we can easily integrate out sn and have the marginal form:

(cid:26)

(cid:27)

(cid:107)x − A−ns−n(cid:107)2

· (1 − λ);

(7)

P (qn = 0 | Θ \ {qn, sn}, x) ∝ exp

−

2σ2

(cid:15)

• likewise for P (qn = 1, sn | Θ \ {qn, sn}, x):

P (qn = 1, sn | Θ \ {qn, sn}, x) ∝ exp

−

(cid:26)

2σ2

(cid:15)

(cid:107)x − A−ns−n − Ansn(cid:107)2

(cid:27)

· λ ·

(cid:112)2πσ2

exp(−

2σ2

)

the quadratic form inside the exponential terms suggests that sn follows a normal distribution

(N (µn, σ2

n)) and can be integrated out analytically to yield :

P (qn = 1 | Θ \ {qn, sn}, x) ∝ exp

−

(cid:26)

2σ2

(cid:15)

(cid:107)x − A−ns−n(cid:107)2

(cid:27)

λσn

σs

· exp

(cid:27)

(cid:26) µ2

2σ2

(8)

where σ2

n =

s (cid:107)An(cid:107)2 and µn = σ2

(cid:15) σ2

σ2

(cid:15) +σ2

σ2

(cid:15)

n(x − A−ns−n).

exp

(cid:111)

(cid:110) µ2

2σ2

By combining eqs. (6)(7)(8), we conclude that qn | Θ \ {qn, sn}, x ∼ Bi(

λσn

σs

The conditional probability of sn | Θ \ {sn}, x is straight forward :

sn ∼ N (µn, σ2

obvious thanks to the conjugate prior laws.

if qn = 0, then sn = 0, else

n) as derived previously. The conditional probabilities of the rest three parameters are

νn

νn+1−λ) for which νn =

Exercise 2.3 From the Gibbs sampling scheme, derive the pseudocode in Tab. 1

Table 1: Pseudocode for the Gibbs sampler

Sample (qi, si) ––––––––––-

s /(σ2

(cid:15) )At

(cid:15) + σ2

iei

s , λ using prior laws

% ––––––––––- Step 1:

e ← x − As

for i = 1 to N do

1: % Initialization

2: q ← 0; s ← 0

(cid:15) , σ2

3: Sample σ2

4: repeat

10:

11:

12:

13:

14:

15:

end for

16:

% ––––––––––- Step 2:

17:

Sample σ2

18:

Sample λ ––––––––––––

% ––––––––––- Step 3:

19:

Sample λ ∼ Be(1 + L, 1 + N − L), with L = (cid:80)

20:

s ––––––––––-

% ––––––––––- Step 4:

21:

Sample σ2

22:

23: until Convergence

ei ← e + Aisi % Ai is the i-th column of A

s (cid:107)Ai(cid:107)2)

(cid:15) σ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

Sample σ2

s ∼ IG(L/2+?, (cid:107)s(cid:107)2 /2+?)

(cid:15) ∼ IG(M/2+?, (cid:107)x − As(cid:107)2 /2+?)

(cid:15) ––––––––––-

Sample σ2

i /(2σ2

i ))

i ) if qi = 1, si = 0 otherwise

n qn

3 Matlab Experiment :

1. launch the TP_inverse_problem.m ﬁle in matlab and visualize the source distribution;

2. understand the direct problem process (code given) of linear source mixture and noise addition;

3. use the Gibbs sampler to perform the "deconvolution" and compare the results.

BSS Applied on EEG Signal

BSS Applied on EEG Signal