CORSO DI MODELLI DI SISTEMI BIOLOGICI LAUREA IN INGEGNERIA CLINICA E BIOMEDICA

1. CORSO DI MODELLI DI SISTEMI BIOLOGICI LAUREA IN INGEGNERIA CLINICA E BIOMEDICA

2. Cenni sulla stima parametrica Definizione del problema Nel processo reale di stima di parametri si presentano quattro diversi tipi di errori: imperfezioni nei segnali di test, errori di misura, disturbi, errori della struttura del modello. Le imperfezioni nei segnali di test sono di norma note e dunque possono essere tenute debitamente in conto o trascurate se non significative. Gli eventuali errori presenti nella struttura del modello non possono essere considerati esplicitamente durante la fase di stima dei parametri (in quanto non noti a priori). Sarà in fase di validazione del modello che l’appropriatezza di una data struttura verrà valutata. In questa sede si analizzeranno esplicitamente gli effetti dovuti al rumore e agli errori di misura inevitabilmente presenti in un processo di acquisizione dati. Nei casi più semplici è sufficiente assumere l’ipotesi di rumore additivo, tale ipotesi dovrà essere rimossa in caso di situazioni sperimentali più complesse.

3. Definizione del problema Dato il modello defnito dalle equazioni: dx(t,q)/dt = f[x(t,q), u(t),q, t] xo = x(to,q) y(t,q) = g[x(t,q),q] h[x(t,q), u(t),q] = 0 Siano: zl(tk, q) = yl (tk, q)+el(tk) per l=1,2,……..,m e per k=1,2,…….,N Le componenti del vettore di misura. Il problema consiste nello stimare il vettore incognito dei parametri q attraverso i dati corrispondenti alle variabili accessibili e i rispettivi valori forniti dal modello. Si noti che il valore misurato zl(tk,q) dell’uscita l-esima all’istante tk è affetta da rumore additivo el(tk). Il caso di misura tempo discreto è dovuto al fatto che non sempre è possibile misurare alcune variabili in modo continuo e che comunque l’applicazione degli algoritmi di stima implica in generale un campionamento delle misure. Per motivi di semplicità si considererà nel seguito un sistema a una sola uscita: y(t,q) = g[x(t,q), q] z(tk,, q) = y (tk, q)+e(tk) per k=1,2,…….,N

4. STIMATORI • Uno stimatore è essenzialmente una funzione opportuna che lega i dati misurati alla stima di una quantità incognita. • Sia z è il vettore dei dati ottenuto misurando l’uscita in N istanti di tempo z=[z(t1), z(t2), …….., z(tN)] • Si definisce stimatore di q un vettore a p dimensioni di funzioni della variabile aleatoria Z così definito:(z) = [1(z), 2(z), …………, p(z)]T • Dove ciascuna funzione è definita nel dominio RN e assume valori nello spazio ammissibile Rp i del parametro qii : RN Rqi per i=1,2,…….,p • ed è tale che tutte le i soddisfino simultaneamente le seguenti condizioni: • 1.   Le funzioni i devono essere non polarizzate, cioè il valore atteso della stima deve coincidere con il valore vero • 2.  Le funzioni i devono essere tali che la varianza della stima ottenuta con lo stimatore i (2i(qi)) sia inferiore alla varianza delle stime ottenute con ogni altro stimatore non polarizzato. • 3.   Le funzioni i devono essere tali che la corrispondente matrice di covarianza i cui elementi sono dati da hk = E[h(z) – Eh(z)]·[ k(z) – Ek(z)] deve essere uguale all’inversa della matrice d’informazione di Fisher che, come stabilito dal teorema di Cramer-Rao, coincide con il più piccolo valore teoricamente ottenibile per la covarianza. • 4.  Le funzioni i devono essere consistenti, cioè tali che al crescere del numero di campioni disponibili forniscano una stima sempre più vicina al valore vero. • Il valore deterministico fornito dalla (z) è detto stima di q basata sui dati z.

