Indice 1 Il metodo della massima verosimiglianza 1.1 1.2 1.3 1.4 1.5 1.6 1.7 Denizione del problema . . . . . . . . . . . . . Variabile campionaria - Statistiche - Stimatori . Correttezza . . . . . . . . . . . . . . . . . . . . Consistenza . . . . . . . . . . . . . . . . . . . . Su cienza - Su cienza minimale - Completezza Stimatori di minima varianza: e cienza . . . . Stime di massima verosimiglianza . . . . . . . . 2 Il metodo dei minimi quadrati . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Formulazione generale del problema (caso lineare) . . . . 2.3 Soluzione del problema di minimi quadrati: stimatori di osservabili e parametri . . . . . . . . . . . . . . . . . . . 2.4 Covarianza degli stimatori e stima di o2 . . . . . . . . . 2.5 Ottimalita degli stimatori m.q.: Teorema di Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Problemi di minimi quadrati con vincoli . . . . . . . . . 2.7 Problemi di stima non lineari . . . . . . . . . . . . . . . 2.8 Applicazione alla regressione lineare . . . . . . . . . . . . 1 1 4 6 12 17 23 35 42 42 47 51 56 62 65 73 81 1 Il metodo della massima verosimiglianza 1.1 Denizione del problema Nel primo quaderno abbiamo denito in modo parallelo le variabili casuali (v.c.) e le variabili statistiche (v.s.) spiegando come queste si comportino in modo formalmente identico una volta che si sia istitui1 ta la corrispondenza frequenze $ probabilita e si sia conseguentemente denito l'operatore di media E fg. Elemento centrale di quella corrispondenza e stato il considerare una v.s. come il riordino di una \popolazione" di valori fx1 : : : xN g, ottenuti ripetendo N volte l'esperimento stocastico E di cui una certa v.c. X descrive il comportamento aleatorio. Diremo in questo caso che fx1 : : : xN g costituisce un campione di tipo bernoulliano tratto dalla v.c. X . L'aggettivo bernoulliano si riferisce all'ipotesi, che qui supporremo sempre soddisfatta, che le ripetizioni dell'esperimento E siano tali da non inuenzarsi stocasticamente tra loro. Osservazione 1.1.1: per comprendere che vi sono casi signicativi in cui invece le ripetizioni si inuenzano reciprocamente, basta pensare all'esempio dell'estrazione di due numeri da un'urna che ne contenga 3, senza che il primo estratto venga riposto nell'urna stessa prima della seconda estrazione. In questo caso la prima estrazione e descritta dalla v.c. X1 11=3 12=3 13=3 ammesso che il primo estratto sia ad esempio x1 = 1, la seconda estrazione ha la distribuzione X2 12=2 13=2 con le ovvie modiche nel caso il primo estratto fosse uno degli altri due numeri. Dunque in questo caso la distribuzione di X2 dipende chiaramente dal valore assunto da X1 . Spesso nell'esperimento E si ha una conoscenza qualitativa del meccanismo stocastico che genera il comportamento aleatorio e poi appunto se ne conosce un campione (bernoulliano) fx1 : : : xN g. La conoscenza qualitativa di E puo allora suggerire che la v.c. X che lo descrive, abbia una distribuzione appartenente ad una certa famiglia fX (x ), dove rappresenta il parametro (mono- o pluri-dimensionale) che specica la particolare densita di probabilita che ci interessa: supporremo che possa assumere valori in un insieme specicato. 2 Qualora questo elemento qualitativo non ci sia, talvolta, se il campione e numeroso, si puo prendere in esame l'istogramma e dalla sua forma arguire l'appartenenza di X ad una certa famiglia. Esempio 1.1.1: sia E il lancio di una moneta: poiche i valori argomentali di X sono solo due, questa v.c. non puo che essere rappresentata da una famiglia ad un parametro X 0p 1 ;1 p (1.1.1) =p: (1.1.2) X = N 2 ] (1.1.3) = 2 : (1.1.4) dove chiaramente analogo ragionamento vale per ogni E che ammetta un numero nito di valori argomentali. Esempio 1.1.2: sia E una misura di precisione di una grandezza continua in tal caso si potra supporre e quindi si potra porre Analogo ragionamento vale per ogni v.c. per cui si possa invocare il teorema centrale della statistica. Il problema che ci poniamo ora e il seguente: supponendo di conoscere la famiglia fX (x ), cui la distribuzione della v.c. X appartiene per un certo = , e conoscendo un campione (bernoulliano) fx1 : : : xN g estratto da X , dare una stima t di , intendendo con cio per il momento almeno che jt ; j e piccolo, con elevata probabilita, rispetto ad un criterio pressato. 3 Poiche t deve essere un numero (od un vettore di numeri) calcolabile sulla base delle conoscenze che abbiamo, risulta ovvio che t, oltre che dalla forma della famiglia f (x ), dovra dipendere dai valori del campione, cioe t = t(x1 : : : xN ) : (1.1.5) Il nostro problema quindi sara di discutere i criteri di scelta della funzione (1.1.5) in base alle proprieta che ne derivano. 1.2 Variabile campionaria - Statistiche - Stimatori La variabile campionaria nel caso bernoulliano e la variabile che descrive N repliche identiche e indipendenti dell'esperimento E , ovvero descrive l'esperimento complesso E N E : : : E . Se la v.c. X e ad 1 dimensione, allora la v.c. campionaria sara N dimensionale X (N ) , e fatta in modo tale che ogni sua componente abbia la stessa distribuzione di X X1 X X (N ) = .. 2 . XN fX (x) = fX (x ) : i (1.2.1) La distribuzione della v. campionaria X (N ) e poi interamente identicata tramite l'ipotesi che le sue componenti siano tutte indipendenti tra loro, cos che fX (x ) = fX (x1 : : : xN ) = Ni=1 fX (xi ) : (1.2.2) La funzione di densita di X (N ) e detta \likelihood" (verosimiglianza) e noi la indicheremo come L(x ) = fX (x ) : (1.2.3) L'importanza della v. campionaria sta nel fatto che un campione ber4 noulliano fx1 : : : xN g di N estrazioni da X puo essere considerato come una singola estrazione della variabile X (N ) . Un esempio chiarira meglio il concetto. Esempio 1.2.1: si consideri il classico gioco a testa e croce (E ) un giocatore lancia due volte una moneta ottenendo il campione (t c). Ma questo stesso puo essere anche considerato come un campione di una sola replica di E 2, cioe di un esperimento costituito da due lanciatori indipendenti che lanciano contemporaneamente due monete identiche. Osservazione 1.2.1: il concetto di variabile campionaria puo essere applicato anche al caso (non bernoulliano) in cui la ripetizione dell'esperimento E non e fatta nelle identiche condizioni, cos che ogni ripetizione sara descritta da una diversa v.c. Xi. Un esempio importante di questo tipo e costituito dal caso in cui una stessa quantita e misurata in condizioni diverse o con strumenti diversi, cos che ogni misura avra una diversa varianza i2 . In questo caso si usera semplicemente la denizione (1.2.3), tenendo pero conto che ora fX (x ) = Ni=1 fXi (xi ) (1.2.4) senza che le fXi siano uguali tra loro. Anzi, per questa strada si potrebbe anche lasciar cadere l'ipotesi di indipendenza stocastica delle Xi, denendo la v. campionaria semplicemente come quella variabile che descrive il complesso degli N esperimenti considerati, indipendenti o no che siano le variabili che li descrivono, e la funzione di likelihood come la densita di probabilita di tale variabile N -pla. Osservazione 1.2.2: nel caso di variabili discrete (ad es. binomiale, poissoniana, ecc.) deniamo come likelihood in un punto x1 = k(1) x2 = = k(2) : : : direttamente il prodotto delle probabilita che ogni componente xi assuma il valore assegnato k(i). Percio se PX fx = k g = pk k = 1 2 : : : si ha LX (x) = Ni=1 PX fxi = k(i)g : 5 (1.2.5) Deniamo statistica S (X (N ) ) una qualunque funzione, a una o piu dimensioni, della v. campionaria X (N ) . Cos ad esempio la v.c. M = N1 X Xi N i=1 (1.2.6) e una statistica che prende ovviamente il nome di media campionaria. E bene osservare che (1.2.6), come statistica e una v.c. e non gia una costante: e chiaro pero che per un certo campione dato fx1 : : : xN g, cioe per una estrazione da X (N ) , si potra calcolare la corrispondente media campionaria, cioe il numero N X 1 m = N xi i=1 (1.2.7) e che tale numero puo essere considerato come un'estrazione da M. Alla luce di queste denizioni possiamo riprendere il problema della stima denito nel x1.1 nel seguente modo: data una certa funzione di likelihood L(x ), cerchiamo un'opportuna statistica T (X (N ) ), con una distribuzione su cientemente concentrata attorno al valore , in modo tale che un'estrazione t da T , calcolata in corrispondenza a un campione estratto x t = T (x), abbia un'elevata probabilita di essere vicina a . Con questa impostazione T si chiama uno stimatore di mentre la sua estrazione t = T (x) si chiama stima. Vedremo nei prossimi paragra le proprieta di maggior interesse degli stimatori, in particolare studieremo la media di T (X (N ) ) e la sua varianza per N ! 1. 1.3 Correttezza Per semplicita rappresentiamo in questo paragrafo la v. campionaria col vettore (X ). Avendo stabilito che uno stimatore T (X ) debba essere distribuito attorno 6 a , la prima questione interessante sara di vedere quale relazione vi sia tra la media di T (X ), parametro di posizione per eccellenza, e . Deniamo il bias (distorsione) di T (X ) come b() = E fT (X )jg ; : (1.3.1) Questa formula richiede un commento: intanto intendiamo con E fT (X )jg la media di T (X ), facendo l'ipotesi che sia il valore corretto del parametro, ovvero E fT (X )jg = Z RN T (x)L(x )dN x : (1.3.2) Dunque, supposto che sia il valore giusto del parametro, b() misura la deviazione della media di T (X ) da . Uno stimatore che abbia bias nullo E fT (X )jg = (1.3.3) si dice unbiased o corretto. Esempio 1.3.1: la media campionaria M = (1=N ) Pi Xi e uno stimatore corretto di , qualunque sia la distribuzione sottostante. Infatti notiamo che per ogni componente Xi si ha 1 E fXijg = (1.3.4) X X E fMjg = N1 E fXi j g = N1 = : (1.3.5) cos che i i Esempio 1.3.2: i momenti (non centrali) campionari di ordine n qualsiasi sono stime unbiased dei corrispondenti momenti della distribuzione sottostante. Sia Cio esprime l'ovvio fatto che ogni componente e di per se uno stimatore unbiased di . 1 7 M n = N1 X Xin ( ) (1.3.6) i la variabile campionaria \momento campionario di ordine n". Poiche chiaramente E fXinjng = n ragionando come per la media si ha (1.3.7) E fM(n)jng = n : (1.3.8) Esempio 1.3.3: la varianza campionaria non e uno stimatore corretto della varianza teorica della popolazione sottostante. Deniamo come varianza campionaria la variabile S 2 X = N1 (Xi ; M)2 = X 2 X ; (M)2 = M(2) ; (M)2 : = 1 i N (1.3.9) Studiamo la media di M2. Si ha tenendo conto che le diverse componenti di X sono tra loro indipendenti, E fM2g ( ) ( ) X = E 12 XiXk = N ik ( ) X X = E 12 XiXk + E 12 Xi2 = N i=k N i 2 = N ;2 N 2 + N22 = (1.3.10) N N = N ; 1 2 + 1 2 : N N 6 Pertanto, tornando alla (1.3.9) si ha 8 E fS 2 g = 2 ; N N; 1 2 ; N1 2 = = N N; 1 2 ; 2 = N N; 1 2 : (1.3.11) Osservazione 1.3.1: nell'esempio (1.3.3) si e dimostrato che la media di S non coincide con 2 , ma e una funzione lineare di 2. In questo caso e particolarmente semplice modicare lo stimatore in modo da ottenerne un altro che sia unbiased. Bastera infatti porre 2 S 2 X = N S2 = 1 (Xi ; M)2 N ;1 N ;1 i (1.3.12) perche si abbia E fS 2g = 2 : (1.3.13) Lo stimatore (1.3.12) e detto varianza campionaria corretta: come si vede esso di!erisce poco (anche numericamente) da S 2 quando N sia grande, mentre la di!erenza puo essere signicativa quando N sia pari a poche unita. Esempio 1.3.4: sia X una variabile doppia X = VU e sia X (N ) la corrispondente variabile campionaria (2N -dimensionale) le cui componenti Xi = UVi sono tra loro indipendenti e identicamente distribuite. La i covarianza campionaria SUV X = N1 (Ui ; MU )(Vi ; MV ) X = N1 Ui Vi ; MU MV (1.3.14) e uno stimatore non corretto di UV . Poiche tanto SUV quanto UV sono indici invarianti per traslazione, possiamo supporre di metterci nella situazione piu comoda, cioe U = 0 V = 0, senza alterare il risultato. 9 Calcoliamo prima ( ) X E fMU MV g = E N12 UiVk = ik X = N12 E fUiVig = i N = NUV = NUV 2 poiche E fUi Vk g = 0 (i 6= k) in base alle ipotesi fatte. Pertanto, dalla (1.3.14) = N ; 1 UV : (1.3.15) N Analogamente a quanto si fa per la varianza, e anche qui semplice correggere lo stimatore denendo una covarianza corretta tramite la formula X E fSUV g = N1 UV ; N1 UV S UV = N 1; 1 X(Ui ; MU )(Vi ; MV ) : (1.3.16) Osservazione 1.3.2: nella denizione (1.3.3) di correttezza si e supposto che E fT g fosse uguale al parametro per tutti i suoi valori: in questo caso si dira che T e uno stimatore uniformemente corretto di . Vi sono casi invece in cui lo stimatore T (X ) puo essere corretto solo per qualche valore di . Si consideri ad esempio il coe ciente di correlazione campionario denito da R = SSUVS : U V (1.3.17) supponiamo che la variabile doppia sottostante sia X = VU una normale doppia con = 0. In tal caso, riscrivendo la (1.3.17) nella forma 10 X Ui ; MU Vi ; MV R = N1 S S U (1.3.18) V si vede che, per l'indipendenza stocastica di Ui Vi, X E fRg = N1 E Ui ;S MU E Vi ;S MV : (1.3.19) U V Ora, essendo le funzioni entro parentesi funzioni dispari, mentre le distribuzioni dei vettori U. i .. UN V.i .. VN sono funzioni pari, si ha chiaramente che ciascuna delle medie di (1.3.19) e nulla, cos E fRg = 0 : (1.3.20) Pertanto R e uno stimatore corretto di = 0 (ed X e normale) e possibile anche vedere che se 6= 0 R non e uno stimatore unbiased. Osservazione 1.3.3: (diminuzione del bias). Se T e uno stimatore corretto di , in generale g(T ) non e uno stimatore corretto di g(), almeno in modo esatto, a meno che g non sia una funzione lineare. Tuttavia se 2(T ) e abbastanza piccola si puo porre in modo approssimato E fg(T )g = g(E fT g) = g() : (1.3.21) Supposto che g sia una funzione liscia, ad esempio con derivata terza continua, (1.3.21) equivale a trascurare, come termine principale, la quantita 1=2 g () T2 , infatti 00 E fg(T )g = E fg() + (T ; )g () + 12 g ()(T ; )2 + o2 g : 0 11 00 Qualora l'approssimazione (1.3.21) non bastasse, si puo ridurre il bias di g(T ), portandolo ad un innitesimo di grado superiore, ponendo G = g(T ) ; 12 T2 g (T ) : 00 (1.3.22) E infatti facile vedere che, per una distribuzione simmetrica di T e per una g abbastanza liscia, (1.3.22) e O4, al contrario di g(T ) che e O2. Si osservi che (1.3.22) e utile quando si sappia calcolare T2 . 1.4 Consistenza Vogliamo ora caratterizzare il comportamento degli stimatori quando la numerosita del campione N tenda all'innito. E ovvio che sia desiderabile che la distribuzione di T (X (N ) ) tenda a concentrarsi sempre piu attorno a , cioe che T (X (N ) ) tenda in probabilita a . Diremo che uno stimatore per cui (N ) ) = in P lim T ( X N !1 (1.4.1) e consistente. In pratica spesso anziche provare la (1.4.1) si fa ricorso alla condizione su ciente della convergenza in media quadratica, ovvero lim E fT (X (N ) g = lim 2fT (X (N ) g = 0 : N N !1 !1 (1.4.2) (1.4.3) Notiamo che la (1.4.2) equivale a dire che il bias b() tende a zero per N ! 1. Esempio 1.4.1: la media campionaria M e uno stimatore consistente di . Infatti M e corretta, cos che b = 0, ed inoltre, usando la propagazione degli errori, 1 2 X 2 1 X 2 (1.4.4) (M) = N 2 X = N =0 N 12 cos che la (1.4.2) e (1.4.3) sono vericate. Esempio 1.4.2: se la v.c. X ha momento quarto limitato ( 4 < +1) allora S 2 e uno stimatore consistente di 2. In e!etti (1.4.2) e chiaramente vericata perche E fS 2 g = (N ; 1)=N ] 2. Quanto alla varianza di S 2 , possiamo osservare che, se 4 < +1, si ha E fS 4 = N N; 1 1 ; N2 + N32 1 4 = +0 N : g 4 + N1 1 ; N1 2 4 = (1.4.5) Percio e anche 2 fS g = EfS g ; 2 4 1; 1 N 2 4 =0 1 N e la consistenza e provata. P 2 Si puo notare che poiche S 2 ! si ha anche chiaramente 2 2 P 2 S2 = N=(N ; 1)S ! , cioe anche S 2 e uno stimatore consistente di . Inne osserviamo che nel caso la X sia una v.c. normale, dalla (1.4.5) e dalla nota relazione 4 = 3 4 si ricava 2 4 (S 2) = N2 ; 1 : (1.4.6) Esempio 1.4.3: sia data la famiglia di distribuzioni uniformi fX (x ) = 1= 0 x 0 x < 0 x > : si puo dimostrare che la seguente statistica campionaria Xi T (X ) = max i 13 (1.4.7) e uno stimatore consistente di . Infatti, per ricavare la distribuzione di T , basta notare che FT (t) = P fT tg = P fX1 t : : : XN tg = = P fX tg]N = FX (t)N cos che fT (t) = N FX (t)N 1 fX (t) ; (T = max Xi) : Usando la (1.4.7) per la fX e per la FX si ha per T = max Xi, 8 tN 1 > <N 0t fT (t) = > N : 0 <t ; cos che Z 1 E fT g = E fmax Xig = N N tN dt = N N+ 1 = + 0 N1 0 : Inoltre Z 1 E fT g = N N tN +1 dt = N N+ 2 2 = 2 + 0 N1 0 2 cos che 2(T ) ! 0 per N ! 1 e la consistenza di T e provata. Osservazione 1.4.1: contrariamente a quanto avviene per la proprieta di correttezza, la consistenza viene mantenuta se si passa da T a g(T ) per una vasta classe di funzioni. In e!etti se g() e continua in si avra che l'immagine inversa di un intorno qualsiasi di g() g() ; " g() + "] contiene un intorno I di . 14 Figura 1.4.1: Pertanto si ha P (jg(T ) ; g()j < ") > P (T 2 I") d'altro canto per la consistenza di T come stimatore di , lim P (T 2 I") = 1 N !1 8 cos che risulta anche lim P (jg(T ) ; g()j < ") = 1 N !1 8" (1.4.8) e g(T ) e uno stimatore consistente di g(). La stessa proprieta vale a piu dimensioni anche se il provare cio richiede un piu attento esame. Premettiamo infatti un lemma. Lemma 1.4.1: sia T = (T1 : : : Tp) uno stimatore di = = (1 : : : p) nel senso che per ogni componente Ti ! i 15 in P : allora T e uno stimatore consistente p-dimensionale di , cioe, preso un qualunque intorno p-dimensionale di I , risulta lim P fT(N ) 2 I g = 1 : N !1 (1.4.9) Infatti notiamo che data una successione P (N ) di distribuzioni di probabilita tali che lim P (N ) (A) = 1 N lim P (N )(B ) = 1 N !1 !1 dovra essere anche necessariamente lim P (N ) (A \ B ) = 1 : N !1 Infatti, ssato ", per N abbastanza grande P (N ) (A) > 1 ; ", e quindi P (N )(B ; A) < ". D'altra parte P (N )(B ) > 1 ; ", e quindi P (N ) (A \ B ) = = P (N ) (B ) ; P (N )(B ; A) > 1 ; 2". Tale proprieta si estende immediatamente al caso di un numero nito di intersezioni. Ne segue che la proprieta (1.4.9) e vera perche ogni intorno di conterra un cubetto abbastanza piccolo f"1g f"2g : : : f"pg = C e perche tale cubo puo essere visto anche come l'intersezione degli \strati" p-dimensionali C = f1 2 "1 82 : : : 8pg \ f2 2 "2 81 : : : 8pg \ : : : : Poiche d'altronde per ipotesi P fT1 2 "1 8T2 : : : 8Tpg ! 1 ed analogamente per le altre componenti, si vede che dovra pure essere 16 P fT 2 C g ! 1 cioe vale la (1.4.9). Una volta dimostrato il lemma, procedendo in modo identico al caso monodimensionale si vede che: se Ti sono stimatori consistenti di i e se g e una funzione continua in = (1 : : : p), allora g(T ) = g(T1 : : : Tp) e uno stimatore consistente di g(). Esempio 1.4.4: il coe ciente di correlazione lineare campionario X MY R = MXY S; M S (1.4.10) MXY = N1 X XiYi (1.4.11) X Y dove si e posto e uno stimatore consistente di . Applicando l'Osservazione 1.4.1 e sufciente supporre che X Y 6= 0, cos che (1.4.10) denisca una funzione continua di tutti i suoi argomenti, e provare che MXY ! XY in P , cosa che viene lasciata al lettore come esercizio. 1.5 Sucienza - Sucienza minimale - Completezza Larga parte della teoria della stima e dominata da un principio che potremmo chiamare principio di verosimiglianza: questo stabilisce che tutte le informazioni che si possono ottenere su un parametro a partire da un campione fx1 : : : xN g, debbono essere ricavate dalla forma della likelihood L(x ). Cio introduce il concetto di su cienza, per tener conto che spesso x entra in L(x ) solo in particolari combinazioni, le quali dunque contengono tutta l'informazione proveniente dal campione, su . Prendiamo un esempio semplice. Esempio 1.5.1: consideriamo la distribuzione 17 fX (x ) = 1 e ; x= x0 che genera la likelihood L(x ) = 1N e 1= ; P xi xi 0 : Poiche i valori campionari entrano in L solo tramite la statistica S= X P x , diciamo che i Xi contiene tutta l'informazione, per la famiglia data, concernente il parametro . Ne segue che, accettando tale principio, si e portati a cercare i possibili stimatori di fra le funzioni di S , in modo tale che campioni diversi che danno lo stesso valore di S , portino anche alle stesse stime di . Questo concetto e posto su basi piu precise con la denizione di statistica suciente: diciamo che S (X (N ) ) e una statistica su ciente per il parametro e la famiglia fX (x ), se la distribuzione della v. campionaria X condizionata a S = s risulta indipendente da . In termini analitici, ricordando la denizione di variabile condizionata, si puo vericare che 1 L(x jS = s) = Z L(x )jgrad S (x)j = h(x) 1 L(x )jgrad S (x)j d ; ; f S =s g 18 (1.5.1) dove x 2 fS = sg, e d e l'elemento d'area di fS (x) = sg. Notando che la distribuzione marginale di S , fS (s ) = Z L(x )j grad S (x))j 1d ; f S =s g (1.5.2) ovviamente non puo che essere funzione di s e di fS (s ) = K (s ) si trova che la condizione (1.5.1) puo essere anche scritta nella forma equivalente L(x ) = K (s ) h(x)jgrad S (x)j K (s )H (x) (s = S (x)) : (1.5.3) Si perviene cos al cosiddetto teorema di fattorizzazione che a!erma appunto che condizione necessaria e su ciente a nche S (X ) sia una statistica su ciente per e per la famiglia fX (x ), e che la likelihood si fattorizzi nel prodotto di due fattori, uno dipendente solo da s = S (x) e , e l'altro dipendente solo da x. Esempio 1.5.2: sia X (N ) un campione di variabili i.i.d. (indipendenti, identicamente distribuite) tratte da una N 2 ]. La likelihood di X (N ) puo essere scritta come 1 L(x 2) = (2)N= e1=(22 )(x 2 N e)+ (x e) (1.5.4) N log 1=(22 )x+x 2e+ x+N2 ] (1.5.5) ; ; dove e = 1 1 : : : 1]+ : Riscrivendo la (1.5.4) come L = (21)N=2 e ; ; ; e confrontando col teorema di fattorizzazione si riconosce che: 19 se = ( ) allora S = (X X e X ) = (P Xi P Xi) se e nota e = allora S = (P Xi) se e nota e = allora S = P(Xi ; ) . + 2 2 + 2 2 2 Osservazione 1.5.1: si noti che l'esistenza di una statistica su ciente S (X ) per una famiglia fX (x ), non signica che S (X ) sia uno stimatore di , ma in generale S (X ) sara uno stimatore di una qualche funzione di . Cos nell'Esempio 1.5.2 quando = ( ) si ha E fS g = fN (2 + 2 ) Ng. Osservazione 1.5.2: il concetto di statistica su ciente si puo illustrare notando che si tratta di una funzione (o piu) le cui superci di livello S (x) = s hanno la particolarita che, condizionando (restringendo) la distribuzione L(x ) ad una di esse, si trova una distribuzione indipendente da . Ci si puo chiedere allora quale sia il sistema piu esteso di tali superci. Infatti la variabile S su ciente che sia costante sul sistema piu esteso di superci, sara anche quella che varia di meno, pur mantenendo tutta l'informazione su : si dira allora che S e suciente minimale. Per comprendere questa problematica si riprenda l'Esempio 1.5.1 e si osservi che chiaramente, per uguale alla media dell'esponenziale, la statistica bidimensionale NX1 S =( 0 ; i=1 Xi XN ) e su ciente tuttavia anche N X S=( i=1 Xi ) e su ciente e quest'ultima compendia, per cos dire, tutta l'informazione su che era contenuta in S . Per riconoscere se una statistica e su ciente minimale o no, si puo fare ricorso al seguente Lemma. 0 20 Lemma 1.5.1: se S (X ) e una statistica su ciente, su ogni supercie di livello Ss = fS (x) = sg il rapporto di likelihood L(x ) L(y ) (x y 2 Ss) e indipendente da . Infatti se S (x) = S (y) = s, per il teorema di fattorizzazione si ha L(x ) = K (s ) H (x) L(y ) = K (s ) H (y) e quindi L(x ) = H (x) : L(y ) H (y) (1.5.6) Ne segue che, se si vuol trovare la statistica su ciente minimale per una certa famiglia fX (x ), basta esaminare il rapporto di likelihood e vedere su quali superci esso diventa indipendente da . Infatti anche per la statistica minimale deve essere S = cost, sulle superci di tipo (1.5.6). E viceversa, se il sistema S = cost coincide con quello su cui si verica (1.5.6), allora S e minimale perche nessun'altra statistica puo essere costante su superci che non siano comprese nel sistema (1.5.6). Esempio 1.5.3: per un campione normale di v.c. i.i.d. con media e P P varianza incognite, S = ( Xi Xi2) e su ciente minimale. Infatti, scritta la likelihood come in (1.5.5) e formando il rapporto, si trova L(x ) = e(1=22 )(x+x L(y ) ; y+ y)+(=2 )(e+ x e+ y) ; e chiaro che tale rapporto risulta indipendente da e 21 : 2 solo se X X x+x = x2i = y+y = yi2 X X e+x = xi = e+y = yi cioe su quelle superci su cui X P x2i = cost X P xi = cost : Pertanto S = ( Xi2 Xi) e su ciente minimale. Esempio 1.5.4: vogliamo dimostrare, mediante il rapporto di likelihood, che per la distribuzione uniforme su (0 ) (cfr. Esempio 1.4.3) la statistica S (X ) = max Xi e su ciente minimale. Infatti notiamo che la funzione di likelihood puo essere scritta come L(x ) = (1=)N 0 t = max xi : 0 < t = max xi Quindi, formando il rapporto di likelihood e convenendo che questo rimanga costante se L(x ) = 0 L(y ) = 0 entrambe, si ha 8 1 0 tx ty > < L(x ) = 0 < tx 0 ty : L(y ) > > 1 < ty 0 tx : 1 < tx ty : Questo rapporto puo risultare costante solo se tx = max xi = ty = max yi il che dimostra appunto che S (X ) = max Xi e su ciente minimale. Chiudiamo questo paragrafo aggiungendo una denizione che ci sara utile in seguito: diremo che la statistica S (X ) e completa per la famiglia fX (x ) se la sola funzione h(S ) per cui vale 22 E fh(S )jg = 0 e 8 (1.5.7) h(S ) 0 (1.5.8) cioe se (1.5.7) implica (1.5.8). Esempio 1.5.5: supposto che 2 sia nota, S = (1=N ) P Xi e completa per N 2 ]. Infatti siccome S = N ( 2 =N )], si ha E fh(S )jg = = Z+ 1 pN p2 ;1 pN h(s) p e N s 2 ( ; e2 N=(22 ) Z +1 ; )2 =(22 ) ds = h(s)e Ns2 =(22 ) ; esN=2 ds ;1 cos che, per un teorema sulla trasformata di Laplace, E fh(S )jg = 0 implica h(s)e s2 N=22 ; 8 0 cioe h(s) 0 : 1.6 Stimatori di minima varianza: ecienza Nei paragra 3) e 4) abbiamo studiato gli stimatori in relazione al loro valore medio e al comportamento asintotico delle loro distribuzioni, cioe quando N ! 1. In questo paragrafo vogliamo studiare gli stimatori in base alla loro dispersione attorno a , poiche piu piccola sara questa, maggiore sara la probabilita che lo stimatore risulti vicino al valore da stimare. 23 In generale deniamo l'errore quadratico medio di stima (e.q.m.s.) come E 2 = E fT (X ) ; ]2 g = 2 fT (X ) g + E fT g ; ]2 = = 2 fT (X ) g + b2 () : (1.6.1) Si potrebbe pensare che minimizzando la (1.6.1) si possano direttamente trovare stimatori con proprieta ottimali: si puo tuttavia vedere che in generale non esiste una T (X ) che minimizzi la (1.6.1) per tutti i valori di , cos che per sfruttare tale approccio si e condotti a costruire criteri di ottimalita piu complessi. Resta tuttavia vero che minimizzare E 2, pero su una classe piu ristretta di possibili stimatori, permette in alcuni casi di ricavare facilmente risultati interessanti e in particolare fornisce esempi in cui si dimostra che uno stimatore biased puo talvolta essere meno disperso attorno a di tutti gli stimatori corretti dello stesso parametro. Esempio 1.6.1: supponiamo di voler stimare la varianza 2 di una popolazione normale, di cui sia incognita anche la media. Se si restringe la classe degli stimatori a quelli aventi le seguenti caratteristiche: siano invarianti per traslazione T (Xi + ) = T (Xi) siano quadratici omogenei nella v. campionaria T= X ik XiXk siano simmetrici rispetto a permutazioni degli indici f : : : N g 1 si puo vedere che questi hanno necessariamente la forma T (X ) = X (Xi ; M)2 = (N ; 1)S 2 : (1.6.2) Scegliendo = 1=(N ; 1), si ha l'unico stimatore corretto di tale classe, mentre, lasciando libero , si ha in generale un bias 24 b = (N ; 1) ; 1] 2 : (1.6.3) Ricordando che 4 (S 2) = 2 N ;1 si ha per E 2 in tale classe l'espressione 2 E 2 = 22(N ; 1) 4 + (N ; 1) ; 1]2 4 : (1.6.4) (1.6.5) Il minimo dell'e.q.m.s. e allora ottenuto per = N 1+ 1 cioe per X T (X ) = N 1+ 1 (Xi ; M)2 : (1.6.6) in corrispondenza il valore di E 2 e E 2 4 = N2 + 1 che come si vede e inferiore a (1.6.4). Volendoci limitare agli stimatori corretti, si ha b() 0 cos che l'e.q.m.s. di T (X ) coincide con la sua varianza. Ci si puo chiedere se esista tra gli stimatori corretti quello di minima varianza, poiche e chiaro che tale stimatore sia da preferirsi tra quelli della sua classe. L'esistenza dello stimatore unbiased di minima varianza puo essere provata come un teorema che, sebbene non di grande uso pratico, serve ad inquadrare bene l'argomento. 25 Teorema 1.6.1: sia data la likelihood L(x ) supponia- mo che esista una statistica S (x) su ciente minimale completa per e per L, e sia V uno stimatore unbiased qualsiasi di . Allora, T = E fV jS g (1.6.7) e uno stimatore di minima varianza tra gli stimatori corretti inoltre esso e unico. Notiamo, prima di dare la dimostrazione, che in conseguenza della denizione (1.6.7) come media condizionata, T risulta una funzione di S T = T (S ) (1.6.8) come e logico aspettarsi proprio perche S e una statistica su ciente minimale. Il teorema si dimostra a passi: a) T e uno stimatore corretto di . Infatti E fT g = E fT (S )g = E fT (S )g = E fE fV jS gg = E fV g = S S b) T e indipendente da V . cioe se V1 V2 sono due qualsiasi stimatori corretti di , e se T1 = E fV1 jS g T2 = E fV2jS g, allora T1 = T2. Infatti, posto h(S ) = T1(S ) ; T2(S ), si ha che E fh(S )g = ; = 0 8 e per la completezza di S , segue h(S ) 0. c) T e di minima varianza. Infatti, sia V un qualunque stimatore unbiased si ha 26 2 (V ) = E f(V = E f(V ; ) g = E f(V ; T + T ; ) g = (1.6.9) ; T ) g + (T ) + 2E f(V ; T )(T ; )g : 2 2 2 2 D'altro canto, usando la relazione E fg = E fE fjS gg S E f(V ; T )(T ; )g = E fT (S ) ; ]E fV jS g ; T (S )]g 0 S cos che dalla (1.6.9) deriveremo 2 (V ) = 2(T ) + E f(V ; T) g 2 (1.6.10) ovvero 2 (V ) 2(T ) (1.6.11) cioe T e di minima varianza. d) T e unico. Infatti la (1.6.11) puo valere col segno \=" solo se E f(V ; T ) g = 0 cioe se V = T con P = 1 : 2 Corollario 1.6.1: come conseguenza del Teorema 1.6.1 si vede subito anche che se uno stimatore T e funzione della statistica su ciente minimale ed e allo stesso tempo corretto, allora esso coincide con lo stimatore corretto di minima varianza. Infatti, se T = T (S ), e inoltre E fT g = , si ha anche T = E fT jS g = E fT (S )jS g = T (S ) : 27 Esempio 1.6.2: la varianza campionaria corretta corretto di 2 S 2 S e lo stimatore di minima varianza per una famiglia normale. Infatti e = 2 1 (X X )2 1 (X X 2 ) ; i i N ;1 N (N ; 1) 2 cioe SP e funzione P della statistica su ciente, per la famiglia normale, S = ( Xi Xi2) inoltre S 2 e corretto e quindi, in base al Corollario 1.6.1, esso e lo stimatore di minima varianza. Si noti che a questo punto resta provato anche che lo stimatore (1.6.6) ha e.q.m.s. inferiore ad ogni simatore corretto. Osservazione 1.6.1: si noti che se = g(') da luogo ad una nuova parametrizzazione della likelihood, ovvero se ' e un nuovo parametro che dipende in modo invertibile da , e se T = h(S ) e uno stimatore corretto di ', allora esso e anche di minima varianza. Infatti, se L(x ) e la likelihood in funzione di , posto L (x ') = L(x g(')) 0 e facile vedere tramite il teorema di fattorizzazione che se S e su ciente per , e anche su ciente per '. Pertanto l'a!ermazione consegue dal Corollario 1.6.1. Anche quando esiste un minimo della varianza tra gli stimatori corretti di un parametro , in genere non e facile calcolare tale minimo, e non e quindi semplice giudicare, per uno stimatore dato, quanto esso sia buono. Si preferisce allora ricorrere ad una maggiorazione, detta limite di Cramer-Rao, che specica un limite inferiore per la varianza di uno stimatore, biased o unbiased che sia, limite che puo essere il minimo od un numero inferiore ad esso. Teorema 1.6.2: (limite di Cramer-Rao). Data una v. campionaria di varabili i.i.d., con una densita di probabilita fX (x ) regolare in , e denita la variabile, dipendente da , 28 @ log L(X ) = X U = X @ log f (X ) U = @ X i i i i @ (1.6.12) per un qualsiasi stimatore T , con bias b(), vale la disuguaglianza 2 fT g 1E+fbU(g)] = ; 1E+f@b(U)]g 2 0 0 2 2 (1.6.13) inoltre (1.6.13) vale col segno di uguaglianza se e solo se esiste una relazione lineare tra T e U . Per dimostrare la (1.6.13) partiamo dall'identita Z RN L(x )dN x = 1 (1.6.14) e supponiamo che L = i fX (xi ) sia cos regolare da poter derivare due volte sotto l'integrale (1.6.14). Alla prima derivazione otteniamo, ricordando la denizione (1.6.12), Z = = = Z RN @ L(x )dN x = X N ZR i @ fX (xi ) fX (xk )dN x = X @ fX (xi ) k=i 6 fX (xi ) L(x )dN x = E fU g = 0 RN i (1.6.15) cioe U e una variabile a media nulla per ogni . Derivando una seconda volta e ragionando come in (1.6.15), si ha Z Z f@U gL(x )dN x + RN = E f@ U g + E fU g = 0 RN 2 29 U@ (L(x ))dN x = (1.6.16) ovvero, ricordando anche (1.6.15), E fU 2 g = Ora si riparte dall'identita + b() = 2 fU g = ;Ef@ U g : Z RN T (x)L(x )dN x : (1.6.17) (1.6.18) supposto di poter derivare sotto integrale, si trova 1 + b () = Z 0 RN T (x)@ L(x )dN x = E fT (X )U g : (1.6.19) Ricordando che E fU g = 0, la (1.6.19) si puo anche scrivere 1 + b () = (T U ) : 0 inne, rammentando che il coe ciente di correlazione tra due variabili soddisfa la disuguaglianza 2 1, e 2 = 1 se e solo se esiste una relazione lineare tra le due variabili, si puo scrivere 1 + b ()]2 2 (T ) 2(U ) 0 (1.6.20) che, tramite la (1.6.17) ci da il risultato richiesto. Si noti che per stimatori corretti la (1.6.20) diventa semplicemente 2 (T ) ; 1 : E f@ U g La funzione I () = ;E f@ U g = ;E f@2 log L(X )g = 2 (U ) (1.6.21) e detta informazione in quanto e una misura del massimo di informazione ottenibile dai dati empirici sul parametro . 30 Osservazione 1.6.2: il limite inferiore di Cramer-Rao e signicativo, nel senso che esistono famiglie di distribuzioni per le quali esso e raggiungibile, cioe e un vero minimo. Infatti la (1.6.20) puo valere col segno di uguaglianza se e solo se U = @ log L(x ) = AT (x) + C : (1.6.22) in tale relazione A e C possono essere funzioni di senza con cio alterare il fatto che 2(U T ) = 2(U ) 2 (T ). Integrando la (1.6.22) rispetto a si ha log L(x ) = aT (x) + c + d dove a= Z A d c = Z (1.6.23) C d e la costante di integrazione d puo essere funzione di x in quanto @ d(x) 0. Riscritta, la (1.6.23) ci da (1.6.24) L(x ) expfa()T (x) + c() + d(x)g una qualsiasi likelihood della forma (1.6.24) e detta appartenere alla famiglia esponenziale. Dunque, il limite minimo di varianza puo essere raggiunto solo per funzioni di densita della famiglia esponenziale, che tuttavia come e facile provare comprende molte importanti variabili, quali ad esempio la normale, la ;, la binomiale e la poissoniana. Osservazione 1.6.3: per uno stimatore corretto che non raggiunga il limite di Cramer-Rao si potra introdurre un indice, che misuri l'ecienza dello stimatore stesso, nei termini = 2(T )I ()] 1 : (1.6.25) 01 (1.6.26) ; E chiaro che 31 e l'e cienza puo essere pari ad 1 se e solo se fX (x ) e della famiglia esponenziale. Osservazione 1.6.4: la trattazione del limite di Cramer-Rao (1.6.13) e essenzialmente monodimensionale. Essa vale ovviamente anche nel caso in cui sia la componente di un vettore e T sia il corrispondente stimatore. Tuttavia nel caso multidimensionale e possibile ottenere un limite inferiore piu stretto, che ci proponiamo qui di dare. Come gia nel caso monodimensionale, osserviamo anche qui che si puo scrivere @k L(x ) = X @k fX (xi ) i fX (xi ) L(x ) : cos, introdotta l'operazione vettoriale @ F () = @1 F () : : : @pF ()]+ (1.6.27) U (x ) = @ log L(x ) (1.6.28) e denendo si puo scrivere (1.6.29) @ L(x ) = U L(x ) : Supporremo L abbastanza regolare in per poter e!ettuare tutte le derivazioni sotto integrale che ci serviranno. Pertanto, partendo da Z e di!erenziando, troviamo E fU g = L(x )dN x = 1 Z U L(x )dN x = 0 : (1.6.30) Se ora applichiamo l'operatore @ a E fU g+ notando anche che @ @+ F e la matrice delle derivate seconde di F , si ha 32 Z Z @ U L(x )dN x + U U +L(x )dN x = 0 : + (1.6.31) Poiche vale la (1.6.30), si vede che E fU U +g = CUU inoltre deniamo come matrice di informazione (di Fisher) I () Z ; @ U L(x )dN x = Z ; f@ @ log LgLdN x = ;E f@ @ log Lg : = = = + + (1.6.32) + In questo modo la (1.6.31) diventa CUU = I () : (1.6.33) sia ora T (x) uno stimatore corretto, p-dimensionale, di : dall'identita = + Z T (x)+L(x )dN x applicando @ si trova I= Z U T +L(x )dN x : dove I e la matrice identita, da non confondersi con l'informazione scalare qui sempre indicata con I () (vedi (1.6.21)). Ricordando la (1.6.30), questa si puo scrivere CUT = E fU (T ; T )+g = I : 33 (1.6.34) Moltiplichiamo a destra e a sinistra per due vettori + e +CUT = E f+U (T ; T )+ g = + : (1.6.35) d'altro canto, come nel caso monodimensionale E f+U (T ; T )+ g2 2 f U g + 2 ( +T ) da cui, applicando la propagazione della covarianza (+ )2 = E f+U (T ; T )+ g2 (+CUU )( +CTT ) : (1.6.36) La (1.6.36) vale per ogni e ogni e se vale col segno di uguaglianza signica che esiste una relazione lineare tra +U e +T . Si noti che, scelto ad esempio + = 1 0 : : : 0] e = , si ritrova lo stesso limite monodimensionale 2 (T1) = +CTT 2(1U ) : 1 tuttavia dalla (1.6.36) e possibile ottenere un limite inferiore piu forte. Infatti, riscritta la (1.6.36) come + 2 +CTT (+C ) (1.6.37) UU si puo cercare per ogni dato l' che rende massimo il secondo membro. Si puo vedere che tale e dato da = CUU1 ; cos che si giunge a 2 ( +T ) = +CTT +CUU1 = +I () 1 ; 34 ; (1.6.38) ovvero CTT I () 1 : (1.6.39) ; Questa relazione dice sostanzialmente che la varianze di +T , che e uno stimatore corretto di +, ha come limite inferiore +I 1 () . In particolare scegliendo ; ::::::i:::::: i+ = 0 : : : 010 : : : 0] si trovano i limiti inferiori delle varianze delle singole componenti Ti. Si puo osservare che I () ha sulla diagonale principale ik () = 2 (Uk ), cos che, quando I () sia diagonale, la (1.6.38) riproduce i singoli limiti monodimensionali. 1.7 Stime di massima verosimiglianza Fino ad ora abbiamo studiato le stime dei parametri in base alle loro proprieta e solo in un caso abbiamo denito un criterio generale per la ricerca di tali stime: nel caso degli stimatori corretti di minima varianza. Come si e visto in quel caso, pero, la soluzione generale del problema e ottenibile quando si conosca gia uno stimatore corretto, cosa non sempre facile. Vedremo in questo paragrafo invece un diverso principio di scelta degli stimatori detto di massima verosimiglianza (maximum likelihood). Questo criterio, inventato da Fisher, stabilisce di scegliere come stimatore di la statistica T (X ) tale che L(x ) : L(x T (x)) = max (1.7.1) Essenzialmente tale principio stabilisce di scegliere per un campione x dato, quel valore di per cui massima era la probabilita di estrarre proprio quell'x. Poiche il logaritmo e una funzione monotona, spesso anziche il massimo di L(x ) si preferisce cercare il massimo di log L(x ). Condizione 35 necessaria2 di massimo e l'annullarsi della derivata (o del gradiente) di log L(x ) cos che (1.7.1) viene tradotto nell'equazione U (x T (x)) = @ log L(x )j = T (x) 0 (1.7.2) da risolversi rispetto a T (x). Esempio 1.7.1: per i campioni normali con = ( 2) si trovano le stime ^ = (^ ^ 2 ) di m.l. delle equazioni X log L(x ) = ; N2 log 2 ; 21 2 (xi ; )2 + cost (1.7.3) 1 P ( x ; ) i 2 (1.7.4) U = ; N + 1 P(x ; )2 i 2 2 24 U (X ^ ^ 2) = 0 : La soluzione ^ di (1.7.5) e (1.7.5) P ^ = N1 PXi (1.7.6) ^ 2 = N1 (Xi ; ^)2 : Si puo notare che la stima di m.l. non e in generale unbiased. Esempio 1.7.2: siano Xi variabili normali con la stessa varianza e medie i funzioni lineari di un numero inferiore di parametri X = N 2 ] (1.7.7) = A (1.7.8) .1 = .. N : : : a1p 1 ... = ... a.11 A = .. aN 1 : : : aNp p < n p Oltre ad imporre la condizione necessaria, occorrerebbe anche vericare che la soluzione ottenuta da (1.7.2) corrisponda ad un vero massimo. Per alcune distribuzioni, come la normale, cio e automaticamente vericato. 2 36 supponiamo inoltre che A sia rango pieno3, cos che A+A sia invertibile. Posto + = (1 : : : p 2) otteniamo la stima m.l. dalle equazioni log L = ; N2 log ; 21 (x ; A) (x ; A) + cost 2 + 2 @ log L @ U = @ log L = @1 + = AN (x ;1 A) 2 ; 2 2 2 + 24 (x ; A)+(x ; A) (1.7.9) (1.7.10) U (x ^ ^ 2) = 0 : (1.7.11) ^ = (A+A) 1 A+X ^ 2 = 1 (X ; A)+(X ; A^) : N (1.7.12) (1.7.13) Il risultato e ; Notiamo che dei due stimatori, il primo e corretto E f^g = (A+A) 1 A+E fX g = (A+ A) 1A+A = ; ; (1.7.14) mentre il secondo non e corretto. Infatti notiamo che E fX ; A^g = 0 per la (1.7.14), cos che si puo scrivere, posto U = X ; A^. E f(X ; A^)+(X ; A^)g = TrCUU : Poiche I(N ) = identita a N dimensioni]. In eetti A di rango pieno equivale a dire che A = 0 ! = 0. D'altro canto se (A A) = 0 allora e anche 0 = + A+ A = (A)+ A = jAj2 ) A = 0 ! = 0: poiche A+ A e quadrata e non ammette soluzioni del sistema omogeneo associato, A+ A e invertibile. 3 + 37 U = I(N ) ; A(A+A) 1A+ ]X ; si ha CUU = I ; A(A+A) 1A+ ] 2I I ; A(A+A) 1A+ ] = = 2 I ; A(A+ A) 1A+] : (1.7.15) ; ; ; Ora notiamo che TrI(N ) = N mentre TrA(A+A) 1A+ = Tr(A+A) 1 A+A = TrI(p) = p : ; ; Riassumendo, dalla (1.7.15) si trova TrCUU = 2(N ; p) : pertanto tornando alla (1.7.13) si vede che ^ ^ E f ^ 2g = E f(X ; AN)(X ; A)g = (N N; p) 2 : Per fortuna e facile correggere tale bias ottenendo come stimatore corretto ^+ ^ ^o2 = (X ; A) (X ; A) : (1.7.16) N ;p Le formule (1.7.12) e (1.7.16) forniscono la soluzione di un problema detto \di modello lineare" nell'ambito delle variabili normali. L'argomento sara ripreso piu in generale trattando del metodo dei minimi quadrati. 38 Gli stimatori di m.l. non hanno in generale nessuna proprieta ottimale, a parte il signicato intuitivo legato alla loro denizione. Tuttavia essi restano importanti perche godono di diverse proprieta ottimali almeno asintoticamente, fatte salve alcune condizioni di regolarita (di!erenziabilita) della funzione di likelihood. In e!etti le equazioni per le stime di m.l. si scrivono come U (X ^) = 0 U = (@=@) log L] (1.7.17) ora notiamo che supposto che lo stimatore ^ = T (X ) ottenuto dalla (1.7.17) abbia un e.q.m.s. piccolo (rispetto ad una zona di variabilita signicativa di U in funzione di ), si puo sostituire l'equazione esatta (1.7.17), con quella linearizzata in corrispondenza al valore esatto di , X ) (^ ; ) 0 : U (X ) + @U (@ = (1.7.18) Dalla (1.7.18), che vale solo in modo approssimato, si possono ricavare interessanti proprieta, approssimate per l'appunto, dello stimatore . Infatti supponiamo ancora che valga q 2 f(@Uj )=(@k)g << E f(@Uj )=(@k )g = Ijk () 8j k nel caso monodimensionale, o analoghe condizioni per il caso multidimensionale in tal caso la (1.7.18) puo essere ulteriormente approssimata con ;I ()(^ ; ) + U (X ) = 0 (1.7.19) ^ = + I () 1 U (X ) : (1.7.20) che ha soluzione ; Ricordando che E fU g = 0 e CUU = I (), dalla (1.7.20) ricaviamo 39 E f^g = C^^ = I () 1 CUU I () 1 = I () 1 : ; ; ; (1.7.21) (1.7.22) Dunque, nell'ambito delle approssimazioni fatte, ^ e quasi unbiased e la sua matrice di covarianza e circa I () 1 : si osservi ancora che piu grande sara I () piu piccola sara la matrice di covarianza di ^ e quindi migliore (piu informativa) la stima. In particolare notiamo che ; U () = @ log L(x ) = X @ log fXj (xj ) = X Uj () (1.7.23) dove le Uj () sono variabili i.i.d. (possibilmente multidimensionali): supposto che per le Uj si possa applicare il teorema centrale della statistica (ad esempio Uj abbia momento 3 nito) si avra che U e asintoticamente normale, e anzi per le (1.7.21), (1.7.22) ^ (1.7.24) = N I () 1 ] : E interessante notare che, se per ogni variabile Uj si pone E fUj Uj+g = = C0 = Io(), si ha ; I () = CUU = N Io() (1.7.25) cos che la (1.7.24) puo anche essere scritta come I () ^ = N o N 1 ; : (1.7.26) Tale relazione mostra che ci si puo aspettare che ^ sia anche consistente, sebbene la dimostrazione accennata sia puramente euristica: in e!etti gode di tale proprieta sotto opportune ipotesi di regolarita per fX (x ). Il ragionamento fatto qui in modo puramente approssimato puo essere reso rigoroso prendendo il lim per N ! 1. Osservazione 1.7.1: si vuole dimostrare che per lo stimatore di m.l. monodimensionale ^, soluzione dell'equazione U (^) = 0, vale la proprieta asintotica, per N ! 1, 40 i( ) ^ = N N con ; 1 i() = E f@ log fX (X )]2g = 2(Uj ) 8j almeno sotto l'ipotesi restrittiva che le variabili i.i.d. Uj = @ log fXj (Xj ) @ Uj = Uj soddisno le condizioni del teorema centrale della statistica e le @2 Uj = Uj siano limitate con P = 1. In e!etti si puo notare che l'equazione di ^ puo essere scritta 0 00 0 = U (^) = U () + U ()(^ ; ) + 21 U ( )(^ ; )2 0 con j 00 ; j j^ ; j. Dividendo per 1 N p P e posto AN = Up() = pUj () = N 0 i()] N PN Uj () cost U ( ) BN = N = N = N ;i() N P U () cost cost 1 U ( ) 1 j CN = 2 N 3=2 = 2 N 3=2 = N p N 2 N p N = N (^ ; ) si ha che N risolve l'equazione 0 0 00 00 0 = AN + BN N + CN N2 cos che per N zione ! 1 N tende in probabilita alla soluzione dell'equap 0 = i()Z ; i() dove Z = Nlim i() ; !1 1=2 AN e una normale standardizzata. 41 2 Il metodo dei minimi quadrati 2.1 Introduzione L'idea di basare la teoria della stima sul principio di massima verosimiglianza, benche attraente, presenta alcuni seri incovenienti: in primo luogo bisogna supporre che sia nota la famiglia di distri buzioni da cui si estrae il campione inoltre il problema di trovare lo stimatore applicando le equazioni del massimo, e risolubile in forma analitica solo in pochi casi particolari inne, passando a soluzioni numeriche (cioe accontentandoci del valore numerico dello stimatore in corrispondenza ai valori numerici del campione), si vede che, a parte casi particolari o con un numero piccolo di incognite, il problema e insolubile. Si pone quindi il problema di ottenere per certe grandezze osservabili yifi = 1 : : : ng, cioe grandezze sulle quali si possono ricavare informazioni empiriche (osservazioni), stime che non implicano necessariamente la conoscenza di tutta la distribuzione ma solo certe sue caratteristiche, ad esempio medie e varianze-covarianze. Esempio 2.1.1: si consideri una maglia di una rete elettrica e si supponga di aver eseguito misure di di!erenze di potenziale ai nodi. In particolare (g. 2.1.1) si supponga di aver osservato: y1 = V12 = V2 ; V1 y2 = V23 = V3 ; V2 y3 = V31 = V3 ; V1 (2.1.1) Le osservazioni sono e!ettuate commettendo errori di misura e quindi vanno 42 descritte da variabili casuali Y1 Y2 Y3 di cui le variabili \osservabili" y1 y2 y3 denite in (2.1.1) sono i valori medi: E fYig = yi : (2.1.2) In virtu della struttura sica del sistema sottoposto ad osservazione, le medie yi non possono assumere tre valori arbitrari, ma a causa della loro forma (2.1.1) (ovvero per la denizione stessa di potenziale) deve essere y1 + y2 + y3 = 0 (2.1.3) viceversa si potrebbe anche notare che dati tre numeri qualsiasi che soddisfano la (2.1.3) e possibile pensare ad una maglia ideale che ha quei numeri come di!erenza di potenziale: in altri termini la (2.1.3) e la condizione necessaria e su ciente perche i tre numeri y1 y2 y3 descrivano le di!erenze di potenziale di una maglia di vertici 1,2,3. Naturalmente la presenza di una condizione come la (2.1.3) da una qualche informazione sugli errori di misura connessi che non saranno conosciuti individualmente, ma in una loro combinazione lineare. Infatti, se indichiamo con "i gli errori, Yi = yi + "i E f"ig = 0 (2.1.4) deve essere X Yi = X yi + X "i = X "i (2.1.5) in conseguenza della (2.1.3). Inoltre, conoscendo i procedimenti di misura, si possono fare ragionevoli ipotesi sulla matrice di covarianza di Y . Ad esempio, adottando particolari precauzioni da studiare caso per caso, si potra supporre che le misure siano indipendenti e cos le componenti di Y + = Y1 Y2 Y3] saranno tra loro incorrelate: inoltre supponendo che le osservazioni abbiano tutte la stessa precisione si potra supporre direttamente CY Y = o2I 43 (2.1.6) ove o2 ha il senso della varianza comune a tutte le osservazioni, e rimane una incognita del nostro problema. Cos il problema relativo a questo esempio puo essere sintetizzato come: conoscendo un vettore Yo1 Yo = Yo2 Yo3 di valori osservati, estratti da una v.c. a 3 dimensioni Y1 Y = Y2 Y3 di cui e denito il vettore dei valori medi, cioe delle osservabili, y1 E fY g = y = y2 y3 si sappia che tale media deve soddisfare la condizione (lineare) X yi = 0 (2.1.7) e sapendo che la covarianza di Y ha la forma CY Y = o2I (2.1.8) si vogliono trovare stime di y e di o2 , nel senso che che si cercano stimatori y^ ^o2 rispettivamente di y e o2, che siano funzioni della variabile campionaria Y , dunque calcolabili conoscendo il vettore delle osservazioni Yo. 44 Inoltre si vorrebbe conoscere anche la dispersione delle nostre stime ovvero Cy^y^ e, dove possibile, anche 2 (^o2). Non manchiamo di osservare, chiudendo questo esempio, che per motivi di simmetria (cioe non avendo motivi di pensare che un errore sia piu grande degli altri) una soluzione intuitiva del nostro probema di stima puo essere data ponendo (per la (2.1.5)) 3 3 X X 1 1 "^1 = "^2 = "^3 = 3 "k = 3 Yok : k=1 k=1 (2.1.9) inoltre si ha anche che, posto "= X Yok = X "k (E f"g = 0) (2.1.10) risulta per la legge di propagazione degli errori E f"2g = 2 (") = 3 o2 cos che ^o2 = 1 "2 3 (2.1.11) e uno stimatore corretto di o2. In base alla (2.1.9), la (2.1.11) puo essere scritta come ^o2 = 31 "2 = 3^"2i (2.1.12) Osservazione 2.1.1: si puo notare che, pur nella sua elementarita, l'E- sempio 2.1.1 e stato svolto prescindendo completamente dalla distribuzione della variabile 3-D, Y . Una maniera per geometrizzare la scelta (2.1.9) e la seguente: si consideri lo spazio R3 in cui ha sede la distribuzione di Y , e in esso il piano V di equazione (2.1.7). 45 Geometria elementare dei minimi quadrati: V= Yo = y= y^ = varieta dei valori ammissibili vettore delle osservazioni vettore delle osservabili vettore delle stime Figura 2.1.2 La distribuzione di Y ha un valore medio y 2 V (cfr. g. 2.1.2), detta varieta dei valori ammissibili, il vettore delle osservazioni Yo in generale sara al di fuori di V si cerca un vettore y^ 2 V che sia \il piu vicino possibile" a y, e sia una funzione lineare di Yo. Se Yo ha una distribuzione isotropa, cioe senza una tendenza a muoversi dalla media y piu in una direzione che in un'altra, cosa che e indicata dalla forma di covarianza (2.1.8), e \intuitivo" scegliere per y^ la proiezione ortogonale di Yo su V , cioe quel vettore y^ = Yo ; "^ (2.1.13) P che tra tutti quelli di V ( y^i = 0) rende minima la distanza "^+"^ = X "^2i = X (Yoi ; y^i)2 = (Yo ; y^)+(Yo ; y^) = min : (2.1.14) In questo caso il quadrato della distanza minimizzata (2.1.13) serve anche a stimare o2 in base alla (2.1.12). Osservazione 2.1.2: il criterio di minimizzare la somma dei quadrati (2.1.14) coincide perfettamente con quello di maximum likelihood quando si supponga che Y N y Infatti in tal caso 46 oI] 2 : (2.1.15) L(Y y) = 1 p ( 2)n ( ne 1=2o2 ) ; o P(Yi yi ) 2 ; (2.1.16) cos che P la condizione di massimo di L rispetto a y, soggetto alla condizione yi = 0, si traduce nella condizione di minimo (2.1.14). In questo senso si puo pensare di generalizzare la (2.1.14) al caso in cui CY Y = o2Q (2.1.17) (Yo ; y^)+Q 1 (Yo ; y^) = min : (2.1.18) con Q 6= I , prendendo, in analogia al caso della distribuzione normale, ; Questo stesso principio sara giusticato anche in base ad una richiesta di invarianza del procedimento di stima per trasformazioni lineari. 2.2 Formulazione generale del problema (caso lineare) Analizzando l'esempio del paragrafo precedente perveniamo alla seguente formulazione generale del problema di stima. E dato un vettore di valori osservati Y.o1 Yo = .. Yon Y.1 tratto da una v.c. n-dimensionale, Y = .. Yn (2.2.1) , di cui non si conosce in generale la distribuzione, ma di cui si sa che il valore medio y = E fY g 47 (2.2.2) detto anche vettore delle variabili osservabili, per motivi sici o geometrici e ristretto a stare su una varieta lineare a m dimensioni (m < n) y2V (2.2.3) detta anche varieta dei valori ammissibili questo costituisce il cos detto modello deterministico. Inoltre sulla base della conoscenza delle modalita di misura si puo a!ermare che la covarianza di Y ha la forma CY Y = o2Q (2.2.4) con o2 incognito e Q una matrice nota, strettamente denita positiva: questo costituisce il cos detto modello stocastico4. Sulla base del vettore Yo, del modello deterministico (2.2.3), del modello stocastico (2.2.4) si vogliono trovare le stime y^ y ^o2 2 o con matrice di covarianza Cy^y^ e, se possibile, 2 f ^o g : 2 Naturalmente per arrivare a denire le stime richieste occorre porre delle condizioni aggiuntive: queste costituiscono il principio dei minimi quadrati, che consiste nell'imporre che y^ sia tale da soddisfare (Yo ; y^)+Q 1 (Yo ; y^) = min y^ 2 V ; (2.2.5) e che ^o2 sia proporzionale alla forma quadratica minimizzata nella (2.2.5), Il modello stocastico (2.2.4) implica che siano note le varianze 2 fY g a meno di un fattore proporzionale e le correlazioni tra le componenti questa ipotesi copre molti casi metrologicamente interessanti. 4 i 48 ^o2 = c(Yo ; y^)+Q 1 (Yo ; y^) ; (2.2.6) con la costante c determinata in modo tale che ^o2 sia uno stimatore corretto di o2 . Osservazione 2.2.1: nel caso Q = I il principio (2.2.5) e interpretato come la ricerca del punto di minima distanza (euclidea), p d(Y y^) = (Y ; y^) (Y ; y^) + tra quelli che appartengono alla varieta V . Come si e visto anche nel paragrafo 2.1 la ricerca del punto piu vicino equivale a riconoscere che non esiste nessuna direzione preferenziale proprio perche Q = I . Vogliamo far vedere che, se Q = I , il principio (2.2.5) puo essere derivato sulla base della seguente richiesta: sia y^ lo stimatore lineare a cui arriveremo in base al criterio scelto sia ora T una trasformazione lineare qualunque che trasforma Y in una variabile isotropa, ad esempio T =Q 1=2 ; 5 (2.2.7) e sia U = TY (2.2.8) CUU = o2 I : (2.2.9) cos che Al vettore osservato Yo corrispondera Uo = TYo, mentre alla varieta V corrispondera la varieta degli u possibili del tipo u = Ty y 2 V . Sulla base di questo problema per u, sappiamo che una stima u^ puo essere ricavata in base al principio Ricordiamo che questa e per denizione l'unica matrice simmetrica, denita positiva per cui T 2 = Q;1. 5 49 (U ; u^)+(U ; u^) = min : u^ = Ty (y 2 V ) (2.2.10) Ora richiediamo che la stima y^ che stiamo cercando sia legata a u^ proprio dalla stessa relazione, cioe che se u^ e la stima di u, allora y^ = T 1u^ sia la stima ricercata di y. Ma in tal caso ispezionando la (2.2.10) si vede che y^ puo essere ricavata direttamente dal principio ; (Y ; y^)+Q 1 (Y y^ 2 V ; ; y^) = min che e appunto quanto volevamo dimostrare. Osservazione 2.2.2: poiche in analisi matematica il concetto di distanza e introdotto in termini generali tramite le sue proprieta formali e in particolare negli spazi euclidei Rn, ogni forma quadratica omogenea, denita positiva in senso stretto d(x y) = (x ; y)+P (x ; y) puo essere presa come nozione di distanza (subordinata alla metrica P = fPikg), si vede che il problema dei minimi quadrati (caso lineare) puo essere denito in termini formalmente molto semplici dicendo che: si vuole trovare il vettore y^ 2 V , varieta dei valori ammissibili, piu vicino ad ogni Y dato, secondo una metrica denita, tramite il modello stocastico (2.2.4), come P = Q 1 . ; Osservazione 2.2.3: dal punto di vista della simbologia notiamo che nelle formule dei paragra precedenti come in quelli che verranno si usera Y quando si vorra mettere in evidenza la variabile casuale che descrive il complesso delle osservazioni si usera Yo quando si vorra indicare lo specico vettore di osservazioni (numeri) trovati nell'eseguire le misure. Con y^ si indichera lo stimatore di y, il quale sara dunque a sua volta una variabile casuale funzione della variabile campionaria Y y^ = y^(Y ) tuttavia indicheremo con lo stesso simbolo anche il vettore delle stime, 50 cioe il vettore dei numeri che si ottengono sostituendo nello stimatore y^(Y ) il vettore delle osservazioni, y^ = y^(Yo). La funzione e ovviamente sempre la stessa, ma in un caso sottolineiamo il carattere di variabile casuale, nell'altro il carattere di vettore numerico. 2.3 Soluzione del problema di minimi quadrati: stimatori di osservabili e parametri Vogliamo risolvere il problema posto col principio dei minimi quadrati (2.2.5), (2.2.6) considerando sia il caso in cui la varieta lineare V sia data in forma parametrica y = Ax + a (2.3.1) (dim y = n dim x = m A = n m] rango A = m) che quello in cui essa e data in forma di equazioni di condizione By = b (2.3.2) (Yo ; y^)+Q 1 (Yo ; y^) = min (2.3.3) (B = c m] c = n ; m = No eq. di condizione, rango B = c BA = = 0 Ba = b). Si noti che nel caso normale con Q = I la soluzione e gia stata trovata nell'Esempio 1.7.2. Il metodo e perfettamente generalizzabile al caso Q 6= I , tuttavia si ritiene importante ricavare qui la soluzione direttamente dal principio variazionale ; sia per l'interesse del metodo in se che per la maggiore generalita dei risultati. Per rendere i conti il piu possibile sintetici conviene radunare le due forme (2.3.1), (2.3.2) in un'unica forma generale della varieta V , cioe Dy = Ax + d 51 (2.3.4) ed una volta ottenuto il risultato specializzare il modello (2.3.4) ai due casi di maggior interesse, fD = I d = ag ovvero fD = B A = 0 d = bg. Prima di procedere specichiamo le dimensioni di (2.3.4): dim y = n dim x = m m < n D = l n] A = l m] m < l n (2.3.5) inoltre supponiamo che il rango di D sia pieno, cioe che tutte le equazioni di condizione siano indipendenti, e che la matrice A sia o nulla (caso in cui mancano i parametri e si hanno pure equazioni di condizione (2.3.2)) oppure abbia rango pieno in modo che si possa sempre pensare di eliminare i parametri, riducendo le (2.3.4) ad un sistema di c = l ; n pure equazioni di condizione. Teorema 2.3.1: nelle condizioni poste i vettori (stimatori) x^ y^ che risolvono il problema variazionale (2.3.3), (2.3.4) (Yo ; y^)+Q 1 (Yo ; y^) = min Dy^ = Ax^ + d ; sono dati da x^ = N 1 A+K 1 (DYo ; d) y^ = Yo ; QD+K 1Uo ; ; ; (2.3.6) (2.3.7) dove 8 < K = DQD+ + 1 normale) : UNo == ADYKo ; AAx^ ; d (matrice (vettore degli scarti) ; (2.3.8) Il problema e risolto cercando i punti di stazionarieta di (Yo ; y^)Q 1 (Yo ; y^) sotto la condizione (2.3.4), mediante la tecnica dei moltiplicatori di Lagrange. Notiamo subito che la forma quadratica (2.3.3) e limitata inferiormente, perche denita positiva, mentre tende ; 52 all'innito quando Yo ; y^ ! 1: percio tale forma deve ammettere almeno un minimo assoluto e se le condizioni di stazionarieta hanno soluzione unica, allora tale soluzione deve necessariamente corrispondere al minimo. Per applicare la tecnica dei moltiplicatori di Lagrange, occorre sommare a (1=2)(Yo ; y^)+Q 1 (Yo ; y^) le componenti di Dy^ ; Ax^ ; d, ognuna moltiplicata per un suo moltiplicatore 1 2 : : : `: ovvero otteniamo la funzione obiettivo ; '(^x y^) = (1=2)(Yo ; y^)+Q 1(Yo ; y^) + (Dy^ ; Ax^ ; d)+ ; (2.3.9) dove + = 1 : : : `]. Per imporre le condizioni di stazionarieta a ', di!erenziamo tale funzione6 d' = ;dy^+Q 1(Yo ; y^) + dy^+D+ ; dx^+A+ : ; Annullando separatamente i coe cienti di dy^+ e dx^+ troviamo (Yo ; y^) + D+ = 0 A+ = 0 ;Q 1 ; (2.3.10) equazione cui va aggiunta la (2.3.4). Dalla prima delle (2.3.10) si ha y^ = Yo ; QD+ che, posta nella (2.3.4), da 6 Si osservi che se S e simmetrica e v un vettore qualunque, d(v+ Sv) = (dv+ )Sv + v+ Sdv = (dv+ )Sv + dv+ S + v = 2(dv+ )Sv : 53 (2.3.11) DYo ; K = Ax^ + d (K = DQD+) : (2.3.12) Notiamo che D Q e quindi anche K sono di rango pieno. Si ricava = K 1 (DYo ; d) ; K 1 Ax^ ; ; (2.3.13) e si usa le seconda delle (2.3.10), trovando N x^ = A+ K 1(DYo ; d) (N = A+K 1A) : ; ; (2.3.14) Ricordando che o A e identicamente nulla, nel qual caso la stima di x^ non si pone e si puo usare direttamente (2.3.13) in (2.3.11), oppure A e di rango pieno, la (2.3.14) puo essere risolta e da la (2.3.6): x^ = N 1A+ K 1(DYo ; d) : ; ; Sostituendo nella (2.3.13) si ha = K 1(DYo ; d ; Ax^) = K 1 Uo ; ; sostituendo nella (2.3.11) si trova la (2.3.7) y^ = Yo ; QD+K 1 Uo : ; Il vettore Uo = DYo ; d ; Ax^ prende il nome di vettore degli scarti delle equazioni. Osservazione 2.3.1: si consideri il caso in cui si abbiano pure equazioni parametriche y = Ax + a, dette anche equazioni di osservazione rispetto al modello generale si ha D=I d=a: 54 In questo caso la matrice A prende il nome di matrice disegno. Inoltre si ha K = Q N = A+ Q 1 A ; e risulta x^ = N 1A+ Q 1(Yo ; a) y^ = Yo ; QQ 1 Uo = Yo ; Uo = Ax^ + a : ; ; ; (2.3.15) Si puo notare che per la seconda delle (2.3.15) il vettore Uo coincide in questo caso anche con il vettore delle correzioni di Yo Yo ; y^. Osservazione 2.3.2: si prenda invece il caso in cui By = b, detto anche delle equazioni di condizione rispetto al modello generale si ha D=B A=0 d=b : dalle (2.3.14) e (2.3.14), posto A = 0, si ricava direttamente y^ = Yo ; QB + K 1(BYo ; b) (K = BQB + ) : ; (2.3.16) In questo caso il vettore Uo = BYo ; b prende il nome di vettore degli errori di chiusura (si veda l'Esempio 2.1.1). Osservazione 2.3.3: si noti che dal punto di vista della stima dei parametri x serve, piu che il vettore Yo, il vettore DYo ; d = Vo: questo signica che Vo contiene tutta l'informazione necessaria alla stima di x: in e!etti il modello generale Dy = Ax + d puo essere ritrasformato in modello puramente parametrico 55 v = Ax con \osservazione" Vo e con modello stocastico CV V = o2DQD+ = o2 K : 2.4 Covarianza degli stimatori e stima di 2 o Gli stimatori x^ y^ sono corretti: infatti dalla (2.3.6), osservando che per i valori medi x y deve valere la (2.3.4), troviamo E fx^g = N 1 A+K 1(DE fY g ; d) = N 1A+ K 1(Ax) = N 1 Nx = x : (2.4.1) Pertanto per il vettore U degli scarti delle equazioni si ha ; ; ; ; E fU g = E fDYo ; Ax^ ; dg = Dy ; Ax ; d = 0 ; (2.4.2) e quindi E fy^g = E fyg ; QD+K 1 E fU g = y : ; (2.4.3) Tra l'altro cio implica che il vettore delle correzioni "^ = Yo ; y^ e a media nulla E f"^g = 0 : (2.4.4) Cerchiamo ora le covarianze di x^ y^ e alcune matrici di covarianza e di cross-covarianza che ci serviranno in futuro. Covarianza di x^: ricordando che per ipotesi 56 CY Y = o2Q ed applicando la propagazione della covarianza, si ha Cx^x^ = o2 N 1 A+K 1DQD+K 1AN 1 = o2N 1 : ; ; ; ; (2.4.5) ; Covarianza di U: dalla denizione (2.3.8) e dalla (2.3.6) si vede che U = DY ; d ; AN 1 A+K 1(DY ; d) = = (I ; AN 1 A+K 1 )(DY ; d) : ; ; Poiche V = DY CUU = = ; ; d ha covarianza CV V = ; AN ; AN o (I 2 o (K 2 1 ; (2.4.6) ; o DQD 2 + = o2 K , si trova A+K 1 )K (I ; K 1 AN 1 A+) = 1 + A ): (2.4.7) ; ; ; ; Cross-covarianza di U e x^: notando che tanto U quanto x^ sono funzioni lineari dello stesso vettore V = DY ; d (cfr. (2.3.6) e (2.4.6)), si ha per denizione CU x^ = E fU (^x ; x)+ g = (I ; AN 1 A+ K 1)CV V K 1 AN 1 = = o2(I ; AN 1 A+K 1)AN 1 = (2.4.8) 2 1 + 1 1 = o (A ; AN A K A)N = 0 : ; ; ; ; ; ; ; ; ; ; Dunque U e x^ sono tra loro non correlati. Si osservi che di conseguenza anche il vettore delle correzioni "^ = Y ; y^, che in base alla (2.3.7) e funzione lineare di U , risultera non correlato ad x^ C"^x^ = 0 : (2.4.9) Covarianza di y^: allo scopo di calcolare questa matrice, conviene prima determinare a livello di servizio C"^"^ e C"y . Dalla denizione "^ = Y ; y^ = = QD+K 1 U e dalla (2.4.7) si trova ; 57 C"^"^ = o2QD+K 1 (K ; AN 1 A+)K 1DQ : ; ; ; (2.4.10) Inoltre per la (2.4.6) C"^Y = QD+K 1(I ; AN 1 A+K 1 )DCY Y = = o2QD+K 1 (K ; AN 1 A+)K 1DQ ; ; ; ; ; ; (2.4.11) cioe C"^Y = C"^"^. Essendo tale matrice simmetrica e anche C"^Y = C"^"^ = CY "^ : Per concludere, dalla relazione y^ = Y ; "^ si ha Cy^y^ = CY Y ; CY "^ ; C"^Y + C"^"^ = CY Y ; C"^"^ = = o2fQ ; QD+K 1 K ; AN 1 A+]K 1DQg : (2.4.12) ; ; ; Osservazione 2.4.1: nel caso di un modello di equazioni parametriche si ha D = I K = Q U = " N = A+Q 1A ; cos 8 < Cx^x^ = o2N 1 = C"^"^ = o2 Q ; AN 1 A+] : CCUU y^y^ = 2 AN 1 A+ : ; ; o ; Si noti come alla decomposizione 58 (2.4.13) Y ; y = Y ; y^ + y^ ; y = "^ + A(^x ; x) corrisponde in questo caso una perfettamente analoga decomposizione CY Y = o2 Q = C"^"^ + Cy^y^ = C"^"^ + ACx^x^ A+ a causa della non correlazione tra i due vettori (Q 1-ortogonali) "^ e y^;y = A(^x ; x). Osservazione 2.4.2: nel caso di un modello di pure equazioni di condizione risulta ; D = B A = 0 K = BQB + U = BY ; b vettore degli errori di chiusura in tal caso si ha 8 < CUU = o2K 2 + 1 C " ^"^ = o QB K BQ : Cy^y^ = 2fQ ; QB + K 1BQg : o ; ; (2.4.14) Anche in questo caso alla decomposizione Y ; y = "^ + y^ ; y corrisponde l'analoga relazione per le covarianze CY Y = C"^"^ + Cy^y^ : Stima di o2 : come gia annunciato nella formulazione generale del problema (cfr. x2.2), cercheremo uno stimatore di o2 nella forma ^o2 = c (Y ; y^) Q + ; 1 (Y ; y^) = c "^ Q + imponendo la condizione di correttezza 59 1 ; "^ (2.4.15) E f ^o2g = 2 o (2.4.16) allo scopo di determinare la costante c. Notando che7 "^+Q 1"^ = TrQ 1"^"^+ ; ; si vede che, (cfr. (2.4.10)): E f"^+Q 1 "^g = TrQ 1E f"^"^+g = TrQ 1C"^"^ = (2.4.17) 2 1 + 1 1 + 1 = o Tr Q QD K K ; AN A ]K DQg = = o2 fTrD+K 1DQ ; TrD+K 1 AN 1 A+K 1DQg : ; ; ; ; ; ; ; ; ; ; ; Per eseguire tale conto e su ciente ricordare che K = DQD+ = l l] (l = noequazioni (cfr:2:3:4)) N = A+K 1A = m m] (m = n parametri) : ; In e!etti TrD+K 1DQ = TrK 1DQD+ = TrK 1K = TrI (l) = l TrD+K 1AN 1 A+ K 1DQ = TrA+K 1DQD+K 1AN 1 = = TrA+K 1AN 1 = TrNN 1 = TrI (m) = m : ; ; ; ; ; ; ; ; ; ; ; ; Cio posto dalla (2.4.17) si deriva semplicemente 7 Ricordiamo una matrice Q = n n], la traccia e denita come P q : perche,taledata funzionale valgono le seguenti due proprieta TrQ = ii TrQ = Trq ] = Trq ] = TrQ+ XX XX Tr(A B ) = ( a b ) = ( b a ) = Tr(B A) : ik ki ik ki i k ki k 60 i ik E f"^+Q 1 "^g = ; 2 o fl ; mg : (2.4.18) Confrontando con (2.4.15), (2.4.16), si vede che uno stimatore corretto di o2 e dato da + 1 " ^ Q (2.4.19) ^o = l ; m"^ : Questo stimatore puo anche essere espresso esplicitamente in funzione degli scarti delle equazioni Uo , anziche delle correzioni "^: infatti, ricordando che ; 2 "^ = QD+K 1 U (Uo = DYo ; Ax^ ; d) ; si ha 2 o + K 1U = Ul ; m ; (2.4.20) Si puo osservare che quando si sia interessati solo alla stima dei parametri x, la (2.4.20) puo essere applicata senza cercare la stima y^, poiche gli scarti delle equazioni U dipendono da DYo ; d = Vo (U = Vo ; Ax^) e non da Yo direttamente. Osservazione 2.4.3: nel caso di un modello puramente parametrico U = "^ K = Q cos che (2.4.19) e (2.4.20) coincidono in modo evidente in piu si puo notare che in questo caso l = n equazioni = n = n osservabili: Osservazione 2.4.4: per il modello con pure equazioni di condizione, la forma (2.4.20) e certamente piu comoda della (2.4.19) in quanto il vettore Uo (vettore degli errori di chiusura) ha dimensioni inferiori ad "^. 61 2.5 Ottimalita degli stimatori m.q.: Teorema di Markov Poiche non abbiamo sin qui fatto alcuna ipotesi sulla distribuzione di Y , non e possibile porci in generale il problema dell'e cienza degli stimatori m.q. Tuttavia e possibile dimostrare una proprieta assai forte di tali stimatori, almeno se li si paragona ad altri stimatori lineari e corretti. In e!etti supponiamo di denire la varieta dei valori ammissibili tramite le equazioni parametriche y = Ax + a : (2.5.1) sappiamo che lo stimatore e poi lo stesso qualunque sia la forma con cui si descrive V . In tal caso gli stimatori m.q. sono dati da x^ = N 1 A+Q 1 (Yo ; a) = M^ (Yo ; a) y^ = AN 1 A+ Q 1(Yo ; a) + a = L^ (Yo ; a) + a : ; ; ; ; (2.5.2) (2.5.3) Piu in generale deniamo la classe degli stimatori lineari corretti tramite le formule x~ = M~ (Yo ; a) + m~ y~ = L~ (Yo ; a) + ~l (2.5.4) (2.5.5) e imponiamo la condizione che quando Yo 2 V , cioe Yo ; a = A per un qualche , allora risulti anche ~ + m~ x~ = MA y~ = Yo = A + a L~ A + ~l : (2.5.6) Le (2.5.6) vanno intese come identita rispetto a , il che comporta le condizioni ~ m~ = 0 A = LA~ ~l = a : I = MA 62 (2.5.7) Dunque la classe degli stimatori lineari e corretti e denita dalle (2.5.4), (2.5.5), con le condizioni aggiuntive (2.5.6), (2.5.7): e facile vedere che gli stimatori m.q. (2.5.2), (2.5.3) appartengono appunto a tale classe. ^ L^ e valida all'interno di tale classe: vale infatti il L'ottimalita di M seguente teorema: Teorema 2.5.1: (di Markov). Siano x~ y~ stimatori lineari corretti qualunque di x y e siano invece x^ y^ gli stimatori m.q. allora si ha Cx~x~ Cx^x^ Cy~y~ Cy^y^ : (2.5.8) Osservazione 2.5.1: le (2.5.8) signicano essenzialmente che Cx~x~ ; Cx^x^ e Cy~y~ ; Cy^y^ sono matrici denite positive, ovvero che +Cx~x~ +Cx^x^ +Cy~y~ +Cy^y^ 8 8 : (2.5.9) (2.5.10) A loro volta le (2.5.9), (2.5.10) possono essere cos interpretate: si voglia stimare una funzione lineare di x e una di y u = +x + uo v = +u + vo (2.5.11) se nelle (2.5.11) si usano gli stimatori m.q. per ottenere u^ = +x^+uo v^ = +y^ + vo si hanno stimatori lineari corretti di varianza minima 2 (^u) = +Cx^x^ 2 (^v) = +Cy^y^ tra tutti gli altri stimatori lineari corretti u~ = +x~ + uo v~ = +y~ + vo 2 (~u) 2(^u) 2(~v) 2(^v) : 63 Per dimostrare il Teorema 2.5.1 notiamo che x~ = M~ (Y x^ = M^ (Y ; a) ; a) con ~ = I MA ^ =I : MA (2.5.12) Si puo allora scrivere x~ = x~ ; x^ + x^ = + x^ Cx~x~ = C + Cx^ + Cx^ + Cx^x^ : (2.5.13) (2.5.14) e chiaro che se Cx^ = Cx^+ = 0 si ha Cx~x~ ; Cx^x^ = C 0 che coincide con la prima delle (2.5.8). D'altronde (2.5.15) Cx^ = (M~ ; M^ )CY Y M^ + = o2(M~ ; M^ )Q Q 1AN 1 = = o2 (M~ ; M^ )AN 1 = ~ ; MA ^ )N 1 = 0 = o2 (MA ; ; ; ; in conseguenza delle (2.5.12): la (2.5.15) e cos provata. Quanto alla 2a delle (2.5.8) essa discende immediatamente dalla osservazione che y~ = Ax~ + a y^ = Ax^ + a 64 cos che Cy~y~ ; Cy^y^ = ACx~x~ ; Cx^x^]A+ 0 : l'ultimo passaggio e facilmente vericato calcolando la forma quadratica associata +(Cy~y~ ; Cy^y^) : 2.6 Problemi di minimi quadrati con vincoli Volendo concentrare l'attenzione sulla stima dei parametri, assumiamo in questo paragrafo il modello parametrico y = Ax + a (2.6.1) cui ci si puo sempre ridurre. Ci sono almeno due casi importanti in cui puo interessare un modello deterministico come in (2.6.1) con l'aggiunta di (k) vincoli sui parametri x, cioe di equazioni (lineari) contenenti solo x: Hx = h (H = k m] k < m) (2.6.2) a) quando non si e sicuri del modello (2.6.1), cioe non si sa se si sono rappresentate tutte le variabili che descrivono un certo fenomeno, puo essere utile provare modelli diversi in cui si parte da uno molto generale che include molti parametri, ad altri piu semplici in cui una parte dei parametri e bloccata (ad esempio vincolata a zero) b) quando la matrice \disegno" A non sia di rango pieno: in questo caso infatti viene meno la corrispondenza biunivoca tra i punti della varieta V e quelli dello spazio dei parametri. Per eliminare questa indeterminazione e necessario introdurre vincoli sulla x. Illustriamo questi casi con due esempi. Esempio 2.6.1: un certo prodotto ottiene delle vendite settimanali crescenti che, a parte uttuazioni casuali, sembrano seguire un andamento lineare nel tempo Yi = x1 + x2 ti + "i (2.6.3) 65 a un certo tempo si introduce l'azione di una pubblicita, misurata ad esempio in termini di numero pi di spot televisivi per settimana: si ritiene che tale fattore possa inuenzare le vendite con una legge del tipo Yi = x1 + x2 ti + x3 pi + "i : (2.6.4) Per valutare l'e!ettiva e cacia della pubblicita, si possono compensare le \osservazioni" Yi (cioe trovare gli stimatori m.q. y^i) prima col modello (2.6.4) e poi col modello (2.6.3) per confrontare tra loro i risultati. Il modello (2.6.3) risulta essere contenuto nel modello (2.6.4), con l'imposizione del vincolo x3 = 0 : Esempio 2.6.2: in una rete di 7 punti si misurano dislivelli tra i punti stessi secondo uno schema come quello in gura: la congruenza geometrica della gura puo essere espressa per mezzo delle quote dei punti stessi, considerate come parametri. Le osservabili, corrispondenti ai lati della rete, sono 10: le equazioni di osservazione si scrivono Qki = Qk ; Qi (2.6.5) per i k opportuni. Le variabili Qki fanno parte del vettore delle osservabili, le quote Q1 : : : Q7 del vettore dei parametri. La matrice \disegno" A implicata dalle (2.6.5) e deciente di rango infatti esiste un vettore Figura 2.6.1 Q. 1 x = .. Q7 1. = .. 6= 0 1 tale per cui Ax = 0 : 66 Questa decienza e ovviamente legata alla mancata ssazione di una quota di riferimento essa puo essere rimossa imponendo un vincolo ai parametri. Ad esempio si potrebbe porre x1 = Q1 = 0 cioe prendere convenzionalmente il punto 1 come riferimento, oppure 7 X i=1 xi = 7 X i=1 Qi = 0 ovvero ssare convenzionalmente a zero la media delle quote. Il principio dei m.q. si scrive in questo caso come 8 < (1=2)(Yo ; y^)+Q 1 (Yo ; y^) = min Ax^ + a : y0^ == H x^ ; h ; (2.6.6) Sostituendo la II nella I e usando dei moltiplicatori di Lagrange per la III, si ha la funzione obiettivo ' = (1=2)(Yo ; Ax^ ; a)+Q 1(Yo ; Ax^ ; a) + (H x^ ; h)+ : ; Imponendo la condizione di stazionarieta rispetto ad x^ si trova ;A + Q 1 (Yo ; Ax^ ; a) + H + = 0 ; cui vanno aggiunti i vincoli, ovvero N x^ + H + = A+ Q 1(Yo ; a) (2.6.7) H x^ =h La (2.6.7) costituisce un sistema normale esteso che puo essere elaborato ulteriormente in modo analitico, solo se N e di rango pieno. ; 67 In questo caso si puo porre x^ = N 1 A+Q 1(Yo ; a) ; N 1 H + ; ; ; (2.6.8) che, con l'uso del vincolo, da HN 1A+Q 1 (Yo ; a) ; HN 1H + = h : ; ; ; (2.6.9) Risolvendo la (2.6.9) rispetto a e sostituendo in (2.6.8) si ottiene la soluzione cercata. Questo procedimento pero non e applicabile quando det N = 0 ed in generale non e conveniente. Piu comunemente (2.6.7) viene risolto con un metodo perturbativo che ha un'interessante interpretazione statistica. Il vincolo 0 = Hx ; h puo essere considerato come una osservazione nulla, Zo = 0, di un vettore casuale Z con media e covarianza nulla E fZ g = 0 : CZZ = 0 (2.6.10) con le condizioni (2.6.10) infatti Z assume il valore zero con probabilita 1. Ora, rilassando un poco tale condizione, possiamo prendere il vincolo come una osservazione Zo = 0 tratta da una variabile Z con media nulla 68 (2.6.11) E fZ g = 0 (2.6.12) CZZ = o2"I : (2.6.13) e con varianza innitesima Specichiamo anche che supponiamo che Z sia incorrelato con Y CZY = CY+Z = 0 : (2.6.14) In questo modo abbiamo un ampliamento del modello che include i vincoli delle osservabili: w = yz osservabili Wo = Y0o osservazioni Q 0 CWW = o 0 "I modello stocastico y ^ A w^ = z^ = H jx^j + ;ah modello deterministico: 2 La soluzione del relativo problema di m.q., fornita dalla (2.3.15) e allora 1 1 A Q 0 x^ = A H ] 0 (1=")I H 1 Y ; a 0 o + + Q A H ] 0 (1=")I h = + ; ; + ; = fA+Q 1 A + (1=")H +H g ; La stima di 2 o 1 ; fA Q + 1 ; (2.6.15) (Yo ; a) + (1=")H +hg : e data poi da + 1 + ^o2 = (Yo ; a ; Ax^) Q (Yo ; a ; Ax^) + 1="(h ; H x^) (h ; H x^) n+k;m (2.6.16) ; 69 in accordo con la (2.4.20), applicata nel presente caso. Osservazione 2.6.1: nel caso in cui il vincolo consista semplicemente nel porre a zero l'i-esimo parametro, si ha 1 1;i i 1+1 m H = 0 : : : : : : 0 1 0 : : : : : : 0] h = 0 : pertanto 1 ::: i ;1 i i+1 ::: m 0 0 0 + H H= 1 0 0 0 e tutta l'operazione (2.6.15) consiste essenzialmente nel modicare la matrice normale aggiungendo il numero \grande" 1=" sull'i-esimo elemento in diagonale. Esempio 2.6.3: si consideri una rete di livellazione (misura di dislivelli) sui tre punti P1 P2 P3 le quantita misurabili sono dunque Q12 y = Q23 Q31 Figura 2.6.2 e il modello deterministico relativo a questo esperimento puo essere espresso parametricamente in funzione del vettore x costituito dalle quote dei tre punti 70 Q1 x = Q2 Q3 : Le equazioni d'osservazione sono Q12 Y = Q23 Q31 Q2 ; Q1 = Ax = Q3 ; Q2 Q1 ; Q3 ;1 = 0 1 1 ;1 0 0 1 ;1 Q1 Q2 Q3 Supposto di avere un vettore di osservazioni Yo tratto da una Y che abbia come media y, il problema della stima di x non puo essere trattato in modo ordinario in quanto la matrice A non e di rango pieno in e!etti, 1 posto e = 1 , risulta 1 Ae = 0 cos che detA = 0, come e ovvio anche da un calcolo diretto. Cio signica che se x^ e una qualsiasi soluzione di minimi quadrati, altrettanto risulta x^ + e qualsiasi sia , in quanto gli scarti delle equazioni risultano invarianti per tale trasformazione V = Yo ; A(^x + e) = Yo ; Ax^ : Fisicamente cio deriva dal fatto che, non avendo ssato un'origine delle quote, le quote possono tutte essere variate di una stessa costante , senza mutare i dislivelli Qik . Proviamo ora a risolvere il corrispondente problema imponendo un vincolo che elimini l'ambiguita dell'origine ad esempio imponiamo 2 3 Q1 4 Q1 = 1 0 0] Q2 5 = Hx = 0 Q3 71 che ha l'ovvio signicato di ssare a 0 la quota di P1. Riscriviamo allora il sistema normale nella forma (2.6.7) tenendo conto che per ipotesi le misure sono indipendenti e con uguale precisione, ovvero 2 ;1 ;1 1 ;1 ;1 2 ;1 ;1 2 0 1 0 0 0 0 Q1 Q2 Q3 ;Q012 + Q031 = Q012 ; Q023 Q023 ; Q031 0 : Dalla quarta equazione si trova ovviamente Q^ 1 = 0, che sostituito nella seconda e terza, da Q^ 2 = Q^ 3 = Q012 ; 3 (" = Q + Q + Q ) 012 023 031 ;Q031 + 3 la prima equazione denisce il valore di ovvero =0 nella fattispecie. Vediamo ora come si sarebbe potuto risolvere, in modo approssimato, il problema applicando la formula (2.6.15) In questo caso la matrice dei pesi e l'unita cos che x^" = (A+A + 1" H +H ) 1A+Yo ; (2.6.17) tenuto conto che a = 0 h = 0. Se si usa il vincolo H = 1 0 0] si puo riscrivere il sistema normale di cui (2.6.17) e soluzione, nella forma 72 (D ;e e+)^x" = A+ Yo 3""+1 0 0 D = 0 3 0 : 0 0 3 (e = (1 1 1)+). Non e di cile vericare che la forma esplicita dell'inversa della normale e data da (D ; e e+) 1 = D 1 + 3(3" + 1)D 1e e+D 1 : ; ; ; ; Posto per brevita = 1 + 3", la forma esplicita della soluzione e allora 0 x^" = "(;Q012 + Q031 ) + 1+3 (Q012 ; Q023 ) + 3 (Q023 ; Q031 ) "(;Q012 + Q031 ) + 3 (Q012 ; Q023 ) + 1+3 (Q023 ; Q031 ) Da qui, preso " ! 0, ovvero ! 1, si vede che 0 lim x ^ = Q ; 012 " 0 " ;Q031 + ! 3 3 = x^ come volevasi dimostrare. 2.7 Problemi di stima non lineari Vogliamo vedere come si modica il metodo dei minimi quadrati quando il modello deterministico, che esprime il legame tra i valori medi delle osservabili, sia non lineare. Da un punto di vista astratto il principio non di!erisce, come enunciazione, da quello lineare: 73 dato un vettore di osservazioni Yo 2 Rn, estratto da una v.c. Y di covarianza CY Y = o2Q (2.7.1) e data una varieta V m-dimensionale, in Rn, cui si sa che deve appartenere la media di Y y = E fY g 2 V (2.7.2) si cerca lo stimatore y^ 2 V , di minima distanza da Yo secondo la matrice Q 1 , ; (Yo ; y^)+Q 1 (Yo ; y^) = min e una stima di o2 proporzionale a tale distanza ; (2.7.3) ^o2 = c (Yo ; y^)+Q 1 (Yo ; y^) (2.7.4) essendo ^o2 corretto, almeno in modo approssimato. ; Esempio 2.7.1: prima di iniziare la trattazione analitica puo essere utile considerare due esempi semplici, particolarmente evidenti dal punto di vista graco. Sia Yo = 012 il vettore delle osservazioni in un primo 1 e supponiamo 2 caso chehesso sia estratto da Y = N 1 (0 8) I , e in un secondo che i 1 Y = N 001 (0 8)2I : i modelli deterministici sono rispettivamente y2 = y12 ed y2 = 0:01y12. La situazione e rappresentata in g. 2.7.1 a) e b), dove sono anche mostrati gli stimatori m.q., nonche l'area circolare contenente circa il 99% delle probabilita. La prima osservazione e che nell'esempio a) lo stimatore m.q. y^ dista di piu da y che non Yo (jy^ ; yj > jYo ; yj), mentre nell'esempio b) le cose vanno meglio. La seconda osservazione e che vi sono punti nel piano in cui non e denita univocamente la proiezione sulla varieta V : ad esempio quelli dell'asse x = 0 8 , e con y > 0 5 per l'esempio a), con y > 50 per l'esempio b), ; Per comprendere questa aermazione basta osservare che se (t at2 ) e l'equazione parametrica di una parabola e se (o y) e un punto ssato dell'asse y, l'equazione d=dtf(y ; at2 )2 + t2 g = 0 ha solo la radice t = 0 quando y < (1=2)a. 8 74 a) b) Figura 2.7.1: ammettono due punti di minima distanza, simmetrici sui due rami della parabola. e anche utile pero osservare che per il caso a) questi punti sono ben dentro la zona P = 99%, mentre per il caso b) essi sono ben al di fuori. La terza osservazione e che 0,5 e 50 non sono altro che i raggi di curvatura delle due parabole nell'origine. Pensando un minuto alla denizione di raggio di curvatura, si comprende che il concetto e generale: in tutta l'area denita da un lato dalla curva, dall'altro dalla curva descritta dai suoi centri di curvatura, i punti ammettono una sola proiezione, a minima distanza, sulla curva stessa. L'esempio si generalizza alle varieta a m dimensioni in uno spazio Rn, prendendo il piu piccolo dei raggi di curvatura di V (m) . Si possono cos trarre alcune conclusioni: Conclusione 1: il metodo dei m.q. non puo, nel caso non-lineare, godere di proprieta di ottimalita come nel caso lineare. Conclusione 2: la soluzione del metodo dei m.q. con modelli non-lineari 75 puo non essere unica. Conclusione 3: l'unicita della soluzione e garantita se il raggio della zona di Y , misurato ad esempio da p 2 di interesse per la distribuzione n o quando Q = I 9 , e assai piu piccolo della curvatura della varieta V . In ultima analisi si puo a!ermare che il metodo dei m.q. funziona bene anche nel caso non lineare quando nella zona di interesse la varieta V non si discosta molto da una varieta lineare: qualora tale condizione fosse seriamente violata, come nell'esempio di g. 2.7.1 a), non ha senso applicare il metodo dei m.q. in quanto l'introduzione dell'informazione data dalla varieta V puo peggiorare quella contenuta nel vettore Yo. Dal punto di vista analitico i problemi non lineari vengono trattati per mezzo di una linearizzazione. Supponiamo ad esempio che la varieta V sia descritta nella forma piu generale, per mezzo di equazioni non lineari contenenti anche parametri aggiuntivi g(y x) = 0 g+(y x) = g1(y x) : : : gl(y x)] y 2 Rn x 2 Rm : (2.7.5) Cio signica quindi che si cercano y^ x^ tali che (Yo ; y^)Q 1(Yo ; y^) = min g(^y x^) = 0 ; (2.7.6) Supponiamo di conoscere dei \valori" approssimati y~ x~ e che nella zona di interesse g sia ben approssimabile con una funzione lineare: poniamo allora y^ = y~ + ^ x^ = x~ + ^ 9 Infatti quando Y e normale, e Q = I , si ha la relazione jY 76 (2.7.7) ; yj 2 = 2 . n e linearizziamo la seconda delle (2.7.6) come @ g~ @ g~ g(~y x~) + @y ^ + @x ^ = 0 (2.7.8) dove @ g~ @g (~y x~) i = @yk @ g@y ~ @g (~y x~) @y = i @xj : La (2.7.8) e il nostro nuovo modello, ora lineare, nelle variabili (parametri) ed (osservabili), per le quali si puo anche porre o = Yo ; y~ valori osservati 2 C = CY Y = o Q : (2.7.9) Notiamo che bastera porre ; @ g~ = D @y @ g~ @x = A g(~y x~) = gi(~y x~)] = d (2.7.10) per trasformare la (2.7.8) nel modello gia trattato nel paragrafo 2.2, D = A + d : cio insieme alle (2.7.9) permette di ricavare le stime cercate ^ ^ applicando la teoria lineare. Una volta ottenute le stime ^ ^ si possono calcolare le corrispondenti y^ x^ dalle (2.7.7). Osservazione 2.7.1: (iterazioni). ci si puo chiedere ora se il vettore y^ cos determinato appartenga veramente a V : infatti sono stati 77 determinati in modo che valesse solo l'equazione approssimata (2.7.8). A tale scopo si possono calcolare gli scarti ;v^ = g(^y x^) = g(~y + ^ x~ + ^) (2.7.11) e si puo valutare l'entita di jv^j decidendo se questa e trascurabile rispetto agli errori di misura.10 Qualora si pensi che v^ non sia abbastanza piccolo, cio signica che i valori approssimati y~ x~ erano troppo grossolani cos che l'approssimazione lineare (2.7.8) non risulta buona. In questo caso e necessario iterare il procedimento prendendo i valori y^1 x^1 determinati al primo passo come nuovi punti di linearizzazione. Si puo notare che nella nuova iterazione il nuovo vettore d risultera proprio uguale a d = g(^y1 x^1) in quanto y^1 x^1 sono i nuovi valori approssimati. Ripetendo piu volte lo stesso ragionamento si trova l'algoritmo iterativo Cio naturalmente riportato alle opportune unita di misura: ad esempio si potrebbe prendere come criterio 10 jv^j 2 << (l ; m)^2 = U +K ;1 U : o 78 Passo i+1 y^i+1 = y^i + ^i+1 x^i+1 = x^hi + ^i+1 i yi x^i ) Di = h; @g(^@y i yi x^i ) Ai = @g(^@x di = g(^yi x^i)] oi+1 = Yo ; y^i C = CY Y Di^i+1 = Ai^i+1 + di # Principio m.q. # ^i+1 ^i+1 Passo 1 y^o = y~ x^o = x~ (valori approssimati iniziali) Il criterio di stop e basato sulla dimensione di jdij quando jdij2 e abbastanza piccolo rispetto (l ; m) o2 , l'iterazione e fermata. Altri criteri di stop possono essere ssati in base a jij e jij. Osservazione 2.7.2: per la determinazione dei valori approssimati non esistono regole generali. Quando non vi sia di meglio, come primi valori approssimati di y si puo sempre scegliere y~ = Yo : (2.7.12) e bene notare che la (2.7.12) va presa come eguaglianza tra due vettori costanti non tra due v.c., senza di che sarebbe falsa la relazione C = CY Y : Per determinare i valori x~ molte volte si cerca di risolvere una parte delle equazioni (2.7.5) rispetto ad x ottenendo x = f (y) 79 e ponendo di conseguenza x~ = f (~y) (2.7.13) Osservazione 2.7.3: nel caso che g(y x) sia gia lineare in una delle due variabili, non serve linearizzare rispetto a quella variabile. Il caso piu importante e quello in cui le equazioni non lineari della varieta V sono in forma puramente parametrica (equazioni d'osservazione) y = g(x) : (2.7.14) Per linearizzare le (2.7.14) basta cercare dei valori approssimati x~, ad esempio invertendo una parte di tali equazioni e usando i corrispondenti valori osservati Yoi come termini noti. Successivamente si pone @ g~ y^ = g(~x + ~) = g(~x) + @x ^ = a + A^ : Trovata la soluzione ^ si puo calcolare y^1 = g(~x + ^1) (2.7.15) (2.7.16) in modo tale che la stima y^1 corrispondente appartiene esattamente a V , al contrario del caso generale (cfr. (2.7.11)). Lo schema puo essere iterato ponendo @g y^i+1 = g(^xi + ^i+1) = g(^xi) + @x (^xi) ^i+1 = ai + Ai^i+1 : (2.7.17) spesso, per semplicare i conti risparmiando tempo di calcolo, la matrice A non viene ricalcolata ogni volta, ma piuttosto viene lasciata ssa al valore iniziale A = @g@x(~x) : 80 (2.7.18) Si puo dimostrare che cio provoca una distorsione del 2 ordine nella stima. Lo schema iterativo e bloccato qundo il ^oi2 non diminuisce piu che di una piccola percentuale ^oi2 ; ^oi2 +1 < " pressato : ^oi2 (2.7.19) 2.8 Applicazione alla regressione lineare Supponiamo di avere un processo descrivibile mediante una (o piu) variabili y su cui compiamo delle osservazioni Yoi che comprendono una parte stocastica e un valore medio yi che si sa essere funzione di altre p variabili t1 t2 : : : tp completamente note per ognuna delle i = 1 : : : n osservazioni Yoi = yi + "i yi = g(t1i t2i : : : tpi) E f"ig = 0 C"" = o2I : (2.8.1) (2.8.2) Cio che importa sottolineare nel modello (2.8.1), (2.8.2) e che la parte stocastica va considerata aggiunta, per ogni osservazione i, alla variabile dipendente y, detta anche variabile criterio, mentre le variabili indipendenti t1 : : : tp, per la stessa osservazione, sono considerate costanti, cioe comunque note con una precisione tale da non modicare signicativmente il modello stocastico di " le variabili t1 : : : tp sono dette variabili concomitanti o di predizione. Spesso siamo in una situazione di analisi puramente empirica dei dati in cui il legame y = g(t1 : : : tp) non e noto perche troppo complesso: tuttavia da una prima indagine dei dati stessi o da una conosenza qualitativa del fenomeno sotto indagine possiamo ritenere che un'approssimazione lineare di g sia ben giusticata almeno nell'ambito dei valori assunti dalle variabili concomitanti in corrispondenza alle nostre osservazioni. Si potra cos ragionevolmente sostituire alla II delle (2.8.1) l'equazione 81 8 p X > @g > (tki ; tk ) < yi = g(t1 : : : tp) + @t k k=1 : n X > 1 > : tk = n i=1 tki (2.8.3) Osservazione 2.8.1: si noti che l'uso delle medie aritmetiche tk per linearizzare g, e naturale per ottimizzare la validita della approssimazione lineare stessa. Per semplicita nei conti successivi, sebbene cio non sia strettamente necessario, supporremo di aver gia scelto l'origine delle variabili concomitanti in modo tale che si abbia il baricentro dei punti d'osservazione nell'origine tk = 0 k = 1 : : : p (2.8.4) qualora tale condizione non fosse gia soddisfatta bastera una semplice traslazione per ottenere nuove variabili indipendenti ki = tki ; tk che invece soddisferanno la (2.8.4). Supponendo dunque che valga la (2.8.4), il modello (2.8.3) diventa yi = co + c1 t1i + : : : + cptpi (2.8.5) i parametri incogniti di questo modello sono co c1 : : : cp che organizzeremo nel vettore xo x c.1 (xo = co x = .. c p ) : In forma matriciale la (2.8.5) puo essere riscritta come y = exo + Tx 82 (2.8.6) dove 1 . t. 11 : : : e = .. T = .. 1 t1n : : : tpi ... : (2.8.7) tpn Il modello deterministico (2.8.6) e quello stocastico (2.8.2) insieme deniscono un tipico problema di minimi quadrati lineari in forma parametrica, la cui compensazione fornira le stime ricercate x^o x^. Osservazione 2.8.2: notiamo che in virtu della condizione (2.8.4) si ha t.1 + T e = n .. tp = 0 cos che (Tx)+(exo) = 0 : (2.8.8) Cio signica che nel modello (2.8.6) la varieta dei valori ammissibili (ovviamente a p +1 dimensioni) e decomposta in due sottovarieta ortogonali fexog fTxg rispettivamente a 1 e a p dimensioni. Se ora formiamo il sistema normale per la (2.8.6), indicato al solito con Yo il vettore delle osservazioni e osservato che nella (2.8.6) non si ha termine noto a, si trova + e e +0 0 T T : (2.8.9) 8 n X > + > (1=n)(T T )kl = (1=n) tkitli = Ctt > > > < n i=1 X (1=n)e+Yo = (1=n) Yoi = Y > > > i=1P > > (1 =n ) t Y 1 i oi > + i = Cty P : (1=n)T Yo = (2.8.10) x^o x^ + = e +Yo T Yo Si ha e+e = n. Poniamo inoltre (1=n) 83 i tpi Yoi e notiamo che queste relazioni sembrano proprio quelle di media, covarianza e cross-varianza di valori campionari, motivo questo della simbologia adottata: tuttavia e necessario sottolineare che ne Ctt Cty possono Figura 2.8.1: essere covarianze, ne Y e una media, perche da un lato ftkig non sono campioni estratti da variabili casuali tk (le variabili concomitanti sono solo costanti note esattamente), dall'altro lato le componenti Yoi non sono estrazioni dalla stessa v.c. Y in quanto le medie di Yoi hanno valori diversi. Con questa simbologia il sistema normale (2.8.9), dividendo entrambi i membri per n, da la soluzione x^o = Y x^ = Ctt 1Cty ; con matrice di covarianza 2 (^xo ) = 2 o n 84 (2.8.11) 2 Cx^x^ = o2 (nCtt ) 1 = no Ctt 1 Cx^ox^ = 0 : ; Per la stima di 2 o ; formiamo dapprima il vettore degli scarti U = Yo ; Y e ; T x^ : ricordando la (2.8.8) e la (2.8.9) si ha jU j 2 (2.8.12) (2.8.13) = U + U = jYo ; Y ej2 ; 2^x+T +Yo + x^+ T +T x^ = = jYo ; Y ej2 ; jT x^j2 = jYo ; Y ej2 ; nx^+ Cty = X X (2.8.14) = (Yoi ; Y )2 ; c^k tkiYoi : i ki Da qui in poi si trova ^o2 secondo la regola generale ^o2 = n ; 1p ; 1 U +U : Osservazione 2.8.3: l'equazione (2.8.14) esprime semplicemente la decomposizione pitagorica di g. 2.8.2 e fornisce uno strumento comodo per la stima di o2. Inoltre, notando che Yoi ; Y puo essere interpretato P \formalmente" come scarto generale di Yoi dallaPmedia Y (T x^)i = pk=1 tkic^k come scarto tra la regressione (^yi = c^o + k tkic^k ) e la media Y (^co = Y ), cioe come lo scarto spiegato dall'equazione di regressione, Ui = Yoi ; yi come scarto residuo tra il valore osservato ed il valore della regressione (cfr. g. 2.8.3), si puo riscrivere la (2.8.14) nella forma SG2 = SR2 + SS2 con P SG2 = P(Yoi ; Y )2 = scarti generali SR2 = P (Yoi ; yi)2 = scarti residui SS2 = (yi ; Y )2 = scarti spiegati 85 (2.8.15) Figura 2.8.2: Figura 2.8.3: Osservazione 2.8.4: si noti che nel caso di una sola variabile concomitante, la retta di regressione stimata y^ = c^o + c^1 t diventa, secondo le (2.8.11) (t = 0) y^ = Y + ty 2t t che coincide formalmente con la retta di regressione per una variabile doppia (Y t), sebbene in questo caso l'interpretazione sia notevolmente diversa. Osservazione 2.8.5: in generale il motivo per cui interessa costruire il modello (2.8.5) non e tanto quello di stimare co c1 : : : cp e quindi 86 anche i valori della regressione y^i per, diciamo cos, regolarizzare i valori osservati Yoi, ma piuttosto quello di poter prevedere pur con un margine di errore, il valore della regressione y^ in corrispondenza di valori delle variabili concomitanti per i quali non si hanno osservazioni. Naturalmente il valore predetto sara quello della regressione y^(t1 : : : tp) = c^o + p X k=1 c^k tk che, introdotto il vettore t.1 t = .. tp puo anche essere scritto come y^(t) = x^o + x^+t : (2.8.16) Ricordando le (2.8.12), possiamo accompagnare il valore predetto con la sua varianza (ovvero l'errore di predizione) 2 o 1 + t+ C 1 ^y(t)] = n tt t] : nel caso di una sola variabile t, la (2.8.17) assume la forma 2 2 6 2 6 o 2 ^y(t)] = n 661 + 4 ; 3 7 t2 77 : n 7 X 25 (1=n) i=1 ti (2.8.17) (2.8.18) Per stimare queste varianze ovviamente occorrera avvalersi del valore stimato ^o2. Dalla (2.8.18) appare in modo piu chiaro che l'errore di predizione si mantiene limitato nche t rimane all'interno dei valori di osservazione 87 Figura 2.8.4: (jtj max jtij), mentre esso aumenta rapidamente quando si tenti di estrapolare la regressione all'esterno dei valori di osservazione. Esempio 2.8.1: nelle operazioni di controllo di una struttura, si misurano le quote di un punto per 6 volte con intervalli di due mesi per vericare gli eventuali movimenti di assestamento del punto stesso. Poiche pero si sa che la forma geometrica della struttura dipende anche dalla sua temperatura, contestualmente alla misura della quota si fa un lettura di temperatura. Si ipotizza che valga un modello di regressione lineare del tipo Q = co + c1 t + c2T (Q = quota t = tempo T = temperatura) Qoi = Qi + "i C"" = o2 I : I dati ottenuti sono: t (in mesi) T (in C ) Q (in mm) (Yo) 2 4,8 127,02 4 9,3 129,22 6 12,5 130,98 8 19,2 133,74 10 10,9 129,98 12 5,2 127,01 Medie: t=7 T = 10 3167 Y = 129 6583 88 Si osservi in primo luogo che ne t ne T sono centrate, percio e utile passare a nuove variabili centrate ponendo t=t+ T =T + e riformando il modello come y = Q = co + c1 + c2. Si ottiene cos la nuova tabella Yo ; Y -5 -5,5167 -2,6383 n=6 -3 -1,0167 -0,4383 -1 2,1833 1,3217 p=2 1 8,8833 4,0817 3 0,5833 0,3217 (m = 3) 5 -5,1167 -2,6483 2 2 2 S = 70 S = 141 6683 SG = 32 6770 grad. lib. S = 13 500 SYo = 4,9900 di ^o2 = 3 SYo = 67 8832 Da cui 11 6 2 2500 0873 1 Ctt = ;2 2500 23 6114 Ctt = ;00 0082 0 8316 : Cty = ; 11 3139 Dunque le stime sono c^o = x^o = 129 6583 c^1 = x^ = ;0 021535 : c^2 0 481227 Inoltre, ricordando le (2.8.14), (2.8.15) 89 ;0 0082 0 0431 SR2 = SG2 ; 6 x^+Cty = 0 1172 : cos che ^o2 = 0 0391 ^o = 0 1977 : In tal modo si ha anche la stima di varianza e covarianza di xo x, cioe (^xo ) = 0 0065 0 0873 Cx^x^ = 0 0065 ;0 0082 2 ;0 0082 0 0431 : in particolare e (^c1) = 0 0238 (^c2) = 0 0167 : Come si vede c^1 e (^c1) sono quasi uguali, il che fa pensare ad una dipendenza da (cioe da t) non molto signicativa resta invece sensibile la dipendenza da (cioe da T ). Supponiamo ora di voler predire il valore di y in corrispondenza a t = 7 T = 15 , cioe = 0 = 4 6833: si ha y(0 4 6833) = c^o + c^1 0 + c^2 4 6833 = 131 9120 : Inne, volendo calcolare l'errore di stima di tale predizione, si ha 2 ovvero 0 (^y) = 6 1 + 0 4 6833]Ctt 4 6833 = 2 o 1 ; 0 0761 = 0 0127 = 0 0391 f 1 + 0 9453 g = 6 6 (^y) = 0 11 mm : 90