1

I Sistemi PositiviX Scuola Nazionale

CIRA di Dottorato "A. Ruberti”

Bertinoro, 10-12 Luglio 2006

Ricostruzione degli ingressi nei sistemi

fisiologici

Giuseppe De Nicolao

Dipartimento di Informatica e Sistemistica, Università di Pavia

2

La misura delle serie temporali delle concentrazioni nel plasma di ormoni, substrati, farmaci, fornisce informazioni preziose per scopi conoscitivi, diagnostici, e terapeutici

Serie Temporali di Concentrazione nel Plasma

ESEMPI

3

Il segnale della concentrazione nel plasma di una sostanza y(t) è interpretabile come l’uscita di un sistema che ha come ingresso il flusso di produzione della stessa sostanza u(t)

INGRESSO u(t)

SISTEMA

USCITA y(t)

concentrazione nel plasma

(massa/volume)

secrezione(massa/tempo)

Interpretazione Causa-Effetto

4

t

d)(u)t(g)t(y

u(t)

SISTEMAg(t) y(t)

Ipotesi: sistema dinamico lineare (tempo-invariante)

Modello del Segnale di Concentrazione

secrezione concentrazione

5

k

t

kkkk vd)(u)t(gv)t(yzk

u(t)

SISTEMAg(t) y(t)

+

+ zk

vk

y(tk)

Ipotesi: sistema dinamico lineare (tempo-invariante)

Modello del Segnale di Concentrazione

k=1, 2. …n

s={t1 , t2 ,..., tn}

secrezione concentrazione serie misurata

6

Ruolo della Cinetica nel Segnale di Concentrazione

0 0.2 0.4 0.6 0.8 10

5

10

SISTEMA LENTO

mas

sa/t

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1/vo

l

0 0.2 0.4 0.6 0.8 10

0.5

1

mas

sa/v

ol

0 0.2 0.4 0.6 0.8 10

5

10

SISTEMA PRONTOm

assa

/t

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1/vo

l

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

mas

sa/v

ol

u(t) u(t)

g(t) g(t)

y(t)=g(t)*u(t) y(t)= g(t)*u(t)

Secrezione

Risposta impulsiva

Concentrazione

Oss. La velocità di g(t) determina la “memoria” del sistema

7

SISTEMA LENTOSISTEMA PRONTO

0 0.2 0.4 0.6 0.8 10

5

10

mas

sa/t

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1/vo

l

0 0.2 0.4 0.6 0.8 10

0.5

1

mas

sa/v

ol

0 0.2 0.4 0.6 0.8 10

0.5

1

t

mas

sa/v

ol

0 0.2 0.4 0.6 0.8 10

5

10

mas

sa/t

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1/vo

l

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

mas

sa/v

ol

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

t

mas

sa/v

ol

{zk} {zk}

u(t) u(t)

g(t) g(t)

y(t)=g(t)*u(t) y(t)= g(t)*u(t)

Ruolo della Cinetica nel Segnale di Concentrazione

Secrezione

Risposta impulsiva

Serie misurata

Concentrazione e campioni rumorosi

E’ importante studiare la secrezione, che però non è misurabile in vivo

8

Formulazione del Problema di Deconvoluzione

k

t

kkkk vd)(u)t(gv)t(yzk

Problema di deconvoluzione: Ricostruire l’ingresso u(t) nota la risposta impulsiva g(t) e n campioni rumorosi dell’uscita su s={t1 , t2 ,..., tn}

u(t)

SISTEMAg(t) y(t)

+

+ zk

vk

z(tk)

?

k=1, 2. …n

9

Esempi: rimozione distorsioni da un segnale o misura di un segnale per via indiretta

•Telecomunicazioni •Image restoration •Geofisica •Sismologia •Spettroscopia •Quantistica •Astronomia•Acustica•Elettromagnetismo•…

Deconvoluzione: un Problema Classico

10

Tramite il sistema http://www.pubmed.com, nel settembre 2004 la parola chiave “deconvolution” individuava 2310 riferimenti bibliografici in decine e decine di ambiti, es:

•Endocrinologia •Metabolismo •Farmacocinetica •Elettrofisiologia / Potenziali evocati •Medicina Nucleare•Radiologia / Tomografia•Medical Imaging / Ultrasuoni •Meccanica dei Materiali•Microscopia•…

Deconvoluzione per l’Analisi di Segnali Fisiologici

11

Deconvoluzione: gli elementi del problema (1 di 2)

I dati

•griglia di campionamento s={t1 , t2 ,..., tn} (uniforme o non uniforme) studiata in relazione alla dinamica attesa dei segnali in gioco

•errore di misura {vk}

E[vk] = 0 E[vk

2] = k2 costante con k o no (es. CV

costante), nota o no

Di solito {vk} è assunto gaussiano a campioni scorrelati, per cui la matrice di covarianza di v=[v1, v2, ..., vn]T:

E[vvT] = v

è diagonale. Inoltre E[v]=0

C-peptide time-series

minutes

pm

ol/m

l

0.6

1.0

1.4

1.8

0 200 400 600 800 1000 1200

C-peptide in plasma

minutes

pm

ol/m

l

0.00.30.60.91.21.51.8

0 40 80 120 160 200 240

12

Il modello della risposta impulsiva

g t A eii

Mti( )

1

Approcci per dare dei valori a M e a [Ai, i]:1) modello di popolazione2) modello individuale

