Indice 1 Il metodo della massima verosimiglianza

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