5. STIMA DI MASSIMA VEROSIMIGLIANZA Il procedimento di stima di massima verosimiglianza richiede un modello che descriva i dati sperimentali, un modello del rumore e la definizione della funzione di verosimiglianza che descrive la probabilità di generazione dei dati sperimentali. Considerando sempre il caso di una sola uscita, sia L(z,q) la funzione di verosimiglianza, il principio del metodo della massima verosimiglianza consiste nel considerare come stima diqquel valore qche rende a posteriori più probabile l’osservazione dei dati z. Cioè deve essere: L[z(t1), z(t2), ……., z(tN), q)  L[z(t1), z(t2), ……., z(tN), q) q Questo problema di massimizzazione è di difficile soluzione e inoltre in molti casi la funzione di verosimiglianza non è nota. Tuttavia la stima di massima verosimiglianza riveste un’importanza fondamentale nelle applicazioni pratiche in quanto essa gode delle quattro proprietà precedentemente pronunciate. Inoltre, nell’ipotesi di campioni di rumore statisticamente indipendenti e normalmente distribuiti, essa coincide con la stima di Markov e se il rumore è bianco coincide con la stima ai minimi quadrati.

6. STIMA LINEARE AI MINIMI QUADRATI • Si faccia riferimento ad un modello a una sola uscita in cui questa dipenda linearmente dai parametri: • y(t,q) = g(x(t))q =  gi(x(t))qi =  ri(t) qi • L’obiettivo del processo di stima è di ottenere i valori numerici di q a partire dai dati rumorosi z(ti) e dalla forma funzionale del modello. La relazione esistente tra la risposta del modello e i dati sperimentali è data da: • z(tl, q) = y(tl, q) + e(tl) l= 1, 2, 3, .... , N   • dove le proprietà statistiche dell’errore sono descritte dal suo valore atteso e dalla matrice di varianza e covarianza. • L’esistenza dello stimatore generalizzato ai minimi quadrati è basata su quattro ipotesi fondamentali: • La variabile y è una variabile aleatoria il cui valore atteso per ogni istante t è descritto da una relazione lineare rispetto a q •   Le variabili indipendenti ri(t) sono quantità deterministiche •  Per ogni vettore rk=(r1(tk), r2(tk), ………., rp(tk)) i corrispondenti valori sono indipendenti. •   Gli errori e(tk) sono statisticamente indipendenti, normalmente distribuiti, a media nulla con matrice di varianza covarianza G

7. STIMA LINEARE AI MINIMI QUADRATI Utilizzando una notazione matriciale si potrà quindi scrivere: e = z-y = z-R·q con Il procedimento di stima ai minimi quadrati generalizzato consiste nell’impostare una funzione di costo J da minimizzare data da: J = e’·W·e dove Wè una matrice peso non singolare diagonale. Sostituendo l’espressione dell’errore e derivando rispetto a q si ottiene con ovvi passaggi il valore di q che minimizza l’indice J q = [RTWR]-1rTWzV(q) = [RTWR]-1RTW G W R [RTW R]-1 Con V(q) matrice di varianza e covarianza della stima Dalla particolare scelta della matrice peso derivano due stimatori di importanza pratica, W = I Stimatore ai minimi quadrati: qL = [RTR]1RTzV(qL) = [RTR]-1RTGR [RTR]-1 W = G-1Stimatore di Markov.: qM = [RTG-1R]1RTG-1zV(qM) = [RTG-1R]-1 r1(t1)...rp(t1) ..................... r1(tN)...rp(tN) R =

8. STIMA LINEARE AI MINIMI QUADRATI Per dimostrare che lo stimatore di Markov è uno stimatore di massima verosimiglianza sotto l’ipotesi di errori normalmente distribuiti, si consideri la funzione di massima verosimiglianza che nel caso considerato assume la forma :   L(z,q) = exp ( (z – Rq)T·G-1·(z – Rq ))/ (2 (2p)N/2G1/2) Il valore q che si ottiene utilizzando il criterio di stima alla massima verosimiglianza è per definizione quel valore per cui la funzione di verosimiglianza è massima. Ne consegue che massimizzare rispetto a qequivale a massimizzare il suo logaritmo: ln[L] = c – [(z – Rq)T·G-1·(z – Rq)]/2 dove c è una costante indipendente dal valore dei parametri. Per risolvere il problema di stima occorre quindi minimizzare la quantità: F(q) = (z – Rq)T·G-1·(z – Rq) che conduce ad un valore del vettore ottimo dei parametri pari a quello fornito dallo stimatore di Markov. Con argomentazioni analoghe si può dimostrare che, in presenza di un rumore gaussiano bianco lo stimatore di Markov è ancora uno stimatore alla massima verosimiglianza e coincide con lo stimatore ai minimi quadrati. Questa circostanza ci permette di usare lo stimatore ai minimi quadrati propriamente detto, ogni volta che si è in presenza di rumore gaussiano bianco.