2a) basato su un esperimento I/O (costoso, ma preciso)

Di solito:

es.

NB. Poiché ci può essere variabilità intraindividuale conviene fare l’esperimento I/O subito prima dell’esperimento di deconvoluzione

2b) basato su informazioni di tipo antropometriche

0 20 40 60 80 100 120 140 160 180

0

5

10

15

C-peptide decay curve

minutes

pmol

/ml

g(t)= 8.52 e-0.287t + 4.65 e-0.081t +2.16 e-0.023t

Deconvoluzione: gli elementi del problema (2 di 2)

13

Nel 1924, Hadamard diede la seguente definizione di problema ben posto

“Un problema è ben posto se esso possiede una e una sola soluzione e questa dipende con continuità dai dati. In caso contrario, viene detto mal posto”

Problemi ben-posti

14

Ipotesi (semplificative) •u(t) e' causale •u(t) costante a tratti sulla griglia di campionamento•il campionamento è uniforme

La deconvoluzione discreta

0 t1 t2 t3 t4 t5 tn-1 tn

u1 u2 u3 un

y1 y2 y3 y4 y5 yn-1 yn

T)1k(

kTk dt)t(gg~ k=0, 1, ..., n-1

k

1iikik g~uy convoluzione discreta

ti=iT i=1, 2, ….nt0=0

ku k=1, 2, ..., n

NOTA

INCOGNITA, DA STIMARE

(NB: n incognite da n dati)

15

Modello delle misure

0121n

012

01

0

g~g~g~...g~...............00g~g~g~000g~g~0000g~

G

k

k

1iikikkk vg~uvyz

z = y+ v= Gu + v

k=1, 2, ..., n

u=[u1, u2, …,un]T

z=[z1, z2, …,zn]T

v=[v1, v2, …,vn]T

y=[y1, y2, …,yn]T

G matrice di Toeplitz

16

Raw-Deconvolution

...

)ug~ug~z(g~1

ug~g~

ug~g~

g~z

u

)ug~z(g~1

ug~g~

g~z

u

g~z

u

R12

R213

0

R1

0

2R2

0

1

0

3R3

R112

0

R1

0

1

0

2R2

0

1R1

0121n

012

01

0

g~g~g~...g~...............00g~g~g~000g~g~0000g~

G

z = Gu + v

G è invertibile

zG u 1R raw-deconvolution (ignora l’esistenza di v)

Implementazione algoritmica

1+2+3+…n = n(n+1)/2 moltiplicazioni

complessità o(n2)

17

Prestazioni della raw-deconvolution

18

Problema reale

0 200 400 600 800 100050

100

150

200

250

300

raw deconvolution

minutes

pm

ol/m

in

0 200 400 600 800 10000.6

0.8

1

1.2

1.4

1.6

reconvolution vs data

minutes

pm

ol/m

l

0 200 400 600 800 1000-1

-0.5

0

0.5

1% residuals

min

%

Problema reale: secrezione insulina

19

0 20 40 60 80 100 120 140 160 180-0.2

0

0.2

0.4

0.6

0.8raw deconvolution

minutes

mU

l/min

0 20 40 60 80 100 120 140 160 1805

6

7

8

9

10

reconvolution vs data

minutes

mU

/ml

0 20 40 60 80 100 120 140 160 180-1

-0.5

0

0.5

1% residuals

min

%

Problema reale: secrezione LH

20

Necessità degli studi di simulazione

Usando solo dati reali è difficile valutare le prestazioni degli algoritmi di deconvoluzione (non posso confrontare û(t) con u(t) vero)

problemi simulati

Esempio:

problemi “fisiologici” (uso generatori di serie temporali)problemi “astratti” (uso funzioni analitiche, es. Hunt)

21

0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

1.2INGRESSO

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

1.2RISPOSTA IMPULSIVA

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 1000

0

50

100

150

200USCITA E CAMPIONI (n=41, SD=3)

t

un

ità d

i mis

ura

Il problema simulato di Hunt (1971)

1025t0 , ee)t(u

22

75600t

75400t

altrove 0502t0 1 )t(g

22

0 100 200 300 400 500 600 700 800 900 1000

-1

-0.5

0

0.5

1

1.5

2

Raw deconvolution e ingresso vero del sistema

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 10000

50

100

150

200riconvoluzione e dati

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 1000-2

-1

0

1

2residui pesati

t

res/

SD

23

0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

1.2INGRESSO

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

1.2RISPOSTA IMPULSIVA

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 1000

0

50

100

150

200USCITA e CAMPIONI (n=102, SD=3)

t

un

ità d

i mis

ura

E se ridefinissi il problema aumentando la frequenza di campionamento ?

24t

0 100 200 300 400 500 600 700 800 900 1000

-1.5

-1

-0.5

0

0.5

1

1.5

2

Raw deconvolution e ingresso vero del sistema

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 1000

0

50

100

150

200riconvoluzione e dati

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 1000-2

-1

0

1

2residui pesati

res/

SD

peggiorerebbero !!! (oscillazioni ancora più marcate)

i risultati ...

25

0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

Raw deconvolution e ingresso vero del sistema

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 1000

50

100

150

riconvoluzione e dati

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 1000-2

-1

0

1

2residui pesati

t

res/

SD

Il problema è dovuto al passaggio da tempo continuo a tempo discreto ?

NO !!!

26

La presenza di oscillazioni spurie è dovuta al fatto che

•i dati sono rumorosi

•il problema di deconvoluzione è mal-condizionato

A causa del mal-condizionamento, anche piccoli errori nei dati possono provocare grandi errori nella soluzione

Il mal-condizionamento

27

Analisi del mal-condizionamento:Ruolo del rumore di misura

zG u 1R

T1 G)vcov(G)ecov(

euvGu

)vGu(G u1

1R

uue R

dove:

è l’errore di stima. Statisticamente

E[e]=0

nella varianza dell’errore c’e’ un fattore di amplificazione

28

Il mal-condizionamento è tanto peggiore quanto più è alta la frequenza di campionamento (vd. ad es. quanto mostrato

tramite il problema simulato)

Per due problemi di deconvoluzione con frequenza di campionamento uguale, il mal-condizionamento è peggiore

per il sistema meno pronto (vd. ad es. quanto illustrato nell’esempio simulato che segue)

Valutazione qualitativa del mal-condizionamento

29

0 0.2 0.4 0.6 0.8 10

5

10

SISTEMA LENTO

mas

sa/t

u(t)

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1/vo

l

g(t)

0 0.2 0.4 0.6 0.8 10

0.5

1

mas

sa/v

ol y(t) {zk}

0 0.2 0.4 0.6 0.8 10

0.5

1

t

mas

sa/v

ol

{zk}

0 0.2 0.4 0.6 0.8 10

5

10

SISTEMA PRONTO

mas

sa/t

u(t)

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1/vo

l

g(t)

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

mas

sa/v

ol

y(t) {zk}

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

t

mas

sa/v

ol

{zk}

Mal-condizionamento: ruolo della cinetica

30

0 0.2 0.4 0.6 0.8 1-10

-5

0

5

10

15raw deconvolution e ingresso vero

t

0 0.2 0.4 0.6 0.8 10

0.5

1

1.5riconvoluzione e dati

t

0 0.2 0.4 0.6 0.8 1-5

0

5

10

15raw deconvolution e ingresso vero

t

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5riconvoluzione e dati

t

Per due problemi di deconvoluzione con frequenza di campionamento uguale, il mal-condizionamento è

peggiore per il sistema meno pronto

SISTEMA PRONTO SISTEMA LENTO

31

Raw-deconvolution: Lezioni

Credere “ciecamente” ai dati è sbagliato perché i dati sono affetti da rumore.

Il mal-condizionamento amplifica il rumore provocando nella stime dell’ingresso degli errori, in particolar modo ad alta frequenza (ovvero oscillazioni veloci e spurie)

zG u 1R

2

u

1v

T

u

R1

vuG-z minarguG-zuG-z minargu

Oss: è la soluzione del problema

La raw-deconvolution determina la stima dell’ingresso servendosi dei soli dati di misura

Ciò porta a stime contaminate da ampie oscillazioni spurie e non conformi alle nostre attese sul segnale

(vd. dimostrazione pagina seguente)

32

)uJ( minarguG-zuG-z minarguu

1v

T

u

R

uGGuuGz2zz)u(J 1v

TT1v

T1v

T

ubuAu21

)u(J TTIII

zΣGG)Σ(GbAu 1v

T11v

T1

)uG(zΣ)uG(z)uJ( 1Tv

zΣGb

GΣGA1

vT

1v

T

La funzione JIII(û) è nella forma

che si dice quadratica. Posso allora sfruttare un noto risultato della teoria del minimi quadrati per cui

Dimostrazione

Poichè il primo addendo a secondo termine non dipende da û, la û che minimizza J (û) è la stessa che minimizza

uGz2uGGu)u(J 1v

T1v

TTII

ed è la stessa che minimizza anche JIII, che differisce da JII per un fattore di scala 1/2.

uGzuGΣGu21

)u(J 1v

T1v

TTIII

Dato che G è per noi quadrata ed invertibile, si ha

zGzΣG)Σ(GGu 11v

T11v

T1

33

Il metodo di regolarizzazione

34

IDEA (Phillips, 1962), (Tikhonov, 1963): Determinare la stima operando un compromesso tra due esigenze: •spiegare i dati sperimentali•ottenere una stima dell’ingresso “regolare” (smooth)

uFFuuG-zuG-z minargu TT1T

u

PTv

aderenza ai dati regolaritàdella stima

Il metodo di regolarizzazione

11......00.........00...1100...01100001

D

F è una matrice quadrata (stessa dimensione di û) costruita in modo tale che ||Fû||2 costituisca un’indice della irregolarità della stima dell’ingresso

Tipicamente, F è costruita in modo che Fû fornisca il vettore delle differenze m-esime di û (m = 1, 2, ...). Ciò si ottiene ponendo:

F=Dm

35

F

1 0 0 0 01 1 0 0

0 1 1 00 00 1 1

...

...... ... ...... ...

m=1

1nn

23

12

1

n

3

2

1

uu...

uuuu

u

u...uuu

1.........0......10...0...1100...0110...001

uD

Per rendersi conto che Fû con F=D fornisce il vettore delle differenze prime di û è sufficiente eseguire l’operazione

Penalizzazione differenze prime

36

Per rendersi conto che Fû con F=D2 fornisce proprio il vettore delle differenze seconde di û ...

2n1nn

123

12

1

1nn

23

12

1

2

uu2u...

uu2uu2u

u

uu...

uuuu

u

1.........0......10...0...1100...0110...001

)uD(DuD

Penalizzazione differenze seconde

121...00.........00...21...0...1210...0120...001

Fm=2

37

zG)FFGG( u 1v

T1T1v

TPT

Osservazione:

Come lo stimatore ûR, anche lo stimatore PT è lineare

uFFuuG-zuG-z minargu TT1v

T

u

PT

22

u

PT uFuG-z minargu 1v

Valori piccoli di privilegeranno stime che rendono piccolo a scapito di stime che rendono piccolo 2uF

Valori grandi di privilegeranno stime che rendono piccola a scapito di stime che rendono piccola

21

vuG-z

2uF

si dice quindi parametro di regolarizzazione

21

vuG-z

38

Risultati problema Hunt con=5000 (m=2)

0 100 200 300 400 500 600 700 800 900 1000

0

0.2

0.4

0.6

0.8

1

Problem Hunt (sd=3) (ST, g=5000; q=7.1031; d=2; Tw)

time

arbi

trar

y un

its

0 100 200 300 400 500 600 700 800 900 10000

50

100

150

200reconvolution against data

time

arbi

trar

y un

its

0 100 200 300 400 500 600 700 800 900 1000-5

0

5weighted residuals

time

res/

SD

39

0 100 200 300 400 500 600 700 800 900 1000

0

0.2

0.4

0.6

0.8

1

Raw deconvolution e ingresso vero del sistema

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 10000

50

100

150

200riconvoluzione e dati

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 1000-4

-3

-2

-1

0

1

2

3

residui pesati

t

res/

SD

=1000 (m=2)

40

0 100 200 300 400 500 600 700 800 900 1000

-0.5

0

0.5

1

1.5

Raw deconvolution e ingresso vero del sistema

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 10000

50

100

150

200riconvoluzione e dati

t

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 1000

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

residui pesati

t

res/

SD

=0.1 (m=2)

41

è un parametro dell’algoritmo e la sua scelta è critica. Nella pratica non posso che affidarmi a dei criteri di regolarizzazione

42

Criterio di discrepanza (Twomey, 1965) uFFuuG-zuG-z minargu TT1

vT

u

PT

Idea (Twomey, 1965): considerare i residui

uGzr (NB: i residui dipendono da

Poiché v= y - Gu, il vettore dei residui r può essere pensato come una stima del vettore degli errori di misura

rv

Si può ritenere ragionevole che debba valere

22 vEr

Dipende da Dipende dalla varianza degli errori di misura

43

Criterio di discrepanza: Regolare in modo che

r r RSSii

n2 2

1

n

1iv

2i

n

1i

2i

n

1i

2i

T2 )(tracevEvEvvEvE

)(trace r v2

Algoritmo da usare in pratica

ripetoscelgo calcolo la stima di u per il sceltocalcolo la riconvoluzionecalcolo il vettore dei residui calcolo

finchè RSS trace v ( )

uGzr RSS ri

i

n

2

1

44

0 100 200 300 400 500 600 700 800 900 1000

0

0.2

0.4

0.6

0.8

1

Problem Hunt (sd=3) (ST, g=104.9976; q=14.4871; d=2; Twomey )

time

arb

itra

ry u

nits

0 100 200 300 400 500 600 700 800 900 10000

50

100

150

200reconvolution against data

time

arb

itra

ry u

nits

0 100 200 300 400 500 600 700 800 900 1000-2

-1

0

1

2

weighted residuals

time

res/

SD

Risultati Hunt con =104 (discrepanza)

45

0 200 400 600 800 100050

100

150

200

250

300

deconvoluted ISR (ST, g=0.0014525; q=64.1331; d=2; Tw

minutes

pmol

/min

0 200 400 600 800 10000.6

0.8

1

1.2

1.4

1.6

reconvolution vs C-peptide data

minutes

pmol

/ml

0 200 400 600 800 1000

-8

-6

-4

-2

0

2

4

6

% residuals

min

%

Risultati Insulina

46

0 20 40 60 80 100 120 140 160 180-0.2

0

0.2

0.4

0.6

0.8DECONVOLUZIONE(g=129.802; q(g)=18.3146; d=1; Tw)

minutes

mU

/min

0 20 40 60 80 100 120 140 160 1805

6

7

8

9

10

RICONVOLUZIONE E DATI

minutes

mU

/ml

0 20 40 60 80 100 120 140 160 180-10

-5

0

5

RESIDUI %

min

%Risultati LH

47

Il criterio di discrepanza di Twomey è intuitivo ma ha una base teorica debole

Si può dimostrare che, in media, il criterio di Twomey conduce ad oversmoothing ( troppo grande)

Esistono in letteratura altri criteri di regolarizzazione, ma più complicati, ad esempio:

•Maximum Likelihood (ML, massima verosimiglianza)•Minimum Risk (rischio minimo)•Generalized Cross-Validation

Il criterio di discrepanza risolve il problema ?

inclusi nel software WINSTODEC

48

Regolarizzazione: problemi ancora aperti

•campionamento rado / non uniforme

0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

Raw deconvolution e ingresso vero del sistema

un

ità d

i mis

ura

0 100 200 300 400 500 600 700 800 900 1000

0

0.2

0.4

0.6

0.8

1

Problem Hunt (sd=3) (ST, g=48.2675; q=16.9019; d=2; Max Lik1)

time

arb

itra

ry u

nits

•vincolo non negatività

49

Campionamento non uniforme / radoLa griglia virtuale

PROBLEMA: essendo partiti dalla raw-deconvolution, abbiamo sempre considerato griglia di campionamento = griglia delle stime

IDEA: Disaccoppiamo la griglia di campionamento (dimensione n) da una griglia "virtuale" (dimensione N>>n) su cui discretizziamo l'ingresso

ESEMPIO:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

griglia “virtuale” (1 min)

griglia campionamento (non uniforme e rada)

1 ? 3 ? ? 6 ? ? ? 10 ? ? ? ? 15

15

4

3

2

1

01121314

0789

0345

012

0

u..............uuuu

g~g~............g~g~g~00g~.........g~g~g~0000g~...g~g~g~000000g~g~g~00000000g~

)15(y)10(y)6(y)3(y)1(y

50

15

14

4

3

2

1

01121314

01213

0123

012

01

0

uu..........uuuu

g~g~g~g~g~0g~g~g~

............000...0g~g~g~g~000000g~g~g~0000000g~g~00000000g~

)15(y)14(y

...

.......

)4(y)3(y)2(y)1(y

15

14

4

3

2

1

01121314

01213

0123

012

01

0

uu..........uuuu

g~g~g~g~g~0g~g~g~

............000...0g~g~g~g~000000g~g~g~0000000g~g~00000000g~

)15(y)14(y

...

.......

)4(y)3(y)2(y)1(y

15

4

3

2

1

01121314

0789

0345

012

0

u..............uuuu

g~g~............g~g~g~00g~.........g~g~g~0000g~...g~g~g~000000g~g~g~00000000g~

)15(y)10(y)6(y)3(y)1(y

In dettaglioSe avessi tutti i campioni sulla griglia virtuale...

Eliminando le righe corrispondenti ai dati “mancanti” ...

Si ottiene ...

51

z = Gu + v

Soluzione in presenza di griglia virtuale

nX1 nXN

NX1

nX1

= +z G u v

zG)FFGG( u 1-v

T1T1-v

T

uFFuuG-zuG-z minargu TT1v

T

u

Soluzione in forma chiusa

invertibile perchè definita positiva

GG 1-v

TFFT

semi-definita positiva

definita positiva

52

0 100 200 300 400 500 600 700 800 900 1000

0

0.2

0.4

0.6

0.8

1

DECONVOLUZIONE ( =1237454.664; q()=13.8563; d=2; Twomey )

time

arb

itra

ry u

nits

0 100 200 300 400 500 600 700 800 900 10000

50

100

150

200RICONVOLUZIONE E DATI

time

arb

itra

ry u

nits

0 100 200 300 400 500 600 700 800 900 1000

-2

-1

0

1

2

RESIDUI PESATI

time

res/

SD

Risultati Hunt con griglia virtuale

53

0 200 400 600 800 100050

100

150

200

250

300

regularized deconvolution (=0.58821, q()=34.07, n=115, Twomey)

minutes

pmol

/min

0 200 400 600 800 10000.6

0.8

1

1.2

1.4

1.6

reconvolution vs data

minutes

pmol

/ml

0 200 400 600 800 1000

-10

-5

0

5

10

% residuals

min

%

Risultati insulina con griglia virtuale

54

Risultati LH con griglia virtuale

0 20 40 60 80 100 120 140 160 180-0.2

0

0.2

0.4

0.6

0.8DECONVOLUZIONE(g=775.619; q(g)=18.4213; d=1; Tw)

minutes

mU

/min

0 20 40 60 80 100 120 140 160 1805

6

7

8

9

10

RICONVOLUZIONE E DATI

minutes

mU

/ml

0 20 40 60 80 100 120 140 160 180-10

-5

0

5

RESIDUI %

min

%

55

Vincolo di non negatività

uFFu)uGz()uGz(minargu TT1v

T

0u

PTV

La regolarizzazione "lineare" puo' produrre stime negative

IDEA: Vincolare il problema

•Per uso un valore che funziona bene (quantità di regolarizzazione adeguata) nel caso lineare•Come algoritmo numerico per risolvere il problema di minimo vincolato utilizzo un metodo iterativo

Problemi

La stima diventa non lineare:•Scelta di ?•Algoritmi ?

Rimedi

56

0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

DECONVOLUZIONE CON VINCOLO DI NON NEGATIVITA’ vs VERO

time

Uni

tà d

i mis

ura

0 100 200 300 400 500 600 700 800 900 10000

50

100

150

200RICONVOLUZIONE e DATI

time

Uni

tà d

i mis

ura

0 100 200 300 400 500 600 700 800 900 1000

-6

-4

-2

0

2

4

6

RESIDUI PESATI

time

res/

SD