METODI MULTIVARIATI PER LA STIMA DELLA …salinari/Didattica/AutoregressiviConnettivita... ·...

49

L. Astolfi - Modelli autoregressivi per la connettività tra segnali biologici - Corso di Modelli di Sistemi Biologici

METODI MULTIVARIATI PER LA STIMA DELLA CONNETTIVITA’

Corso di Modelli di Sistemi Biologici – Prof. S. Salinari

Università di Roma “La Sapienza” - Facoltà di Ingegneria

L. Astolfi - Modelli autoregressivi per la connettività tra segnali biologici - Corso di Modelli di Sistemi Biologici 50

Cos’è un metodo multivariato?

MVARMVARx2[n]

x1[n]

x N[n]

MX1X1

X2X2

X3X3

X4X4

X5X5

METODI MULTIVARIATI: Il pattern di connettività si ottiene basandosi su un unico modello stimato sull’intero set di segnali, che tiene conto di tutte le loro relazioni reciproche

METODI BIVARIATI: un modello per ogni possibile coppia di segnali all’interno dell’insieme considerato


Modelli Multivariati AR

• Dato un set di N segnali:

• Si definisce il modello MULTIVARIATO

AUTOREGRESSIVO (MVAR) di ordine p:

[ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ]∑ ∑∑

∑ ∑∑

∑ ∑∑

= ==

= ==

= ==

+−−−−−−−=

+−−−−−−−=

+−−−−−−−=

p

1kN

p

1kNNN22N

p

1k11NN

p

1k2

p

1kNN2222

p

1k1212

p

1k1

p

1kNN1212

p

1k1111

neknxkaknxkaknxkanx

neknxkaknxkaknxkanx

neknxkaknxkaknxkanx

L

M

L

L

[ ] [ ] [ ][ ]TN21 1x1x1xX L=



• In cui:

• aij[k] = parametro autoregressivo relativo alla

coppia di segnali ji e al lag k

• ei[k] = residuo relativo al segnale i

• p = ordine del modello

• Hp:

• processo stazionario in senso lato.

• residui bianchi, a media

nulla, scorrelati tra loro e con il segnale.

[ ]nX[ ] [ ] [ ][ ]1e1e1eE N21 L=



• I parametri del modello sono N·N·p:

• A cui si aggiungono le N varianze dei residui:

[ ][ ] [ ]

[ ] [ ][ ]

[ ] [ ]

[ ] [ ][ ]

[ ] [ ]

[ ] [ ]⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

=⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

=⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

=

papa

papa

pa

2a2a

2a2a

2a

1a1a

1a1a

1a

NN1N

N111

NN1N

N111

NN1N

N111

L

MOM

L

L

L

MOM

L

L

MOM

L

⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢

⎣

⎡

=

N

2

1

ES

σ

σ

σ

M

Totale dei parametri da identificare:

N·N·p+N=N(N·p+1)



• Il modello MVAR di ordine p descrive l’evoluzione

del generico segnale xj, appartenente al set

considerato X, in funzione dei campioni agli p

istanti passati di tutti gli altri segnali x

• Secondo il concetto di causalità di Granger, ciascun

coefficiente aij[k] può essere nullo oppure diverso

da zero, a seconda che l’effetto di j su i sia causale

o meno

• => un modello MVAR costituisce ancora uno

stimatore della causalità di Granger



X1[n-3]X1[n-2]

X1[n-1]X[n]

X4[n-3]X4[n-2]X4[n-1]

X[n]

[ ]nX̂ 4X2[n-3]

X2[n-2]X2[n-1]

X3[n-3]X3[n-2]

X3[n-1]

a41 [2]

a41 [1]

a41 [3]

a42 [1]a42 [2]

a42 [3]

a43 [1]

a43 [2]

a43 [3]


MVAR nel dominio della frequenza

• Portando le x al primo membro:

Con:

a11[0]=1, a22[0]=, … , aNN[0]=1 =>

[ ] [ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ] [ ]neknxkaknxka

neknxkaknxka

neknxkaknxka

N

p

0kNNN

p

0k11N

2

p

1kNN2

p

1k221

1

p

1kNN1

p

0k111

=−++−

=−++−

=−++−

∑∑

∑∑

∑∑

==

==

==

L

M

L

L

[ ]⎥⎥

⎦

⎤

⎢⎢

⎣

⎡=

10

010a

LMOM

M



• Passando al dominio della frequenza otteniamo:

• Dove:

( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( ) ( ) ( ) ( )fEfXfAfXfAfXfA

fEfXfAfXfAfXfA

fEfXfAfXfAfXfA

NNNN22N11N

2NN2222121

1NN1212111

=+++

=+++=+++

L

M

L

L

( ) [ ]∑=

−=p

0k

fTk2jijij ekafA π

è la trasformata di Fourier di [ ] [ ] [ ] [ ][ ]pa,,1a,0aka ijijijij L=



• In forma matriciale si ha ancora:

( ) ( ) ( )fEfXfA =

( )( ) ( )

( ) ( )⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

fAfA

fAfA

fA

NN1N

N111

L

MOM

L

( )( )

( )⎥⎥

⎦

⎤

⎢⎢

⎣

⎡=

fX

fXfX

N

1

M ( )( )

( )⎥⎥

⎦

⎤

⎢⎢

⎣

⎡=

fE

fEfE

N

1

M

Con:


Predittore MVAR

L’elemento ij della matrice A è la funzione di trasferimento tra l’i-esimo ingresso e la j-esima uscita del predittore lineare MVAR.

In generale

( ) ( )fAfA jiij ≠



• L’equazione precedente può essere riscritta come segue:

dove H(f) è detta MATRICE DI TRASFERIMENTO del filtro MVAR:

( ) ( ) ( ) ( ) ( )fEfHfEfAfX 1 == −

( ) ( )( ) ( )

( ) ( )⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡== −

fHfH

fHfH

fAfH

NN1N

N111

1

L

MOM

L


Filtro generatore MVAR

L’elemento ij della matrice H è la funzione di trasferimento tra l’i-esimo ingresso e la j-esima uscita del filtro generatore MVAR.

In generale

( ) ( )fHfH jiij ≠



• Come già dimostrato per i modelli bivariati, per le densitàspettrali di potenza si ottiene:

( ) ( ) ( ) ( ) ( ) ( ) ( )fHfEfEfHfXfXfS HHH ==

• Dove è la matrice spettrale di e1[n], e2[n]

• Ma, per Hp, e1[n], e2[n], … eN[n] sono bianchi, quindi:

1) la loro densità spettrale di potenza è costante e pari alla

loro potenza σ112 , … , σNN

2

2) La loro densità di mutua potenza spettrale è costante e

pari a σij2= σji

2

( ) ( ) ( )fSfEfE EH =



• La matrice spettrale di E è indipendente dalla frequenza e pari a:

( )⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

==

2NN

21N

2ij

222

2N1

211

fE

σσσ

σσσ

Σ

L

O

MM

L

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=2NN

2ii

211

000

000

σσ

σ

Σ

L

M

MOL

Inoltre, se il modello fornisce una buona stima del processo X, Eij e Eii

sono scorrelati => la loro mutua correlazione è zero =>matrice spettrale DIAGONALE



• Quindi la matrice spettrale S espressa in funzione dei parametri del modello AR èdata da:

( ) ( ) ( ) ( ) ( ) ( ) ( )fHfHfHfEfEfHfS HHH Σ==

Con:

( ) ( ) ( )fHfHfS HΣ=

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=2NN

2ii

211

000

000

σσ

σ

Σ

L

M

MOL



e2[n]

e1[n]

e N[n]

M

x2[n]

x1[n]

x N[n]

M

( ) ( )

( ) ( )⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

fHfH

fHfH

NN1N

N111

L

MOM

L

x2[n]

x1[n]

x N[n]

M

e2[n]

e1[n]

e N[n]

M

( ) ( )

( ) ( )⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

fAfA

fAfA

NN1N

N111

L

MOM

L

Filtro generatore MVAR

Predittore lineare MVAR


Directed Transfer Function (DTF)

• Sulla base della matrice di trasferimento H(f) si definisce la DIRECTED TRANSFER FUNCTION (DTF) da j verso i:

• Più spesso usata nella sua forma normalizzata:

( ) ( ) 2ijij fHf =θ

( )( )

( )∑=

= N

1m

2

im

2

ij

ij

fH

fHfθ ( )∑

=

=N

1nin 1fθCon:



• Normalizzazione:

La somma di tutti i valori di DTF diretti verso uno stesso segnale è 1.

• Quindi il valore di DTFij normalizzato appartiene all’intervallo [0, 1] ed indica la percentuale di potenza del segnale i dovuta a j:

xi ∑=1

%303.0DTFij =>=



• Poiché , si ha che:

Il valore di θij ad una certa frequenza f0 indica l’esistenza a quella frequenza di un legame di causalità diretto dal segnale j verso il segnale i

( ) ( )ff jiij θθ ≠( ) ( )fHfH jiij ≠



21

3

Set di segnaliSet di segnali

MVARMVAR

( ) ( ) ( )fEfHfX =

[ ][ ]

[ ]⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=

nx

nx

nx

X

N

2

1

M

DTFij (f)

2

2

1

ij

N

m

H

Him(f)=

=∑

In ascissa: frequenza (in Hz)

In ordinata: valori di DTF normalizzata [0,1]

Curve in blu: funzioni di DTF dirette da j (colonne) ad i (righe)

Curve in verde (sulla diagonale) : PSD dei segnali


Partial Directed Coherence (PDC)

• Sulla base della matrice di trasferimento A(f) si definisce la PARTIAL DIRECTED COHERENCE (PDC) da j verso i:

• Anche della PDC esiste una forma normalizzata:

( ) ( ) 2ijij fAf =π

( )( )