9. STIMA LINEARE AI MINIMI QUADRATI • Nella pratica, purtroppo, la matrice di covarianza dell’errore è generalmente incognita e ciò che è possibile fare è soltanto avanzare delle ipotesi a priori sulle sue proprietà, ipotesi che dovranno essere verificate a posteriori. I casi pratici che più frequentemente si presentano durante i processi di misura sono sostanzialmente due: • Si è in presenza di un rumore in cui i campioni sono incorrelati, a media nulla, con varianze diverse ma note a meno di un fattore di scala s2 (G = W*s2). Si può dimostrare in modo analogo a quanto già fatto, che in queste condizioni lo stimatore che si ottiene ponendo W = W* è ancora uno stimatore alla massima verosimiglianza. • Si è in presenza di un rumore bianco e stazionario. In questo caso la matrice di varianza e covarianza G sarà una matrice diagonale con gli elementi della diagonale principale uguali tra loro (G = s2·I) cons2 fattore di scala incognito. Si ricade quindi in un caso particolare dell’ipotesi precedente, ne consegue quindi che in queste condizioni lo stimatore che si ottiene ponendo W = I è ancora uno stimatore alla massima verosimiglianza.

10. STIMA NON LINEARE AI MINIMI QUADRATI I modelli più frequentemente usati hanno, in generale, uscite non lineari nei parametri. Con riferimento al caso particolare di una sola uscita si avrà y(t)= g(x(t),q) Per quanto detto nel paragrafo precedente, se le informazioni sul rumore sono tali da far supporre la presenza di rumore gaussiano, è corretto applicare una procedura di stima ai minimi quadrati al problema in questione. Il problema della sima dei parametri è così ricondotto alla minimizzazione della solita funzione costo dove in questo caso risulta e = z – g(x(t),q) Poiché il modello è non lineare ai parametri, la funzione di costo non sarà in generale quadratica nei parametri e quindi non sarà possibile, nella maggioranza dei casi, ottenere una soluzione analitica esplicita del problema di ottimizzazione. Sarà quindi necessario adottare un approccio di tipo iterativo.

11. Stima di parametri di popolazione Esistono vari metodi per la stima dei parametri di una popolazione. Questi metodi ottengono tali stime mediante la minimizzazione di una funzione obiettivo, derivata dall’applicazione dei minimi quadrati pesati o estesi a sviluppi di Taylor dell’equazione di uscita del modello arrestati al primo o al secondo termine (NONMEM). Metodi più recenti (NPML) forniscono una stima per la distribuzione di probabilità dei parametri da cui è possibile stimare le grandezze statistiche di interesse.