( )∑=

= N

1m

2

im

2

ij

ij

fA

fAfπ ( )∑

=

=N

1nin 1fπCon:


Partial Directed Coherence (PDC)

• Anche per la PDC, poichési ha che:

Il valore di πij ad una certa frequenza f0 indica l’esistenza a quella frequenza di un legame di causalità diretto dal segnale j verso il segnale i

( ) ( )ff jiij ππ ≠( ) ( )fAfA jiij ≠


Confronto DTF - PDC

• Consideriamo ad esempio il seguente modello MVAR a tre canali:

[ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ][ ] [ ] [ ] [ ] [ ] [ ][ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ]ne3nx3a2nx2a1nx1anx

ne3nx3a1nx1anx

ne3nx3a2nx2a1nx1anx

32322322323

21211212

12122122121

+−⋅+−⋅+−⋅=+−⋅+−⋅=

+−⋅+−⋅+−⋅=

• Caratterizzato dai seguenti parametri:

[ ] [ ][ ]

[ ][ ]

[ ][ ]

[ ][ ]

[ ][ ]

[ ] ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

03a0

003a

03a0

3a

02a0

000

02a0

2a

01a0

001a

01a0

1a

100

010

001

0a

32

21

12

32

12

32

21

12


Confronto DTF - PDC

• Per la matrice A(f) si avrà:

[ ]∑=

−=p

0k

fk2jijij eka)f(A π( )

( ) ( )( ) ( )

( ) ( )⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

fAfA0

0fAfA

0fAfA

fA

3332

2221

1211

( )( ) ( ) ( ) ( )( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )⎥

⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−−−

=fAfAfAfAfAfAfAfA

0fAfAfAfA

0fAfAfAfA

)Adet(

1fH

1221221111323221

33113321

33123322

• Per la matrice H(f):

Con:


Come mai?

Confronto DTF - PDC

• Osserviamo in particolare la coppia di segnali 1-3. Per la causalità diretta da 1->3 si avrà:

( ) ( ) ( ) ( ) !0fAfAfHfDTF2

3221

2

3131 ≠==

( ) ( ) 0fAfPDC2

3131 ==

• Quindi DTF e PDC, sullo stesso set di dati, restituiscono due pattern di connettività diversi:

21

3

PDCPDC

21

3

DTFDTF


Confronto DTF - PDC

• Rivediamo le formule per un generico modello con N=3:

( ) ( ) ( )( ) ( ) ( )( ) ( ) ( )

( )( )( )

( )( )( )⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

fE

fE

fE

fX

fX

fX

fAfAfA

fAfAfA

fAfAfA

3

2

1

3

2

1

333231

232221

131211

( ) ( ) ( ) ( ) ( ) ( )( )2

2

2231322123131

fA

fAfAfAfAfHfDTF

−==

( ) ( )fAfH 1−=

( ) ( ) 23131 fAfPDC =

• Si avrà:

!0≠

Se a31[k]=0 per ogni k =>A31(f)=0

[ ]∑=

−=p

0k

fk2j3131 eka)f(A π

In cui:


Confronto DTF - PDC

• La DTFij descrive l’effetto complessivo di j su i, attraverso tutti i cammini sia diretti che indiretti (= mediati da altri canali intermedi):

21

3

=>− 31DTF )32()21( >−⋅>−+31 >−


Confronto DTF - PDC

• La PDCij descrive il solo effetto del cammino diretto (= non mediato da nessun canale intermedio) tra j e i:

=>− 31PDC 31 >−

21

3

L. Astolfi - Modelli autoregressivi per la connettività tra segnali biologici - Corso di Modelli di Sistemi Biologici 78VALUTAZIONE VALUTAZIONE PRESTAZIONIPRESTAZIONI

--

PDC/DTFPDC/DTFPDC/DTF

Valoristimati

Studio di simulazione

Segnale indipendente

(fisiologico)

Modello di connettività imposto

50 REPLICHE50 REPLICHE50 REPLICHE

BANDABANDA

a11 a12 a13

a21 a22 a23

a31 a32 a33

SNRSNR

LUNGHEZZALUNGHEZZA

GENERAZIONE SEGNALI SIMULATI

GENERAZIONE SEGNALI SIMULATI

a42

a43

a32

a31

a21

X3

X2

X1X4

GENERAZIONE INDICI TEORICI

Valori teorici


Analisi statistica dei risultati

LENGTH*SNR; Weighted Means

Current effect: F(12, 588)=1,6055, p=,08574Effective hypothesis decomposition

1 3 5 10

SNR

0,02

0,04

0,06

0,08

0,10

0,12

0,14

0,16

Rel

ativ

e E

rror

10 s 27 s 45 s 60 s 80 s

METODI MULTIVARIATI PER LA STIMA DELLA …salinari/Didattica/AutoregressiviConnettivita... ·...

Documents

Transcript of METODI MULTIVARIATI PER LA STIMA DELLA …salinari/Didattica/AutoregressiviConnettivita... ·...