12. IL METODO NPD ( Naive Pooled Data ) Il principio su cui è basato tale metodo è quello di considerare tutti i dati di popolazione come provenienti da un unico individuo: da qui il nome del metodo: Metodo “ingenuo” dei dati raggruppati. Si consideri per semplicità il caso in cui l’uscita del sistema sia monodimensionale e sia yi il valore di tale uscita misurato nell’istante di tempo ti. Sia inoltre N il numero complessivo delle osservazioni su tutti gli individui del campione. Se ymodi = h(x(ti);q) è l’uscita del modello calcolata nello stesso istante ti in funzione dello stato e dei parametri, la funzione obiettivo ai minimi quadrati per tale tipo di approccio sarà: J(q) = i=1...N [yi-h(x(ti;q)]2 Il minimo della J rispetto a q fornisce la stima q* dei valori medi dei parametri di popolazione. La matrice di covarianza della stima può essere ottenuta in modo approssimato calcolando l’inversa della matrice delle derivate parziali del secondo ordine di J(q) rispetto a qvalutata in q*.Gli errori standard sono proporzionali agli elementi della diagonale principale dell’inversa di tale matrice. Il metodo NPD non fornisce stime relative alle varianze dei parametri nella popolazione.

13. IL METODO TS (Two-Stage) In tale metodo i dati relativi a ciascun individuo sono analizzati separatamente (Primo Stadio), ottenendo così stime individuali dei parametri. Nel Secondo Stadio, da tali stime si determinano le caratteristiche statistiche della popolazione (Media e Varianza dei parametri sulla popolazione). Per i parametri individuali si utilizza in genere una stima ai minimi quadrati pesata con la varianza dell’errore di misura. Se si hanno M individui nella popolazione e sono p i parametri da stimare, al termine del primo stadio si avranno p valori dei parametri ottimi per ognuno degli M individui: q*kj con j = 1...M e k = 1...p . Le stime dei parametri di popolazione si ottengono mediante le relazioni che forniscono media e varianza del campione. qkjMedio = 1/M qkjs2qk = 1/(M-1)  (qkj-qkjMedio )2 con j = 1,..., M e k = 1, ..., p. I problemi principali di questo metodo derivano dalla scarsa affidabilità delle stime individuali. Se queste sono errate (ad es. a causa delle poche misure disponibili)gli errori saranno elevati anche nelle stime di popolazione. Inoltre la correlazione degli errori sulle stime individuali potrà simulare una correlazione tra soggetti della popolazione anche quando questa non sussiste

14. IL METODO NONMEM (Non Linear Mixed Effect Model) Il metodo NONMEM proposto nel 1980 da Sheiner e Bealfornisce stime accurate dei parametri di popolazione e si basa sulle seguenti assunzioni: Si consideri per l’individuo j-esimo della popolazione (j = 1, ..., M), il modello statistico: yj = f(tij, xj; qj) + eij in cui i = 1, ... Nj sono le misure sul soggetto j. Il vettore dei parametri risulterà essere una variabile aleatoria, supposta avente la stessa distribuzione per ogni individuo e normale, con valore atteso qmedio e matrice di covarianza Sq : (qj  N(qmedio, Sq) Sviluppando in serie di Taylor intorno a qmedio e arrestandosi al primo termine si ottiene: yij = f(tij, xij; qmedio) +Sk=1,...,p f/ qktij;qmedio(qkj – qmedio) + eij Indicando con hj = (qj – qmedio) si avrà in forma più compatta: yj = f(qmedio,xj) + Gj(qmedio,xj)hj + ej Ove Gj è la matrice Jacobiana f/qktij;qmedio di dimensioni Nj x p. Per ottenere un’espressione della cov(yj), poichè la Gj + ej è combinazione lineare di vettori aleatori, gaussiani e a media nulla, la cov(yj) avrà media nulla e covarianza data da: cov(yj) =Gj(qmedio,xj) Sq Gj(qmedio,xj) T + se2 INj Con INj matrice identità di dimensione Nj e se2 varianza dell’errore di misura.

15. IL METODO NONMEM (Non Linear Mixed Effect Model) Si osserviche le covarianze sulle osservazioni, relative ad un dato individuo, cioè gli elementi fuori della diagonale principale di cov(yj) sono non nulle , e sono anche funzioni della variabilità interindividuale. Il NONMEM riconosce quindi che le componenti di ogni vettore yj devono essere correlate le une alle altre tramite il set comune di differenze dei parametri di ciascun individuo dai valori medi di popolazione. La funzione obiettivo fornita dall’applicazione della massima verosimiglianza è quindi: J(qmedio, Sq ) =  j=1,...,M[ lncov(yj) + (yj - f(qmedio,xj) )T (cov(yj))-1 (yj - f(qmedio,xj) )] Dove cov(yj) indica il determinante di cov(yj). Se il numero di parametri è elevato ed esiste possibilità di correlazione fra essi ed il numero dei soggetti è ridotto, si può utilizzare un approccio in cui si stima il valore medio dei parametri sulla popolazione , costruendo la matrice di covarianza a partire dai dati (Metodo MACOV) . Assumendo che per ognuno degli M soggetti si abbiano N misure (N uguale per tutti i soggetti) si determina il vettore m = [ Sj=1,...,M y1j/M, ..., Sj=1,...,M yNj/M]Te il vettore Y = [Y1 -m ......YN -m] . Una stima della matrice di covarianza delle misure è quindi fornita dalla matrice NxN:  =1/(M-1) YTY. La funzione obiettivo da minimizzare rispetto a qmedio è: J(qmedio) = Sj=1,...,M (yj - f(qmedio )T  -1(yj - f(qmedio )