Appunti di teoria della stima Parte prima: approccio...

Universita degli Studi di SienaDipartimento di Ingegneria dell’InformazioneVia Roma 56, 53100 Siena, Italy

Elaborazione Numerica dei Segnali

Appunti di teoria della stimaParte prima: approccio classico

Prof. Andrea Garzelli

Ing. Luca Capobianco

Appunti tratti dal testo S.M.Kay ‘Fundamentals of Statistical Signal Processing: Esti-mation Theory’, al quale lo studente dovra fare riferimento per una trattazione piu com-pleta.

Indice

Introduzione 1

1 Stimatore non polarizzato a minima varianza (MVU) 81.1 Stimatore non polarizzato . . . . . . . . . . . . . . . . . . . . 81.2 Prestazioni: criterio a minima varianza . . . . . . . . . . . . . 101.3 Esistenza dello stimatore MVU . . . . . . . . . . . . . . . . . 111.4 Ricerca dell’MVU . . . . . . . . . . . . . . . . . . . . . . . . . 13

2 Cramer-Rao Lower Bound 152.1 Considerazioni sull’accuratezza di uno stimatore . . . . . . . . 152.2 Cramer Rao Lower Bound . . . . . . . . . . . . . . . . . . . . 16

Teorema Cramer-Rao . . . . . . . . . . . . . . . . . . . . . . . 16- Esempio 2.1 - Stimatore non polarizzato per livello di corrente

continua in rumore WGN . . . . . . . . . . . . . . . . . 172.3 CRLB per un segnale immerso in rumore bianco Gaussiano . . 19

- Esempio 2.2 - -Stima della frequenza di un segnale sinusoidale 202.4 Trasformazione dei Parametri . . . . . . . . . . . . . . . . . . 202.5 Stima di piu parametri - caso vettoriale . . . . . . . . . . . . . 22

- Esempio 2.3 - DC level in WGN . . . . . . . . . . . . . . . . 22- Esempio 2.3 - Line fitting . . . . . . . . . . . . . . . . . . . 23Teorema Cramer-Rao - Caso Vettoriale . . . . . . . . . . . . . 24

3 MVU per modelli lineari 253.1 Definizione e proprieta . . . . . . . . . . . . . . . . . . . . . . 25

Teorema MVU per modello lineare . . . . . . . . . . . . . . . . 25- Esempio 3.1 - Curve fitting . . . . . . . . . . . . . . . . . . 26- Esempio 3.2 - Sistema di identificazione . . . . . . . . . . . . 27

4 General MVU Estimation 324.1 Statistiche sufficienti . . . . . . . . . . . . . . . . . . . . . . . 324.2 Come trovare la Statistica Sufficiente . . . . . . . . . . . . . . 34

i

Dipartimento di Ingegneria dell’Informazione ii

Teorema di Neymann-Fisher . . . . . . . . . . . . . . . . . . . 34- Esempio 4.1 - DC level in WGN . . . . . . . . . . . . . . . . 34- Esempio 4.2 - Power of WGN . . . . . . . . . . . . . . . . . 35- Esempio 4.3 - Fase di una sinusoide . . . . . . . . . . . . . . 35Teorema Rao-Blackwell-Lehmann-Scheffe . . . . . . . . . . . . 36Procedura per trovare lo stimatore MVU . . . . . . . . . . . . 37- Esempio 4.4 - DC level in WGN . . . . . . . . . . . . . . . . 37

4.3 Legame fra l’unicita di g e la completezza della statistica . . . 39- Esempio 4.5 - Incompleta statistica sufficiente . . . . . . . . . 40- Esempio 4.6 - Media di un rumore bianco . . . . . . . . . . . 42Teorema Neymann-Fisher - Caso Vettoriale . . . . . . . . . . . 45Teorema Rao-Blackwell-Lehmann-Scheffe - Caso Vettoriale . . . 45

5 Best Linear Unbiased Estimator (BLUE) 465.1 Definizione e Vincoli per il BLUE . . . . . . . . . . . . . . . . 46

- Esempio 5.1 - DC level in white noise . . . . . . . . . . . . . 51- Esempio 5.2 - DC level in uncorrelated zero-mean noise . . . . 52

5.2 BLUE - Caso Vettoriale . . . . . . . . . . . . . . . . . . . . . 53Teorema di Gauss-Markov . . . . . . . . . . . . . . . . . . . . 54

6 Stima a massima verosimiglianza (MLE) 556.1 MLE - Stimatori Consistenti (caso scalare) . . . . . . . . . . . 55

- Esempio 6.1 - DC level in white Gaussian Noise . . . . . . . . 556.2 Trovare lo stimatore MLE . . . . . . . . . . . . . . . . . . . . 57

Teorema - Proprieta asintotiche dell’MLE . . . . . . . . . . . . 59- Esempio 6.2 - MLE per la stima di una fase di una sinusoide . 59Teorema - Proprieta di invarianza dell’MLE . . . . . . . . . . . 62

6.3 Valutazione numerica dell’MLE . . . . . . . . . . . . . . . . . 62- Esempio 6.3 - Metodo di Newton-Raphson - Metodo di Scoring 63

7 Metodo dei minimi quadrati (Least Squares) 677.1 Generalita sul problema dei minimi quadrati . . . . . . . . . . 67

- Esempio 7.1 - Segnale sinusoidale . . . . . . . . . . . . . . . 697.2 Minimi quadrati lineari - caso scalare . . . . . . . . . . . . . . 697.3 Minimi quadrati lineari - caso vettoriale . . . . . . . . . . . . 707.4 Interpretazione Geometrica . . . . . . . . . . . . . . . . . . . . 72

- Esempio 7.2 - Sovrapposizione di sinusoidi . . . . . . . . . . . 727.5 Principali tecniche iterative LLS . . . . . . . . . . . . . . . . . 757.6 Problema LS non lineare . . . . . . . . . . . . . . . . . . . . . 79

- Esempio 7.3 - Digital Filter Design . . . . . . . . . . . . . . 80- Stima AR di parametri per un modello ARMA . . . . . . . . 84

Introduzione

La teoria della stima rappresenta l’argomento chiave per molti dei piu at-tuali sistemi di signal processing, che vengono utilizzati per l’elaborazione el’estrazione di informazione da un insieme di dati.Esempi di tali sistemi possono essere:

• Radar

• Sonar

• Riconoscimento vocale

• Analisi di immagini

In molte applicazioni e necessario stimare il valore di uno o piu parametrisignificativi da un set di dati, i quali possono dipendere da un insieme difattori il piu delle volte sconosciuti o caratterizzabili solo per via statistica.Matematicamente, il problema puo essere caratterizzato nel modo seguente:si ha un vettore di dati di N elementi dipendenti da un parametro sconosciutoθ. Vogliamo determinare θ utilizzando i dati a disposizione, che in altritermini significa definire uno stimatore

θ = f(x[0], x[1], . . . , x[N − 1]) (1)

dove f indica una funzione e x e il vettore dei dati.Ad esempio, nel caso del Radar, l’interesse potrebbe essere rivolto a deter-minare la posizione dell’oggetto in esame, a partire dal ritardo τ0 al quale siriceve l’eco dell’oggetto, mediante lo studio dell’equazione

τ0 =2R

c.

Esistono pero delle difficolta, dovute a fenomeni aleatori, come perdite dipropagazione e disturbi di diversa natura dovute a sorgenti di rumore, ri-tardi aggiuntivi dovuti alle distorsioni di canale. La teoria della stima puo

1

Dipartimento di Ingegneria dell’Informazione 2

essere d’aiuto, permettendo l’approssimazione τ0 mediante un valore τ0.Allo stesso modo, in un sistema sonar di tipo passivo, in cui il parametrodi interesse e ancora la posizione di un oggetto in esame, si stima il valoredell’angolo di vista β (definito da β = arccos(vτ0

d)) mediante il valore β.

Per determinare uno stimatore che produca risultati affidabili, e necessariocome primo passo effettuare una buona modellizzazione matematica dei datiin esame. Tale modellizzazione deve essere effettuata mediante una descri-zione probabilistica, e piu precisamente attraverso una densita di probabilita(PDF), ovvero:

p(x; θ) = p(x[0], x[1], . . . , x[N − 1]; θ) . (2)

Dalla 2 si vede come la PDF venga parametrizzata tramite la variabile sco-nosciuta θ; cio che si ottiene e una famiglia di PDF, ognuna per un datovalore di θ.

Per illustrare un esempio tipico ed abbastanza esplicativo, consideriamol’andamento della curva in figura 1 (estrapolata da una media dell’indiceDow-Jones): nonostante sia molto fluttuante, l’andamento segue un aumento

Figura 1: Media Dow Jones

medio piuttosto costante, apparentemente lineare.Per valutare la veridicita di questa ipotesi di linearita, possiamo assumereche i dati possano essere descritti da una somma di una retta e di un rumore,ovvero:

x[n] = A + Bn + w[n] n = 0, 1, . . . , N − 1 . (3)


Un assunto ragionevole che riguarda il rumore e che esso sia WGN, bianco eGaussiano, ovvero che ogni campione di w[n] segua una PDF Gaussiana, conmedia nulla e varianza σ2, e che sia scorrelato con tutti gli altri campioni.Dunque, i valori sconosciuti sono a questo punto A e B, che possono essereespressi in un unico vettore θ = [AB]T . Sia x = [x[0], x[1], . . . , x[N − 1]],allora la PDF puo essere espressa come

p(x; θ) =1

(2πσ2)N2

exp[− 1

2σ2

N−1∑n=0

(x[n]− A−Bn)2] . (4)

L’assunzione di rumore WGN e giustificata dalla necessita di formulare unmodello matematico in modo che possa essere trovata un’espressione dellostimatore in forma chiusa. Ovviamente l’affidabilita di qualsiasi stimatore efortemente dipendente dalla scelta della PDF.Una volta specificata la PDF, il problema si sposta sulla possibilita di speci-ficare uno stimatore ottimo per la curva dei dati: uno stimatore puo infattidipendere da altri parametri, a condizione pero che questi siano noti. E’ beneporre l’attenzione sulla differenza che intercorre fra lo stimatore e una stima:una stimatore puo essere pensato come la regola che assegna un valore a θper ogni realizzazione di x: invece, la stima di θ, e il valore di θ ottenuto peruna data realizzazione di x.Vediamo adesso, con l’aiuto di un esempio, come sia possibile valutare leprestazioni di uno stimatore. Consideriamo i dati visibili in figura 2. Gia daun primo sguardo, appare come il segnale x[n] sia composto da un segnale incontinua a cui si sovrappone un certo livello di rumore. Si ha cioe:

x[n] = A + w[n] (5)

dove w[n] e un processo di rumore a media nulla. Basandoci su questa ipotesi,il nostro interesse e rivolto ad una stima di A. Intuitivamente, dato che Arappresenta il livello medio di x[n], una scelta ragionevole per stimare A puoessere:

A =1

N

N−1∑n=0

x[n] . (6)

Occorre subito chiedersi:

• Con che precisione A approssima A?

• Esistono stimatori migliori di quello proposto?

Per i dati di figura 2, si ha A = 1.1, che e molto vicino al valore reale utiliz-zato A = 1.


0 10 20 30 40 50 60 70 80 90 100−2

−1

0

1

2

3

4

n

x[n]

Figura 2: Corrente DC con rumore


Un altro stimatore per lo stesso problema potrebbe anche essere A = x[0];tuttavia e lecito aspettarsi che non sia affidabile, in quanto non fa uso ditutti i dati e della informazione in essi contenuta, ed inoltre non utilizzauna operazione di media che riduca l’effetto del rumore. In particolare, inquesto caso, A = 2.1, valore molto lontano dalla media reale. Inoltre, purammettendo che A abbia un valore vicino al valore reale, non potremmocertamente concludere che esso sia uno stimatore affidabile: uno stimatoredeve essere una funzione dei dati in esame, che sono variabili casuali, ed epertanto anch’esso una variabile casuale. Pertanto, per stimare l’affidabilitadi uno stimatore, dobbiamo utilizzare metodi di natura statistica. Per ren-dere piu chiaro questo concetto, immaginiamo di ripetere l’esperimento dacui sono stati tratti i dati, e di applicare i due stimatori proposti ai vari setdi dati ottenuti. Supponiamo inoltre di fissare A = 1 e di sommare differentirealizzazioni di rumore w[n]; riportando su un istogramma i valori dei duestimatori proposti per ogni insieme di dati, possiamo ricostruire una appros-simazione della PDF che descrive il numero di volte che uno stimatore haprodotto un certo valore. I grafici nelle figure 3 e 4 riportano gli istogram-

Figura 3: Istogramma di A

mi relativi a 100 realizzazioni. E’ evidente che A e uno stimatore migliore,perche i valori ottenuti sono piu concentrati attorno al valore A = 1. Perprovare che cio sarebbe vero per qualsiasi numero di realizzazioni, possiamopassare alla valutazione statistica degli stimatori, ovvero ad una valutazionedei valori di media e varianza dello stimatore.Assumendo che il rumore abbia media nulla e varianza σ2, si ha

E(A) = E(1

N

N−1∑n=0

x[n]) =1

N

N−1∑n=0

E(x[n]) = A (7)

E(A) = E(x[0]) = A , (8)


Figura 4: Istogramma di A

dunque in media gli stimatori producono il valore A. La varianza deglistimatori e, essendo w[n] scorrelato:

var(A) = var(1

N

N−1∑n=0

x[n]) =1

N2

N−1∑n=0

var(x[n]) =1

N2Nσ2 =

σ2

N(9)

var(A) = var(x[0]) = σ2 > var(A) . (10)

L’esempio sottolinea due aspetti importanti della teoria della stima che de-vono essere tenuti ben chiari in mente:

• Uno stimatore e una variabile casuale e pertanto puo essere completa-mente descritto solo statisticamente o dalla sua PDF

• L’uso di una simulazione al computer puo non essere completa: la stimadelle prestazioni puo al massimo essere valutata con un certo gradodi accuratezza o, nel peggiore dei casi, si possono ottenere risultatisbagliati, dovuti ad un insufficiente numero di esperimenti o ad unacattiva caratterizzazione del problema.

Abbiamo anticipato che dalla scelta della PDF dipende l’accuratezza el’affidabilita dello stimatore: vogliamo sottolineare inoltre che tale scelta delmodello deve ovviamente essere guidata dall’analisi di alcune caratteristichefondamentali che la PDF deve presentare:

• compatibile con i vincoli del problema

• compatibile con la conoscenza a priori dei dati

• matematicamente maneggevole.


La teoria della stima si fonda su due tipi di approcci fondamentali, classico eBayesiano. Nel primo, i parametri di interesse vengono assunti sconosciuti,ma deterministici, e come tali vengono ricercati.Nell’esempio relativo all’indice Dow Jones, si vede dai dati come la mediasi aggiri all’incirca intorno al valore 3000(figura 1). Sarebbe pertanto inu-tile cercare un valore di A inferiore a 2000 o superiore a 4000. Possiamoalora restringere la ricerca all’intervallo [2800, 3200], introducendo questa in-formazione nota nella PDF di A: ovvero possiamo assumere che A non siadeterministica, ma una variabile con una sua PDF, possibilmente uniformenell’intervallo [2800, 3200]. In un approccio del genere, detto Bayesiano, ilparametro che stiamo tentando di stimare, viene visto come una realizzazionedi una variabile casuale, le cui informazioni a priori sono contemplate nellarelativa PDF:

p(x, θ) = p(x|θ)p(θ) (11)

dove la p(θ) e la PDF che riassume le informazioni riguardo θ, note prima del-l’osservazione dei dati, mentre la p(x|θ) e la PDF che riassume la conoscenzaa posteriori, offerta dall’osservazione dei dati condizionati da θ.

Capitolo 1

Stimatore non polarizzato aminima varianza (MVU)

Il capitolo presenta una prima analisi degli stimatori che, in media, produ-cono il valore vero del parametro in esame. Fra questi, cercheremo quelloa minore variabilita. Vedremo come uno stimatore cosı ottenuto, produrravalori prossimi al valore vero nella maggior parte delle applicazioni.Gli stimatori che presentano la caratteristica appena descritta e riassuntamatematicamente dall’equazione 1.1, vengono denominati stimatori non po-larizzati. Non sempre esiste lo stimatore non polarizzato a minima varianza(MVU, minimum variance unbiased). Qualora esso esista possiamo affidar-ci a diversi metodi per trovarlo: fra questi citiamo il metodo lower boundCramer-Rao, esaminato nel corso del capitolo; qualora invece non esista, oin caso di fallimento del metodo, si possono utilizzare dei vincoli lineari perottenere una semplice implementazione subottima dello stimatore.

1.1 Stimatore non polarizzato

Perche uno stimatore sia non polarizzato, si richiede che esso in media pro-duca il valore vero del parametro sconosciuto in esame. In altri termini, lacaratteristica di non polarizzazione afferma che non ha importanza quale siail valore di θ, in quanto il nostro stimatore lo produrra in media, per qualun-que θ ∈ [a, b].Uno stimatore e non polarizzato se:

E(θ) = θ ∀ θ ∈ [a, b]. (1.1)

Cominciamo ad affrontare l’argomento partendo da un semplice esempio.

8


Esempio 1.1 - Stimatore non polarizzato per livello dicorrente continua in in rumore WGN

Si consideri l’osservazione:

x[n] = A + w[n] n = 0, 1, ..., N − 1 (1.2)

dove A e il parametro che dev’essere stimato, mentre w[n] e un vettore dirumore WGN. Il valore di A puo essere compreso nell’intervallo −∞ < A <∞. Dunque, una scelta ragionevole per la valutazione del valore medio e

A =1

N

N−1∑n=0

x[n] . (1.3)

Per le proprieta di linearita dell’operatore E[·], come gia illustrato nell’intro-duzione abbiamo:

E(A) = E

[1

N

N−1∑n=0

x[n]

]=

1

N

N−1∑n=0

E(x[n]) = A (1.4)

per qualsiasi A.Consideriamo invece, per la stessa osservazione in 1.2, lo stimatore

A =1

2N

N−1∑n=0

x[n] . (1.5)

Applicando l’operatore di media, si ha

E(A) = E[1

2A]

= A se A = 0

6= A se A 6= 0

Lo stimatore modificato funziona solo se A = 0, pertanto esso e polarizzato.

La proprieta di non polarizzazione di uno stimatore ha una importante im-plicazione quando diversi stimatori vengano combinati fra loro: puo suc-cedere infatti che per uno stesso parametro siano disponibili piu stimatori{θ1, θ2, ..., θn}. Si puo ottenere da questi stimatori, uno presumibilmentemigliore, effettuando una media su di essi cioe vale a dire:

θ =1

n

n∑i=1

θi . (1.6)


L’importanza dell’assunzione che gli stimatori siano tutti non polarizzati, constessa varianza, e scorrelati a due a due, ci porta a poter concludere che

E(θ) = θ (1.7)

var(θ) =1

n2

n∑i=1

var(θi) =var(θi)

n. (1.8)

Come si puo vedere, se n cresce, la varianza dello stimatore cosı ottenutodecresce e passando al limite, se n →∞, θ → θ, con varianza nulla.Tuttavia, se gli stimatori sono polarizzati, ovvero E(θi) = θ+b(θ) - dove b(θ)e il valore di polarizzazione (bias) dello stimatore allora

E(θ) =1

n

n∑i=1

E(θi)

= θ + b(θ)

θ non converge al valore vero, indipendentemente dal numero di stimatoriutilizzati nella media, cioe

limn→∞

θ 6= θ

1.2 Prestazioni: criterio a minima varianza

Nella ricerca di uno stimatore dobbiamo ovviamente stabilire un criterio di ot-timalita per valutarne le prestazioni. Possiamo utilizzare l’errore quadraticomedio (MSE), definito come

mse(θ) = E[(θ − θ)2] . (1.9)

Questo indice misura la deviazione quadratica media dello stimatore dal va-lore vero. In generale l’adozione di questo criterio porta a stimatori che nonpossono essere scritti esclusivamente in funzione dei dati. Per capire cosa ciosignifichi, riscriviamo l’mse come

mse(θ) = E{[θ − E(θ) + E(θ)− θ]2}= var(θ) + [E(θ)− θ]2

= var(θ) + b2(θ) (1.10)

da cui si evidenzia come l’MSE contenga un errore dovuto alla varianza deglistimatori cosı come al valore di polarizzazione (bias).


Esempio 1.2

Consideriamo, nel problema dell’esempio 1.1, lo stimatore modificato

A = a1

N

N−1∑n=0

x[n]

dove cercheremo di determinare la costante a in modo che risulti minimol’mse.Siccome E(A) = aA e var(A) = a2σ2/N , abbiamo dalla 1.10,

mse(A) =a2σ2

N+ (a− 1)2A2 .

Derivando rispetto ad a otteniamo,

mse(A)

da= 2

aσ2

N+ 2(a− 1)A2

che uguagliata a zero fornisce un valore di minimo con

aopt =A2

A2 + σ/N.

Per valori finiti di N , il valore ottimo di a dipende dal valore sconosciuto diA, pertanto lo stimatore non e realizzabile.

1.3 Esistenza dello stimatore MVU

Vogliamo a questo punto stabilire un criterio per valutare se esista uno sti-matore non polarizzato a minima varianza (MVU) per ogni θ appartenenteall’intervallo di esame.Illustriamo, nel prossimo esempio, come in generale l’MV U non sempre esi-sta. Per capire in che modo questo sia possibile, osserviamo i due casi visibiliin figura 1.1. Nel caso della figura 1.1a, risulta ovvio che lo stimatore MVU equello individuato dalla curva θ3; tuttavia, nel caso in figura 1.1b non si puoindividuare uno stimatore ottimo, in quanto per valori di θ inferiori a θ0, θ2

e migliore, mentre per valori θ > θ0, lo stimatore θ3 risulta migliore secondoil criterio MVU. Spesso, per riferirsi a situazioni simili al caso in figura 1.1adove non si prsentano incertezze, θ3 viene chiamato stimatore non polarizzatoa varianza uniforme minima.


Figura 1.1: Esempio di inesistenza di uno stimatore MVU

Esempio 1.3 - Inesistenza dell’MVU

Se la funzione densita di probabilita (PDF) cambia al variare di θ, possiamoaspettarci che con θ possa cambiare anche lo stimatore da utilizzare.Date due successive osservazioni, x[0] e x[1], con PDF

x[0] ∼ N (θ, 1)

x[1] ∼{ N (θ, 1) se θ ≥ 0N (θ, 2) se θ < 0

dove N indica una PDF normale.Si puo mostrare che gli stimatori

θ1 =1

2(x[0] + x[1])

θ2 =2

3x[0] +

1

3x[1]

non sono polarizzati. La varianza risulta essere rispettivamente

var(θ1) =1

4(var(x[0]) + var(x[1]))

var(θ2) =4

9var(x[0]) +

1

9var(x[1])

per cui,

var(θ1) =

{1836

= 12

se θ ≥ 02736

= 34

se θ < 0

var(θ2) =

{2036

= 59

se θ ≥ 02436

= 23

se θ < 0

La figura 1.2 mostra il risultato. Si puo dimostrare che in generale per θ ≥ 0il valore minimo possibile per la varianza di un stimatore non polarizzato e18/36, mentre per θ < 0 e 24/36.


Figura 1.2: Esempio di inesistenza di uno stimatore MVU

1.4 Ricerca dell’MVU

A questo punto risulta piu chiaro che la ricerca di uno stimatore non e unproblema di semplice soluzione, in quanto non sempre esiste lo stimatoreMVU. Inoltre, anche qualora esista, non esiste un metodo che dia la certezzadi poterlo trovare.In generale, possiamo dire che abbiamo tre possibilita, illustrate in dettaglionei capitoli successivi:

• Determinare il limite inferiore col metodo Cramer-Rao (CRLB, Cramer-Rao lower bound) e verificare l’esistenza di uno stimatore che lo soddisfi;

• Applicare il teorema Rao-Blackwell-Lehmann-Scheffe, (RBLS);

• Ricavare una classe di stimatori lineari e non-polarizzati dallo stimatoreiniziale, e successivamente trovare in questo insieme l’MVU.

Ovviamente, mentre il primo ed il secondo metodo possono produrre lo sti-matore MVU, il terzo lo produrra solo se questo e lineare nei dati.Il metodo CRLB ci permette di affermare che dato un insieme di dati, non

Figura 1.3: Esempio di CRLB


e possibile trovare uno stimatore non polarizzato con varianza strettamenteminore di un valore che rappresenta il limite inferiore per quel dato problemadi stima. Pertanto qualora trovassimo uno stimatore con varianza uguale alvalore dato dal metodo CRLB, potremmo affermare con certezza che quelloe l’MVU che stiamo cercando.

Capitolo 2

Cramer-Rao Lower Bound

La possibilita di determinare un limite inferiore (lower bound) per la varianzadi un qualsiasi stimatore non polarizzato risulta essere estremamente utilenella pratica. Questo infatti ci permette di poter analizzare le performancedello stimatore in esame da un confronto diretto fra i valori teorici e i valoriottenuti. Nel migliore dei casi, questo metodo ci permette di identificare lostimatore MVU.Nonostante esistano diverse tecniche per determinare il lower bound, il meto-do Cramer-Rao (CRLB) analizzato in questo capitolo e il piu utilizzato perla sua semplicita di calcolo.

2.1 Considerazioni sull’accuratezza di uno sti-

matore

In generale l’informazione necessaria alla stima e rappresentata dai dati os-servati e dalla PDF che li caratterizza; pertanto e abbastanza ovvio comel’accuratezza della stima dipenda direttamente dalla PDF. Ad esempio, sela PDF dipende solo debolmente dal parametro che stiamo cercando di sti-mare o addirittura non dipende affatto da esso, non possiamo certamenteaspettarci una corretta stima di tale parametro. In generale, quanto piu laPDF e influenzata dal parametro sconosciuto, maggiore e la precisione concui possiamo effettuare la stima.Se la PDF e vista come una funzione del parametro sconosciuto - con x fissato- viene denominata funzione di verosimiglianza (likelihood function). Intui-tivamente, la precisione con cui possiamo stimare il parametro dipende dallapresenza di variazioni brusche della funzione di verosimiglianza: maggiore el’irregolarita della funzione di verosimiglianza (intuitivamente: grandi valoridi curvatura), maggiormente sara accurata la stima.

15


Per meglio formalizzare questo concetto, osserviamo che l’irregolarita dellafunzione puo essere misurata dal valore massimo, cambiato di segno, delladerivata seconda del logaritmo della funzione di verosimiglianza, ovvero dallacurvatura media della log-likelihood function. La misura della curvatura puoessere fornita dalla funzione

− E

[∂2 ln p(x; A)

∂A2

](2.1)

che e appunto una misura della curvatura media della log-likelihood function.Maggiore risulta la quantita espressa dalla 2.1, minore risultera la varianzadello stimatore.

2.2 Cramer Rao Lower Bound

Possiamo adesso enunciare il teorema che e alla base del metodo CRLB.

Teorema 1 (Cramer Rao Lower Bound - parametro scalare) Si as-suma che la PDF p(x; θ) soddisfi la condizione di regolarita, ovvero

E

[∂ ln p(x; θ)

∂θ

]= 0 ∀θ.

La varianza di un qualsiasi stimatore θ non polarizzato deve soddisfare ladisequazione

var(θ) ≥ 1

−E[

∂2 ln p(x;θ)∂θ2

] ∀θ. (2.2)

dove la derivata e valutata sul valore vero di θ, e il valore atteso su p(x; θ).Inoltre e possibile trovare uno stimatore non polarizzato per cui valga l’ugua-glianza per ogni θ, se e solo se la derivata rispetto a θ della log-likelihoodfunction, ovvero la score function soddisfa l’uguaglianza

[∂ ln p(x; θ)

∂θ

]= I(θ)(g(x)− θ) (2.3)

per qualche coppia di funzioni g e I. Tale stimatore, che risulta essere propriol’MVU, e θ = g(x) e la funzione che esprime il valore di varianza minima e1/I(θ).

L’espressione esplicita del valore atteso presente nell’espressione 2.2 e datada:

E

[∂2 ln p(x; θ)

∂θ2

]=

∫∂2 ln p(x; θ)

∂θ2p(x; θ)dx (2.4)


in quanto la derivata seconda e una variabile casuale dipendente da x. Essorappresenta, come gia anticipato, il valor medio della curvatura della log-likelihood function.

Esempio 2.1 - Stimatore non polarizzato per livello dicorrente continua in rumore WGN

Riprendendo l’esempio 1.1, consideriamo l’osservazione

x[n] = A + w[n] n = 0, 1, ..., N − 1

dove w[n] rappresenta un rumore WGN con varianza σ2. Per determinare ilCRLB per la costante A, abbiamo, essendo i campioni scorrelati,

p(x; A) =N−1∏n=0

1√2πσ2

exp

[− 1

2σ2(x[n]− A)2

]

=1

(2πσ2)N2

exp

[− 1

2σ2

N−1∑n=0

(x[n]− A)2

]

Effettuando la derivata prima

∂ ln p(x; A)

∂A=

∂

∂A

[− ln

[(2πσ2)

N2

]− 1

2σ2

N−1∑n=0

(x[n]− A)2

]

=1

σ2

N−1∑n=0

(x[n]− A)

=N

σ2(x− A) (2.5)

dove x rappresenta il valor medio. Effettuando nuovamente la derivata

∂2 ln p(x; A)

∂A2= −N

σ2

e notando che la derivata seconda e costante, si ottiene dalla 2.2 il CRLB

var(A) ≥ σ2

N. (2.6)

Dimostriamo adesso che quando viene assunto il valore CRLB, si ha

var(θ) =1

I(θ)


dove

I(θ) = −E

[∂2 ln p(x; θ)

∂θ2

].

Dalle equazioni 2.2 e 2.3 si ottiene

var(θ) =1

−E[

∂2 ln p(x;θ)∂θ2

]

∂ ln p(x; θ)

∂θ= I(θ)(θ − θ)

∂2 ln p(x; θ)

∂θ2=

∂I(θ)

∂θ(θ − θ)− I(θ)

per cui il valore atteso cambiato di segno diventa

− E

[∂2 ln p(x; θ)

∂θ2

]= −∂I(θ)

∂θ

[E(θ)− θ

]+ I(θ)

= I(θ)

ed infine

var(θ) =1

I(θ). (2.7)

E’ importante sottolineare che il CRLB non e sempre soddisfatto.Uno stimatore che contemporaneamente sia non polarizzato e rispetti il cri-terio CRLB, viene denominato efficiente, in quanto utilizza efficientemente idati, sfruttando gran parte dell’informazione in essi contenuta.Intuitivamente, maggiore e l’informazione, minore sara il valore del bound inquanto minore e l’incertezza sulla misura. Si noti che non e detto che unostimatore MVU sia anche efficiente, in quanto pur essendo lo stimatore avarianza minima, puo non rispettare il criterio CRLB.La quantita I(θ) e detta Informazione di Fisher, ed e una quantita non-negativa in quanto puo essere dimostrata l’equivalenza

− E

[∂2 ln p(x; θ)

∂θ2

]= E

[(∂ ln p(x; θ)

∂θ

)2]

(2.8)

Un’altra importante caratteristica della I(θ) e la proprieta di additivita perosservazioni indipendenti. Questa proprieta porta direttamente alla conclu-sione che il CRLB per N osservazioni IID (Indipendenti e Identicamentedistribuite) e pari a 1/N volte il CRLB per una singola osservazione. Ovverosi ha:

I(θ) = Ni(θ)


dove

i(θ) = −E

[∂2 ln p(x[n]; θ)

∂θ2

](2.9)

e la Informazione di Fisher per un campione. Per osservazioni non indi-pendenti, ci aspettiamo che l’informazione sia minore della quantita Ni(θ);per osservazioni completamente dipendenti, in cui per esempio x[0] = x[1] =x[2] = ... = x[N − 1], avremo I(θ) = i(θ).

2.3 CRLB per un segnale immerso in rumore

bianco Gaussiano

Risulta molto utile ricavare le condizioni per il CRLB nel caso di rumorebianco Gaussiano, in quanto e molto comune nell’analisi di un segnale assu-mere un rumore con tali caratteristiche.Assumiamo di avere un’osservazione di un segnale deterministico, immersoin un segnale WGN e dipendente da un parametro θ

x[n] = s[n, θ] + w[n] n = 0, 1, ..., N − 1

dove abbiamo esplicitamente espresso la dipendenza da θ. La funzione diverosimiglianza e

p(x; θ) =1

(2πσ2)N2

exp

[− 1

2σ2

N−1∑n=0

(x[n]− s[n; θ])2

].

Effettuando la derivata prima della log-likelihood function, si ha:

∂ ln p(x; θ)

∂θ=

1

σ2

N−1∑n=0

(x[n]− s[n; θ])∂s[n; θ]

∂θ,

mentre effettuando la derivata seconda

∂2 ln p(x; θ)

∂θ2=

1

σ2

N−1∑n=0

[(x[n]− s[n; θ])

∂2s[n; θ]

∂θ2−

(∂s[n; θ]

∂θ

)2]

.

L’operatore di media produce

E

(∂2 ln p(x; θ)

∂θ2

)= − 1

σ2

N−1∑n=0

(∂s[n; θ]

∂θ

)2


per cui in definitiva abbiamo:

var(θ) ≥ σ2

N−1∑n=0

(∂s[n; θ]

∂θ

)2. (2.10)

L’espressione del bound dimostra la forte dipendenza del segnale da θ. Talemetodo produce stimatori molto accurati nel caso dell’analisi di segnali chevariano rapidamente al variare di θ (in quanto c’e una forte dipendenza dalladerivata seconda).I casi in cui si abbia

s[n; A] = A

s[n; ϕ] = A cos(2πf0n + ϕ)

risultano essere dei casi particolari del caso appena studiato. Vediamo unesempio.

Esempio 2.2 - Stima della frequenza di un segnale sinu-soidale

Sia il segnale in esame

s[n; f0] = A cos(2πf0n + φ) 0 < f0 <1

2

con A e φ note. Dall’equazione 2.10 si ottiene

var(f0) ≥ σ2

A2

N−1∑n=0

[2πn sin(2πf0n + φ)]2. (2.11)

In figura 2.1 e visibile l’andamento del CRLB in funzione della frequenza, conun SNR A2/σ2 = 1, lunghezza del record dati N = 10 e fase nulla. Si notiche se f0 → 0, il CRLB tende ad infinito, in quanto per valori di f0 prossimi azero, piccole variazioni di frequenza non alterano significativamente il segnale.

2.4 Trasformazione dei Parametri

Nella pratica accade frequentemente che il parametro che interessa stimaresia una funzione di altri parametri fondamentali. Ad esempio, piuttosto che


Figura 2.1: CRLB per l’esempio 2.2

essere interessati alla stima dell’ampiezza A di una sinusoide o del suo segno,potremmo essere interessati alla potenza del segnale, legata al quadrato diA, A2. Se conosciamo il CRLB per A, possiamo facilmente risalire al CRLBper A2 o piu in generale, per qualsiasi funzione di A.Il problema che vogliamo formalizzare e la stima di α = g(θ) invece che di θ.Si puo dimostrare che, in generale,

var(α) ≥

(∂g

∂θ

)2

−E

[∂2 ln p(x; θ)

∂θ2

] (2.12)

Per l’esempio appena riportato, si avrebbe α = g(A) = A2 e

var(A2) ≥ (2A)2

N/σ2=

4A2σ2

N.

Se α = g(θ) e una funzione non lineare e θ e uno stimatore efficiente, allorag(θ) non e uno stimatore efficiente, mentre esso conserva questa caratteristicanel caso in sui α = g(θ) sia una funzione lineare (o affine), come si puodimostrare facilmente dalla 2.12.

Tuttavia l’efficienza e mantenuta approssimativamente per trasformazio-ni non lineari con N molto grande. Infatti, quando N cresce, la PDF dig(θ) diventa piu concentrata intorno alla media θ = E(θ). Dunque possia-mo linearizzare g intorno alla media con una buona approssimazione, tanto


migliore quanto piu stretta e la PDF, ovvero quanto piu lungo e il recorddati:

g(θ) ≈ g(θ) +dg(θ)

dθ(θ − θ) (2.13)

Questo porta alle equazioni approssimate

E[g(θ)

]≈ g(θ) (2.14)

var[g(θ)

]≈

[dg(θ)

dθ

]2

var(θ) (2.15)

2.5 Stima di piu parametri - caso vettoriale

Andiamo ad estendere i risultati ottenuti nelle precedenti sezioni al caso diuna stima di piu parametri, o di un vettore di parametri θ = [θ1, θ2, ..., θp]

T .Assumeremo che θ non sia polarizzato.Si puo dimostrare che il CRLB, che ci permette di assegnare un bound per lavarianza dell’elemento i-esimo di θ, puo essere trovato come l’elemento [i, i]della matrice inversa

var(θi) ≥ [I−1(θ)]ii , (2.16)

dove I(θ) e la p× p matrice dell’informazione di Fisher p× p. Quest’ultimae definita come

[I(θ)]ij = −E

[∂2 ln p(x; θ)

∂θi∂θj

](2.17)

per i = 1, 2, ..., p, j = 1, 2, ..., p. Si noti che per p = 1, si ritorna al casoscalare. Per il calcolo della 2.17 viene utilizzato il valore vero di θ.

Esempio 2.3 - DC level in WGN

Generalizziamo l’esempio 2.1, supponendo che oltre ad A sia sconosciutaanche la varianza σ2. Il vettore, con p = 2, diventa dunque θ = [A σ2]. Lamatrice 2× 2, detta della Fisher Information, e

I(θ) =

−E

[∂2 ln p(x; θ)

∂A2

]−E

[∂2 ln p(x; θ)

∂A∂σ2

]

−E

[∂2 ln p(x; θ)

∂σ2∂A

]−E

[∂2 ln p(x; θ)

∂σ22

]

(2.18)

La matrice e simmetrica e si puo dimostrare che essa e semidefinita positiva.La funzione di verosimiglianza e

ln p(x; θ) = −N

2ln 2π − N

2ln σ2 − 1

2σ2

N−1∑n=0

(x[n]− A)2


Effettuando le derivate e l’operazione di media, la matrice diventa

I(θ) =

N

σ20

0N

2σ4

(2.19)

Nonostante non sia vero in generale, in questo caso la matrice e diagonalee puo essere facilmente invertita, trovando

var(A) ≥ σ2

N(2.20)

var(σ2) ≥ 2σ4

N. (2.21)

Notiamo che il CRLB per A e lo stesso dell’esempio 2.1, ma questo non puoessere generalizzato, come vediamo nell’esempio successivo.

Esempio 2.4 - Line fitting

Consideriamo il problema di line fitting (o regressione lineare) per una dataosservazione

x[n] = A + Bn + w[n] n = 0, 1, ..., N − 1

in cui w[n] e un rumore WGN; vogliamo determinare il CRLB per la coppiaθ = [A B]T . La matrice di Fisher e

I(θ) =

−E

[∂2 ln p(x; θ)

∂A2

]−E

[∂2 ln p(x; θ)

∂A∂B

]

−E

[∂2 ln p(x; θ)

∂B∂A

]−E

[∂2 ln p(x; θ)

∂B2

]

. (2.22)

La funzione di verosimiglianza e

p(x; θ) =1

(2πσ2)N2

exp

{− 1

2σ2

N−1∑n=0

(x[n]− A−Bn)2

}.

Effettuando le operazione di derivate e di media, si ottiene la matrice

I(θ) =1

σ2

N

N−1∑n=0

n

N−1∑n=0

n

N−1∑n=0

n2

=1

σ2

NN(N − 1)

2N(N − 1)

2

N(N − 1)(2N − 1)

2

. (2.23)


Invertendo la matrice,

I−1(θ) = σ2

2(2N − 1)

N(N + 1)− 6

(N(N + 1))

− 6

(N(N + 1))

12

(N(N2 − 1))

IL CRLB e allora

var(A) ≥ 2(2N − 1)σ2

N(N + 1)

var(B) ≥ 12σ2

N(N2 − 1).

Teorema 2 (Cramer Rao Lower Bound - Caso vettoriale) Si assumache la PDF di p(x; θ) soddisfi la condizione di regolarita

E

[∂ ln p(x; θ)

∂θ

]= 0 per ogni θ

dove la media e realizzata su p(x; θ) rispetto ad x. La matrice di covarianzadi un qualsiasi stimatore θ soddisfa l’equazione

Cˆθ− I−1(θ) ≥ 0

dove la disequazione indica che la matrice e semidefinita positiva. La Fisherinformation matrix I(θ) e data da

[I(θ)]ij = −E

[∂2 ln p(x; θ)

∂θi∂θj

](2.24)

dove le derivate sono effettuate rispetto al valore vero di θ e l’operazione dimedia e realizzata rispetto a p(x; θ). Inoltre lo stimatore raggiunge il boundC¯θ = I−1(θ) se e solo se

∂ ln p(x; θ)

∂θ= I(θ)(g(x)− θ) (2.25)

per qualche funzione p-dimensionale g, e qualche matrice I di dimensionip×p. Lo stimatore MVU cosı ottenuto e θ = g(x), e la sua matrice dicovarianza I−1(θ).

Capitolo 3

MVU per modelli lineari

La valutazione dello stimatore MVU e spesso agevolata dall’esistenza di mo-delli che formalizzano e semplificano una larga classe di problemi di stimanell’ambito dell’elaborazione dei segnali e che permettono di determinarepiuttosto facilmente lo stimatore: a questa classe appartengono i modelli li-neari, per i quali, con precisi artifici, e immediato trovare lo stimatore unavolta identificato il modello. Per questa classe la soluzione al problema dellavalutazione dello stimatore ottimo consiste nella costruzione del modello li-neare, dalle cui proprieta segue direttamente una soluzione.

Cominciamo dallo studio del teorema che permette di identificare e ca-ratterizzare la procedura di stima per un problema lineare.

3.1 Definizione e proprieta

Teorema 3 (MVU per un Modello Lineare) Se i dati osservati posso-no essere modellati come

x = Hθ + w (3.1)

dove x e il vettore di osservazione N × 1, H e la matrice di osservazioneN × p (N > p) e di rango p, θ e un vettore di parametri p× 1 che dev’esserestimato, w e un vettore di rumore N × 1 con una PDF N (0, σ2I), allora lostimatore MVU e

θ = (HTH)−1HTx (3.2)

e la matrice di covarianza di θ e

C ˆθ= σ2(HTH)−1 (3.3)

Inoltre per i modelli lineari lo stimatore MVU e efficiente in quanto la suavarianza soddisfa il CRLB.

25


Si noti che l’inversa di (HTH) esiste in quanto il rango di H e p, le colon-ne sono linearmente indipendenti. Inoltre non solo la media e la varianza,ma tutta la statistica di θ e completamente specificata, in quanto θ e unatrasformazione lineare di un vettore Gaussiano x e dunque

θ ∼ N (0, σ2(HTH)−1).

Passiamo all’applicazione del teorema, dapprima con un esempio che euna naturale estensione dell’esempio 2.4 al caso di una curva sperimentale,successivamente con l’analisi ingresso-uscita di un sistema di elaborazionenumerica.

Esempio 3.1 - Curve fitting

In molte situazioni sperimentali si cerca di determinare delle relazioni empi-riche fra due o piu variabili. Ad esempio, in figura 3.1, sono visibili i risultatidi un esperimento di misura di voltaggio agli istanti t = t0, t1, t2, ..., tN−1. I

Figura 3.1: Dati Sperimentali

valori della misura mostrano un andamento quadratico in funzione del tempo.Pertanto una ipotesi ragionevole per il modello dei dati puo essere

x(tn) = θ1 + θ2tn + θ3t2n + w(tn) n = 0, 1, 2, ..., N − 1.

Assumiamo che w(tn) siano variabili casuali gaussiane indipendenti e identi-camente distribuite (iid), con media nulla e varianza σ2 o che siano campionidi un rumore WGN. Il modello utilizzabile e

x = Hθ + w

dovex = [x(t0)x(t1)...x(tN−1)]

T ,


lo stimatore eθ = [θ1θ2θ3]

T

ed infine

H =

1 t0 t201 t1 t211 t2 t22...

......

1 tN−1 t2N−1

.

Piu in generale, se vogliamo adattare un polinomio di ordine (p − 1) aidati sperimentali,

x(tn) = θ1 + θ2tn + θ3t2n + . . . + θpt

p−1n + w(tn) n = 0, 1, 2, ..., N − 1,

la soluzione e ancoraθ = (HTH)−1HTx

x = Hθ + w

dovex = [x(t0)x(t1)...x(tN−1)]

T ,

H =

1 t0 · · · tp−10

1 t1 · · · tp−11

1 t2 · · · tp−12

......

. . ....

1 tN−1 · · · tp−1N−1

.

La matrice di osservazione in questo caso ha la forma di una matrice diVandermonde, e la risultante curva e

s(t) =

p∑i=1

θiti−1.

Esempio 3.2 - Sistema di identificazione

E’ molto frequente essere interessati ad identificare un modello di un sistemaattraverso l’analisi dei dati in ingresso e uscita. Osserviamo l’esempio in figu-ra 3.2(a), che implementa un filtro FIR: il sistema e pilotato da un ingressou[n], che serve per testare il sistema.Idealmente, all’uscita la sequenza

∑p−1k=0 h[k]u[n− k] ci permette di stimare

la risposta all’impulso del filtro. Nella pratica, tuttavia, l’uscita e corrotta


(a) Tapped delay line (TDL)

(b) Modello per dati corrotti da rumore

Figura 3.2: Esempi

dal rumore, pertanto assumere un modello come quello illustrato nella figura3.2(b) risulta piu corretto. Avremo

x[n] =

p−1∑

k=0

h[k]u[n− k] + w[n] n = 0, 1, ..., N − 1

dove si assume che u[n] = 0 per n < 0. In forma matriciale avremo

x =

u[0] 0 · · · 0u[1] u[0] · · · 0

......

. . ....

u[N − 1] u[N − 2] · · · u[N − p]

︸︷︷︸H

h[0]h[1]

...h[p− 1]

︸︷︷︸θ

+w (3.4)

Se assumiamo inoltre che w[n] sia WGN, allora la 3.4 e nella forma di unmodello lineare, per cui lo stimatore MVU per la risposta all’impulso e

θ = (HTH)−1HTx (3.5)

e la matrice di covarianza di θ e

C ˆθ= σ2(HTH)−1 (3.6)


Ci chiediamo adesso quale debba essere il segnale pilota u[n] affinche la va-rianza dello stimatore venga minimizzata: la risposta, come vedremo frabreve, e che tale segnale deve essere un rumore PRN (pseudorandom) perchetale condizione sia soddisfatta.La varianza di θi e

var(θi) = eTi Cθei

dove ei = [00...010...00]T dove ‘1’ occupa la i-esima posizione e C−1

θpuo essere

fattorizzato come DTD con D matrice invertibile p × p; possiamo usare ladisuguaglianza di Cauchy-Schwartz come segue.Siccome

(eTi DTDT−1

ei)2 = 1

possiamo imporre ξ1 = Dei e ξ2 = DT−1ei, che porta alla disequazione

(ξT1 ξ2)

2 ≤ ξT1 ξ1ξ

T2 ξ2.

Dall’uguaglianza ξT1 ξ2 = 1, si ha

1 ≤ (eTi DTDei)(e

Ti D−1DT−1

ei)

= (eTi C−1

θei)(e

Ti Cθei)

e in definitiva

var(θi) ≥ 1

(eTi C−1

θei)

=σ2

[HTH]ii.

Vale l’uguaglianza (e dunque l’MVU viene raggiunto) se e solo se ξ1 = cξ2

per una data costante c, ovvero se

Dei = ciDT−1

ei

o, equivalentemente, se

DTDei = ciei i = 1, 2, ..., p.

Siccome si ha

DTD = C−1

θ=

HTH

σ2

abbiamoHTH

σ2ei = ciei.


Si ottiene, dalla combinazioni delle precedenti espressioni in forma matriciale,che la condizione per il raggiungimento dell’MVU e esprimibile in funzionedella matrice H come

HTH = σ2

c1 0 · · · o0 c2 · · · 0...

.... . .

...0 0 · · · cp

. (3.7)

La condizione a cui deve soddisfare l’ingresso pilota u[n] per minimizzare lavarianza dello stimatore MVU e dunque quella espressa dalla 3.7. Siccome[H]ij = u[i− j],

[HTH] =N∑

n=1

u[n− i]u[n− j] i = 1, 2, ..., p j = 1, 2, ..., p (3.8)

e per N molto grande si ha

[HTH] ≈N−1−|i−j|∑

n=0

u[n]u[n+ | i− j | ] (3.9)

in cui puo essere riconosciuta una funzione di correlazione di una sequenzadeterministica u[n]. Inoltre con questa approssimazione HTH diventa unamatrice simmetrica di autocorrelazione di Toeplitz

HTH = N

ruu[0] ruu[1] · · · ruu[p− 1]ruu[1] ruu[0] · · · ruu[p− 2]

......

. . ....

ruu[p− 1] ruu[p− 2] · · · ruu[0]

(3.10)

con

ruu[k] =1

N

N−1−k∑n=0

u[n]u[n + k]

che puo essere vista come l’autocorrelazione di u[n]. Perche HTH sia diago-nale e necessario che

ruu[k] = 0 k 6= 0,

che e una condizione approssimativemente vera se usiamo una sequenza PRNin ingresso. In definitiva, sotto questa condizione, HTH = Nruu[0]I, edunque

var(h[i]) =1

Nruu[0]/σ2i = 0, 1, ..., p− 1. (3.11)


Scegliendo una sequenza PRN allora abbiamo ottenuto lo stimatore MVUcome

θ = (HTH)−1HTx (3.12)

in cui HTH = Nruu[0]I. Si ottiene

h[i] =1

Nruu[0]

N−1∑n=0

u[n− i]u[n] (3.13)

=

1

N

N−1−i∑n=0

u[n]x[n + i]

ruu[0](3.14)

in quanto u[n] = 0 se n < 0. Il numeratore nella 3.14 e la crosscorrelazionerux[i] fra le sequenze di ingresso e uscita, per cui, se usiamo una sequenzaPRN per identificare il sistema, lo stimatore MVU (per alti valori di N) e

h[i] =rux[i]

ruu[0]i = 0, 1, ..., p− 1 (3.15)

dove

rux[i] =1

N

N−1−i∑n=0

u[n]x[n + i]

e

ruu[0] =N−1∑n=0

u2[n]

Capitolo 4

General MVU Estimation

Finora abbiamo visto come i modelli lineari siano uno strumento utile per lavalutazione del CRLB e dunque dello stimatore MVU. Tuttavia, qualora nonesista uno stimatore efficiente o non ci si possa ricondurre a un modello linea-re, e ancora nel nostro interesse vagliare l’esistenza di uno stimatore MVU.Questo richiede la conoscenza del teorema Rao-Blackwell-Lehmann-Scheffe edel concetto di statistica sufficiente. Grazie a queste nozioni risulta possibilein molti casi determinare lo stimatore MVU da una semplice ispezione dellaPDF.

4.1 Statistiche sufficienti

Nel capitolo 1 abbiamo affrontato il problema della stima di un livello A dicorrente DC immerso in rumore WGN (si veda esempio 1.1): la media ditutti i campioni

A =1

N

N−1∑n=0

x[n]

e risultata essere proprio lo stimatore MVU, con varianza minima σ2/N .Tuttavia e ragionevole chiedersi quali campioni portino maggiore informa-zione al problema della stima, ovvero se esista un insieme o sottoinsiemedi campioni che sia sufficiente ai fini della stima. Nel caso del problemadell’esempio 1.1, avendo trovato come stimatore A, i seguenti set di dati ri-sulterebbero essere sufficienti, in quanto al fine della stima siamo interessati

32


solo ad una somma dei valori dei campioni:

S1 = {x[0], x[1], ..., x[N − 1]}S2 = {x[0] + x[1], ..., x[N − 1]}

S3 =

{N−1∑n=0

x[n]

}.

Il set di dati che contiene il minor numero di campioni pur essendo sufficien-te, viene chiamato set minimo o minima statistica sufficiente:in questo casodunque, S3 risulta dunque essere la minima statistica sufficiente.Per la stima di A infatti, una volta nota S3, non abbiamo piu bisogno diconoscere i valori dei singoli campioni, in quanto tutta l’informazione neces-saria e contenuta nella somma di essi.Per estendere questo concetto, si consideri la PDF dei dati

p(x; A) =N−1∏n=0

1√2πσ2

exp

[− 1

2σ2(x[n]− A)2

]

=1

(2πσ2)N2

exp

[− 1

2σ2

N−1∑n=0

(x[n]− A)2

](4.1)

in cui si assume che T (x) =∑N−1

n=0 x[n] = T0 sia stata osservata. Proprio laconoscenza del valore di questa statistica cambia la PDF dell’equazione 4.1in una PDF condizionata

p(x|N−1∑n=0

x[n] = T0; A),

che restituisce la PDF dopo l’osservazione della statistica sufficiente. Essendoquesta statistica sufficiente per la stima di A, questa PDF condizionata nondeve dipendere dal valore di A. Se cio accade, significa che potremmo dedurredai dati una maggiore informazione su A, oltre a quella gia fornita dallastatistica che riteniamo sufficiente, e cio vorrebbe dire che la statistica non eaffatto sufficiente. Si consideri come esempio la figura 4.1(a): se x = x0 perun x0 arbitrario, allora i valori di A prossimi ad A0 saranno piu probabili.Questo e in contrasto col fatto che

∑N−1n=0 x[n] e una statistica sufficiente.

In figura 4.1(b), invece, qualunque valore di A ha la stessa probabilita, epertanto l’osservazione di T (x) fornisce tutta l’informazione necessaria.


(a) Statistica T (x) non sufficiente (b) Statistica T (x) sufficiente

Figura 4.1: Definizione di statistica sufficiente

4.2 Come trovare la Statistica Sufficiente

Verra adesso affrontato il teorema della fattorizzazione di Neymann-Fisher,un utile strumento per trovare la statistica sufficiente, a cui seguiranno esempiapplicativi.

Teorema 4 (Neymann-Fisher Factorization) T (x)e una statistica suf-ficiente per la variabile θ se e solo se e possibile fattorizzare la PDF p(x; θ)come

p(x; θ) = g(T (x), θ)h(x) (4.2)

dove g e una funzione dipendente da x soltanto tramite la funzione T (x) e he una funzione dipendente solo da x.


Dimostriamo adesso che e possibile effettuare una fattorizzazione per l’equa-zione 4.1, assumendo che σ2 sia noto.A tal fine riscriviamo l’esponente della PDF come

N−1∑n=0

(x[n]− A)2 =N−1∑n=0

x2[n]− 2AN−1∑n=0

x[n] + NA2

in modo che la PDF si possa scrivere come

p(x; A) =1

(2πσ2)N2

exp

[− 1

2σ2

(NA2 − 2A

N−1∑n=0

x[n]

)]

︸︷︷︸g(T (x), A)

exp

[− 1

2σ2

N−1∑n=0

x2[n]

]

︸︷︷︸h(x)

.


T (x) =∑N−1

n=0 x[n] e una statistica sufficiente per A. Si noti inoltre che an-

che T ′(x) = 2∑N−1

n=0 x[n] e una statistica sufficiente per A, infatti qualsiasi

funzione biettiva di∑N−1

n=0 x[n] e una statistica sufficiente.

Esempio 4.2 - Power of WGN

Consideriamo ancora l’equazione data dalla 4.1, stavolta con σ2 sconosciutoe A = 0. Avremo

p(x; σ2) =1

(2πσ2)N2

exp

[− 1

2σ2

N−1∑n=0

x2[n]

]

︸︷︷︸g(T (x), A)

· [1]︸︷︷︸h(x)

.

Esempio 4.3 - Fase di una sinusoide

Si supponga di dover stimare la fase di una sinusoide in un rumore WGN

x[n] = A cos(2πf0n + φ) + w[n] n = 0, 1, ...N − 1.

L’ampiezza A della sinusoide e la frequenza f0 sono note, cosı come lavarianza σ2. La PDF e

p(x; φ) =1

(2πσ2)N2

exp

{− 1

2σ2

N−1∑n=0

[x[n]− A cos(2πf0n + φ)]2

}.

Si puo espandere l’esponente come

N−1∑n=0

x2[n]− 2AN−1∑n=0

x[n] cos(2πf0n + φ) +N−1∑n=0

A2 cos2(2πf0n + φ).

Sfruttando la proprieta

cos(α + β) = cos α cos β − sin α sin β

si ottieneN−1∑n=0

x2[n]− 2A

(N−1∑n=0

x[n] cos(2πf0n + φ)

)cos φ

+2A

(N−1∑n=0

x[n] sin(2πf0n + φ)

)sin φ +

N−1∑n=0

A2 cos2(2πf0n + φ).


Il problema non sembra che presenti una fattorizzazione del tipo del teo-rema di Neyman-Fisher. Dunque non esiste una sola statistica sufficiente.Possiamo pero fattorizzare la PDF come

p(x; φ) = g(T1(x), T2(x), φ) · h(x)

in cui

g(T1(x), T2(x), φ) =1

(2πσ2)N2

exp

{− 1

2σ2

[N−1∑n=0

A2 cos2(2πf0n + φ)−

−2AT1(x) cos φ + 2AT2(x) sin φ]} .

e

h(x) = exp

[1

(2πσ2)

N−1∑n=0

x2[n]

]

dove si e posto

T1(x) =N−1∑n=0

x[n] cos 2πf0n

e

T2(x) =N−1∑n=0

x[n] sin 2πf0n.

Con una lieve generalizzazione del teorema, possiamo affermare che T1(x)e T2(x) sono congiuntamente una statistica sufficiente per φ. Tuttavia, lageneralizzazione del teorema alle statistiche sufficienti congiunte non ha unagrande rilevanza a livello pratico, e pertanto verra tralasciata.

Teorema 5 (Rao-Blackwell-Lehmann-Scheffe) Se θ e uno stimatore nonpolarizzato di θ e T (x) e una statistica sufficiente per θ, allora θ = E(θ|T (x))e

• uno stimatore valido per θ (non dipendente da θ)

• non polarizzato

• uno stimatore con varianza minore o uguale a θ, per ogni θ.

Inoltre, se la statistica T (x) sufficiente e completa1, allora θ risulta essere lostimatore MVU.

1Una statistica si dice completa se esiste una sola funzione g della statistica che sia nonpolarizzata:

E[g(T (x))] = θ ∀θo equivalentemente, data v(T)=g(T)-h(T),

∫∞−∞ v(T )p(T ; θ)dT = 0 ∀θ, e soddisfatta solo

per v(T ) = 0, ∀T


Procedura per trovare lo stimatore MVU

• Applicare il teorema di Neymann-Fisher e trovare una singola statisticaT (x) sufficiente per θ

• Se T (x) e anche completo, si proceda, altrimenti ci si puo fermare

• Trovare una funzione g in modo tale che θ = g(T (x)) sia uno stimatorenon polarizzato

=⇒ θ e lo stimatore MVU.

In generale, l’alternativa all’ultimo passo e calcolare lo stimatore come θ =E(θ|T (x)), ma questo approccio e molto difficile da implementare, comedimostra il prossimo esempio.


Riprendiamo l’esempio 4.1. Sappiamo che A = x e lo stimatore MVU, inquanto raggiunge il limite CRLB. Tuttavia applicheremo il teorema RBLSche puo essere usato anche quando non esista uno stimatore efficiente, edunque quando non sia attuabile il metodo CRLB. Esistono due strade pertrovare A, lo stimatore MVU: entrambe sono basate sulla statistica sufficienteT (x) =

∑N−1n=0 x[n].

i. Trovare un qualsiasi stimatore non polarizzato di A, ad esempio A =x[0], e determinare A = E(A | T ). La media dev’essere effettuatarispetto a p(A | T )

ii. Trovare una funzione g tale che A = g(T ) sia uno stimatore nonpolarizzato di A

Per quanto riguarda il primo metodo possiamo assumere che lo stimatore siaA = x[0] e determinare A = E(x[0] | ∑N−1

n=0 x[n]).Per un vettore [x y]T , realizzazione di una distribuzione Gaussiana, aventecome media un vettore µ = [E(x)E(y)]T e matrice di covarianza

C =

[var(x) cov(x, y)

cov(x, y) var(y)

]


si puo mostrare che

E(x | y) =

∫ ∞

−∞xp(x | y)dx

=

∫ ∞

−∞xp(x | y)

p(y)dx

= E(x) +cov(x, y)

var(y)(y − E(y)). (4.3)

Applicando questo risultato a x = x[0] e y =∑N−1

n=0 x[n], si ha

[xy

]=

x[0]N−1∑n=0

x[n]

=

[1 0 0 . . . 01 1 1 . . . 1

]

︸︷︷︸L

x[0]x[1]

...x[N − 1]

.

Segue che la PDF di [x y]T e N (µ, C) in quanto rappresenta una trasforma-zione lineare di un vettore Gaussiano, in cui

µ = LE(x) = LA

11...1

=

[A

NA

]

C = σ2LLT = σ2

[1 11 N

]

Dalla 4.3, otteniamo lo stimatore MVU, dato da

A = E(x | y) = A +σ2

Nσ2

(N−1∑n=0

x[n]−NA

)

=1

N

N−1∑n=0

x[n].

Come anticipato, questo approccio, che richiede il calcolo della media condi-zionata, e generalmente difficile da trattare matematicamente.


4.3 Legame fra l’unicita di g e la completezza

della statistica

Il metodo prevede di trovare una qualche funzione g in modo che

A = g

(N−1∑n=0

x[n]

)

sia uno stimatore non polarizzato di A: possiamo scegliere g(x)∀T= x/N , che

ci porta ad ottenere

A =1

N

N−1∑n=0

x[n]

come stimatore MVU.Sappiamo che se g(x) e unica abbiamo E[g(

∑N−1n=0 x[n])] = A e T (x) risul-

ta completo. Tuttavia, supponiamo che esista un’altra funzione h tale cheE[h(

∑N−1n=0 x[n])] = A; questo vorrebbe dire che

E[g(T )− h(T )] = A− A = 0 ∀A.

Possiamo formalizzare il problema in maniera del tutto equivalente, impo-nendo che ∫ ∞

−∞v(T )p(T ; θ)dT = f(A) ∀θ

dove v(T ) = g(T )− h(T ) con E[g(∑N−1

n=0 x[n])] = E[h(∑N−1

n=0 x[n])] = A. Perdimostrare che T (x) e completo, e sufficiente allora dimostrare che

f(A)∀A= 0 ⇐⇒ v(T )

∀T= 0.

Dal fatto che T ∼ N (NA,Nσ2), imponiamo f(A)∀A= 0, ottenendo

∫ ∞

−∞v(T )

1√2πNσ2

exp

[− 1

2Nσ2(T −NA)2

]dT = 0 ∀A.

Imponendo τ = T/N e v′(τ) = v(Nτ),

∫ ∞

−∞v′(τ)

N√2πNσ2

exp

[− N

2σ2(A− τ)2

]dτ = 0 ∀A (4.4)

in cui si puo riconoscere la funzione di convoluzione fra v′(τ) con l’impulsoGaussiano w(τ) (si osservi la figura 4.2.)Affinche il risultato sia nullo, e necessario che v′(τ) sia identicamente nullo per


(a) Integrale nullo (b) Integrale diverso da zero

Figura 4.2: Condizione di completezza per una statistica sufficiente soddisfatta

tutti i valori di A. Un segnale e sempre nullo se e solo se la sua trasformatadi Fourier e identicamente nulla; possiamo utilizzare questa proprieta perfissare la condizione

V ′(f)W (f) = 0 ∀fdove V ′(f) = F [v′(τ)] e W (f) e la trasformata di Fourier dell’impulso Gaus-siano nella 4.4. D’altra parte W (f) e ancora Gaussiana e positiva per ognivalore di f , per cui la condizione e soddisfatta solo se V ′(f) = 0. In altreparole dev’essere v′(τ) = 0 ∀τ , ovvero g = h. Essendo la funzione g unica, lastatistica T (x) e completa e

A =1

M

N−1∑n=0

x[n]

e lo stimatore MVU.

Esempio 4.5 - Statistica sufficiente non-completa

Consideriamo il problema di stima di A per il dato

x[0] = A + w[0],

dove w[0] ∼ U [−12, 1

2]. Una statistica sufficiente e x[0], essendo l’unico da-

to disponibile, ed inoltre x[0] e uno stimatore non polarizzato. Possiamoaffermare che g(x[0]) = x[0] rappresenta un buon candidato per essere lostimatore MVU: dobbiamo pero verificare che questo stimatore soddisfi leproprieta di statistica sufficiente. Come nel paragrafo precedente, suppo-niamo che esista un’altra funzione h con la proprieta di non polarizzazione


h(x[0]) = A e proviamo a dimostrare che h = g.Sia dunque v(T ) = g(T ) − h(T ), ed esaminiamo le possibili soluzioni per vdell’equazione ∫ ∞

−∞v(T )p(x; A)dx = 0 ∀A.

Per questo problema tuttavia x = x[0] = T ,

∫ ∞

−∞v(T )p(T ; A)dT = 0 ∀A.

Inoltrep(T ; A)

p(T ; A) =

{1 A− 1

2≤ T ≤ A + 1

2

0 altrimenti(4.5)

e dunque la condizione si riduce a

∫ A+ 12

A− 12

v(T )dT = 0 ∀A.

Possiamo scegliere come funzione non nulla v(T ) = sin 2πT , che soddisfaquesta condizione come illustrato nella figura 4.3. La soluzione diviene

(a) Integrale nullo (b) Integrale diverso da zero

Figura 4.3: Condizione di completezza per una statistica sufficiente non soddisfatta

v(T ) = g(T )− h(T ) = sin 2πT

ovveroh(T ) = T − sin 2πT.


In definitiva, lo stimatore

A = x[0]− sin 2πx[0]

e basato su una statistica sufficiente per A e risulta essere non polarizzato.Ma abbiamo trovato almeno un altro stimatore non polarizzato con la stessastatistica sufficiente, pertanto possiamo concludere che la statistica non ecompleta; non e possibile applicare il teorema RBLS per trovare lo stimatoreMVU.

Esempio 4.6 - Media di un rumore bianco uniformemen-te distribuito

Supponiamo di osservare l’insieme di dati

x[n] = w[n] n = 0, 1, ..., N − 1,

in cui w[n] siano i campioni di un rumore IID (indipendente e identicamentedistribuito), con PDF U [0, β] con β > 0, dove U [a, b] indica una distribuzioneuniforme sull’intervallo [a, b]. Siamo interessati a trovare lo stimatore MVUper la media θ = β/2. Non possiamo usare il teorema CRLB in quanto la di-stribuzione non soddisfa le condizioni di regolarita. Sembrerebbe abbastanzanaturale scegliere per la media lo stimatore

θ =1

N

N−1∑n=0

x[n].

La varianza e

var(θ) =1

Nvar(x[n])

=β2

12N. (4.6)

Per determinare se questo sia lo stimatore MVU seguiremo una differenteprocedura. Definiamo il gradino unitario come

u(x) =

{1 per x > 00 per x < 0.

Possiamo riscrivere

p(x[n], θ) =1

β[u(x([n])− u(x[n]− β)],


in cui β = 2θ, e pertanto la PDF dei dati e

p(x, θ) =1

βn

N−1∏n=0

[u(x([n])− u(x[n]− β)].

Questa PDF e non nulla solo se 0 < x[n] < β per ogni x[n], pertanto si puoriscrivere

p(x, θ) =

1

βn0 < x[n] < β n = 0, 1, ..., N − 1

0 altrimenti.

In alternativa possiamo scrivere

p(x, θ) =

1

βnmax x[n] < β, min x[n] > 0

0 altrimenti

per cui

p(x, θ) =1

βnu(β −max x[n])

︸︷︷︸g(T (x), θ)

u(min x[n])︸︷︷︸h(x)

Per il teorema di Neyman-Fisher T (x) = max x[n] e una statistica sufficienteper θ, ed inoltre, omettendo la prova, possiamo dire che essa e anche unastatistica completa.Determiniamo il valore atteso di T = max x[n], calcolando dapprima la suafunzione di distribuzione cumulativa, tenendo presente che le variabili casualisono IID:

Pr{T ≤ ξ} = Pr{x[0] ≤ ξ, x[1] ≤ ξ, ..., x[N − 1] ≤ ξ}

=N−1∏n=0

Pr{x[n] ≤ ξ}

= Pr{x[n] ≤ ξ}N .

Effettuando la derivata per avere la PDF

pT (ξ) =d Pr{T ≤ ξ}

dξ

= N Pr{x[n] ≤ ξ}N−1d Pr{x[n] ≤ ξ}dξ

.


Mad Pr{x[n] ≤ ξ}

dξe la PDF di x[n]

px[n](ξ, θ) =

1

β0 < ξ < β

0 altrimenti

Integrando otteniamo

Pr{x[n] ≤ ξ} =

0 ξ < 0ξ

β0 < ξ < β

1 ξ > β

che in definitiva produce

pT (ξ) =

0 ξ < 0

N

(ξ

β

)N−11

β0 < ξ < β

0 ξ > β

L’operazione di media produce

E(T ) =

∫ ∞

−∞ξpT (ξ)dξ

=

∫ β

0

ξN

(ξ

β

)N−11

βdξ

=N

N + 1β

=2N

N + 1θ.

(4.7)

Affinche questo valore sia non polarizzato, dev’essere θ = [(N + 1)/2N ]T ,ovvero in definitiva lo stimatore MVU e

θ =N + 1

2Nmax

nx[n].

Contrariamente a quanto l’intuito puo suggerire, per un rumore uniforme-mente distribuito, la media campione non rappresenta lo stimatore MVU; siinvita a verificare quanto appena trovato con una simulazione in ambienteMATLAB.


Teorema 6 (Neymann-Fisher Factorization (Caso Vettoriale)) Se ri-sulta possibile fattorizzare la PDF p(x;θ) come

p(x; θ) = g(T(x),θ)h(x) (4.8)

dove g e una funzione dipendente da x attraverso T(x), statistica r × 1, eda θ, e h e una funzione dipendente solo da x, allora T(x) e una statisticasufficiente per θ.Viceversa, se T(x) e una statistica sufficiente per θ, allora la PDF puo esserefattorizzata come nella 4.8.

Teorema 7 (Rao-Blackwell-Lehmann-Scheffe (Caso Vettoriale)) Se θrisulta essere uno stimatore non polarizzato per θ e T(θ) e una statisticasufficiente r × 1 per θ, allora θ = E(θ | T(θ)) e

• uno stimatore valido per θ (non dipendente da θ)

• non polarizzato

• uno stimatore con varianza minore o uguale a θ (ogni elemento di θha varianza minore o uguale)

Inoltre, se la statistica sufficiente e anche completa, allora θ e lo stimatoreMVU.

Capitolo 5

Best Linear UnbiasedEstimator (BLUE)

E’ molto frequente nella pratica che i metodi descritti nei capitoli precedenti,basati sul teorema CRLB e sulle teorie della statistica sufficiente, non produ-cano dei risultati soddisfacenti per la mancanza di un modello appropriatoper la PDF. Inoltre, anche quando la PDF sia nota, l’applicazione di questimetodi non garantisce di trovare lo stimatore MVU. Risulta molto utile inquesti casi cercare uno stimatore subottimo: in questa operazione tuttavia,non possiamo essere sicuri della quantita di informazione che viene persa equanto vengano ridotte le performance, in quanto non conosciamo la mini-ma varianza dello stimatore MVU. Possiamo pero valutare la varianza dellostimatore subottimo, e, una volta specificati i vincoli del problema che sista cercando di risolvere, trovare uno stimatore che rispetti le specifiche delproblema.Come abbiamo gia avuto modo di constatare nel capitolo 3, un approccio chesemplifica il problema della stima consiste nel vincolare una classe di stimato-ri ad avere un comportamento lineare in funzione dei dati, e successivamenteestrarre da questa classe lo stimatore non polarizzato e con varianza minima.Questo stimatore, chiamato Best Linear Unbiased Estimator(BLUE), puoessere valutato semplicemente dalla conoscenza dei momenti di primo e se-condo ordine della PDF, proprieta che rende il BLUE molto piu maneggevolenei casi di interesse pratico.

5.1 Definizione e Vincoli per il BLUE

Si supponga di osservare un insieme di dati {x[0], x[1], ..., x[N − 1]} carat-terizzato da una PDF p(x, θ) dipendente da un parametro θ sconosciuto.

46


L’approccio BLUE vincola lo stimatore ad essere lineare, ovvero il problemaconsiste nel trovare delle costanti an tali che

θ =N−1∑n=0

anx[n]. (5.1)

Ovviamente da questa definizione non possiamo dedurre un solo stimato-re, in quanto la scelta delle costanti an puo determinare un gran numero distimatori per θ: fra questi, il BLUE e definito come lo stimatore a mini-ma varianza e non polarizzato. Dipendentemente dai vincoli del problemadella stima, il BLUE puo essere ottimo o subottimo: e ottimo solo quandolo stimatore MVU presenti un comportamento lineare. Chiariamo questoconcetto: nell’esempio 2.1, lo stimatore MVU era

θ = x =N−1∑n=0

1

Nx[n]

che e una combinazione lineare dei dati. Pertanto, se nella ricerca dellostimatore restringiamo la nostra attenzione solo agli stimatori lineari, nonintroduciamo nessuna perdita di informazione, in quanto l’MVU appartienea questa classe (si osservi la figura 5.1(a)). Nell’esempio 4.6, avevamo invece

(a) DC level in WGN; BLUE e lostimatore ottimo

(b) Media di un rumore uniforme;BLUE e subottimo

Figura 5.1: Ottimalita del BLUE


trovato uno stimatore definito come

θ =N + 1

2Nmax

nx[n]

che non e lineare. Pertanto, se costringiamo lo stimatore ad un vincololineare, il BLUE in questo caso diventa subottimo, come illustrato in figura5.1(b).Esistono anche problemi di stima per i quali il BLUE risulta completamenteinappropriato, come ad esempio nel problema di stima dell’esempio 4.2, peril quale lo stimatore MVU e

σ2 =1

N

N−1∑n=0

x2[n]

che e evidentemente non lineare. Se vincoliamo lo stimatore a ad esserelineare, in modo che

σ2 =N−1∑n=0

anx[n],

il valore atteso dello stimatore diventa

E(σ2) =N−1∑n=0

anE(x[n]) = 0,

in quanto E(x[n]) = 0 per tutti i valori di n. Non possiamo trovare unostimatore lineare che sia non polarizzato. In questo caso tuttavia, nono-stante il BLUE non sia adatto al tipo di problema affrontato, utilizzandouna trasformazione sui dati del tipo y[n] = x2[n], si produce uno stimatoreefficace.

Come gia anticipato, dopo aver vincolato θ ad essere lineare e non pola-rizzato, si calcolano i coefficienti an in modo tale da minimizzare la varianza.Per la proprieta di non polarizzazione, il vincolo puo essere formalizzatomatematicamente come

E(θ) =N−1∑n=0

anE(x[n]) = θ (5.2)

che e un vincolo sul primo ordine della distribuzione di θ; cosicche il vincolosulla varianza - momento di secondo ordine - risulta essere

var(θ) = E

(N−1∑n=0

anx[n]− E

(N−1∑n=0

anx[n]

))2 . (5.3)


Usando la 5.2 e ponendo a = [a0a1...aN−1]T possiamo riscrivere

var(θ) = E[(aTx− aT E(x))2

]

= E[(aT (x− E(x)))2

]

= E[aT (x− E(x))(x− E(x))Ta

]

= aTCa, (5.4)

dove C e la matrice di covarianza. Il vettore a dei pesi puo essere trovato mi-nimizzando la 5.4 rispetto al vincolo espresso dalla 5.2. Perche sia soddisfattoil vincolo di non polarizzazione, E(x[n]) dev’essere lineare in θ,

E(x[n]) = s[n]θ, (5.5)

dove i coefficienti s[n] devono essere noti, altrimenti il vincolo non puo esseresoddisfatto.Si noti che se scriviamo x[n] come

x[n] = E(x[n]) + [x[n]− E(x[n])]

e osservando che x[n]−E(x[n]) puo essere interpretato come una dispersionedei dati, e quindi come un rumore w[n], abbiamo

x[n] = θs[n] + w[n].

La 5.5 indica che il BLUE e utile nella stima dell’ampiezza di un segnale notoimmerso in un rumore.Partendo dalla stessa assunzione espressa dalla 5.5, andiamo adesso a rias-sumere i passi fondamentali del problema della stima. Per trovare il BLUE,dobbiamo minimizzare la varianza

var(θ) = aTCa,

che grazie alle 5.2 e 5.5, puo essere riscritta come

N−1∑n=0

anE(x[n]) = θ

N−1∑n=0

ans[n]θ = θ

N−1∑n=0

ans[n] = 1


ovveroaT s = 1

in cui si e posto s = [s[0], s[1], ..., s[N − 1]]T .Per minimizzare la quantita var(θ) = aTCa rispetto al vincolo aT s = 1 usiamoil metodo dei moltiplicatori di Lagrange, con la funzione Lagrangiana J

J = aTCa + λ(aT s− 1).

Il gradiente della funzione Lagrangiana rispetto ad a e

∂J∂a

= 2Ca + λs.

Uguagliando a zero e risolvendo, si ottiene

a = −λ

2C−1s.

La costante λ, moltiplicatore di Lagrange, dev’essere trovata imponendol’equazione del vincolo; si ha

aT s = −λ

2sTC−1s = 1,

da cui si deduce

−λ

2=

1

sTC−1s

e pertanto le due condizioni - gradiente nullo e vincolo di non polarizzazione- sono soddisfatte per

aopt =C−1s

sTC−1s. (5.6)

Pertanto utilizzando il risultato 5.6, si ottiene per il BLUE l’espressione

θ = aTx =sTC−1x

sTC−1s. (5.7)

L’espressione della varianza per lo stimatore BLUE che si trova in corrispon-denza di aopt e

var(θ) = aToptCaopt

=sTC−1CC−1s

(sTC−1s)2

=1

sTC−1s(5.8)


Si noti che, essendo E(x) = θs, il BLUE e non polarizzato:

E(θ) =sTC−1E(x)

sTC−1s

=sTC−1θs

sTC−1s= θ. (5.9)

Risulta piu chiaro a questo punto il motivo per cui il BLUE necessita solodella conoscenza delle stime di primo e secondo ordine. Piu precisamente, enecessaria la conoscenza di

1. s, che e equivalente, a meno di una costante, alla media;

2. C, matrice di covarianza.

Esempio 5.1 - DC level in white noise

Supponiamo di osservare un insieme di dati

x[n] = A + w[n] n = 0, 1, ..., N − 1

in cui w[n] e un rumore bianco con varianza σ2 e la PDF e sconosciuta (o nondefinita); vogliamo stimare A. Siccome il rumore w[n] non e necessariamenteGaussiano, i campioni possono essere statisticamente dipendenti anche sescorrelati 1.Siccome E(x[n]) = A, abbiamo dalla 5.5 che s[n] = 1, ovvero s = [1, 1, ..., 1].Pertanto, dalla 5.7, il BLUE e:

A =sT 1

σ2Ix

sT1

σ2Is

=1

N

N−1∑n=0

x[n] = x (5.10)

ed ha varianza minima, dalla 5.8

var(A) =1

sT1

σ2s

=σ2

N(5.11)

1Si ricorda che, fatta eccezione per le variabili Gaussiane, l’indipendenza implica lascorrelatezza, ma non viceversa.


Dunque la media dei campioni risulta essere proprio il BLUE, e qualora ilrumore si possa assumere con distribuzione Gaussiana, il BLUE e anche lostimatore MVU.

Esempio 5.2 - DC level in uncorrelated zero-mean noise

Sia adesso w[n] un rumore non correlato e a media nulla, con varianza σ2n, e

riprendiamo l’esempio 5.1. Sia ancora s = [1, 1, ..., 1], e

A =sTC−1x

sTC−1s

var(A) =1

sTC−1s.

La matrice di covarianza e

C =

σ20 0 . . . 00 σ2

1 . . . 0...

.... . .

...0 0 . . . σ2

N−1

,

e la sua inversa e

C−1 =

1

σ20

0 . . . 0

01

σ21

. . . 0

......

. . ....

0 0 . . .1

σ2N−1

.

Dunque

A =

N−1∑n=0

x[n]

σ2n

N−1∑n=0

1

σ2n

(5.12)

var(A) =1

N−1∑n=0

1

σ2n

.

Il denominatore nella 5.12 e il fattore di scala necessario affinche sia soddi-sfatta la condizione di non polarizzazione. Notiamo inoltre che in generale


possiamo dire che la presenza della matrice C−1 nello stimatore BLUE ha l’ef-fetto di una operazione di pre-whitening sui dati che vengono successivamentemediati.

5.2 BLUE - Caso Vettoriale

Se il parametro da stimare e un vettore p×1, allora la condizione da soddisfarediventa

θi =N−1∑n=0

ainx[n] i = 1, 2, ..., p (5.13)

dove gli ain rappresentano la matrice dei pesi A.In forma matriciale dunque

θ = Ax,

in cui A e una matrice p×N . Affinche θ sia non polarizzato, dev’essere

E(θi) =N−1∑n=0

ainE(x[n]) = θi i = 1, 2, ..., p (5.14)

o in forma matricialeE(θ) = AE(x) = θ. (5.15)

Bisogna ricordare che uno stimatore lineare puo essere appropriato solo perproblemi nei quali possa essere soddisfatto il vincolo di non polarizzazione.Dalla 5.15 risulta che e necessario avere

E(x) = Hθ, (5.16)

con H matrice N × p nota, che generalizza il ruolo del vettore s incontratonel caso scalare.Sostituendo la 5.16 nella 5.15, otteniamo il vincolo per la condizione di nonpolarizzazione

AH = I. (5.17)

Se definiamo il vettore ai = [ai0ai1...ai(N−1)], in modo da avere θi = aTi x , il

vincolo puo essere riscritto per ogni ai, con

A =

aT1

aT2...

aTp


eH = [ h1 h2 . . . hp ] .

Introducendo queste definizioni infatti, il vincolo descritto dalla 5.17 si riducealla condizione

aTi hi = δij i = 1, 2, ..., p ; j = 1, 2, ...p. (5.18)

Per analogia con il caso scalare, la varianza e

var(θi) = aTi Cai (5.19)

Il BLUE dev’essere trovato proprio minimizzando la 5.19 rispetto al vincolo5.18, e ripetendo la minimizzazione per ogni indice i. Si puo dimostrare chela minimizzazione si ha per

θ = (HTC−1H)−1HTC−1x (5.20)

e la matrice di covarianza

Cθ = (HTC−1H)−1. (5.21)

La forma del BLUE appena descritta estende i risultati trovati nel Capitolo 3e negli esempi 5.1 e 5.2. Possiamo affermare che quando i dati siano realmenteGaussiani, allora il BLUE corrisponde allo stimatore MVU: in questo casoinfatti, vincolare lo stimatore ad essere lineare non porta al compromessodi uno stimatore subottimo, in quanto l’MVU appartiene proprio all’insiemedegli stimatori lineari.

Teorema 8 (Gauss-Markov) Se i dati osservati sono nella forma generaledi un modello lineare

x = Hθ + w

in cui H sia una matrice N × p nota, θ e il vettore p × 1 dei parametri dastimare, e w e un vettore di rumore N × 1, con media nulla e covarianza C,allora il BLUE di θ e

θ = (HTC−1H)−1HTC−1x (5.22)

e la minima varianza di θi e

var(θi) =[(HTC−1H)−1

]ii

. (5.23)

Inoltre, la matrice di covarianza di θ e

Cθ = (HTC−1H)−1. (5.24)

Capitolo 6

Stima a massimaverosimiglianza (MLE)

Il metodo di stima MLE (Maximum Likelihood Estimation) si pone comemetodo alternativo all’MVU, e viene in aiuto qualora l’MVU non esista ocomunque non possa essere trovato. Basato sul principio della massima ve-rosimiglianza, il metodo MLE e sicuramente il piu utilizzato nella praticagrazie alla sua semplice implementazione anche nel caso di complicati pro-blemi di stima e grazie anche alle buone prestazioni su grandi quantita didati.

6.1 MLE - Stimatori Consistenti (caso scala-

re)

L’MLE per un parametro scalare θ e il valore di θ che massimizza la funzionedi verosimiglianza p(x, θ) di tutti i valori possibili di θ per x fissato.Prima di descrivere i metodi di stima MLE, passiamo ad illustrare un esempioper comprendere l’utilizzo e l’attuabilita degli stimatori a massima verosimi-glianza.

Esempio 6.1 - DC level in white Gaussian Noise

Supponiamo di osservare i dati

x[n] = A + w[n], n = 0, 1, . . . , N − 1

dove si assume che A sia un livello di segnale sconosciuto (ma positivo, A >0), e w[n] sia un rumore WGN con varianza A sconosciuta. In questo caso

55


dunque, nonostante il problema sia simile ai precedenti, si ha che non solo Ae il valore di media, ma anche della varianza dei dati. Nel cercare lo stimatoreMVU di A prima cercheremo di determinare il CRLB (si faccia riferimentoal capitolo 3). La PDF

p(x; A) =N−1∏n=0

1√2πA

exp

[− 1

2A(x[n]− A)2

]

=1

(2πA)N2

exp

[− 1

2A

N−1∑n=0

(x[n]− A)2

](6.1)

Effettuando la derivata della log-likelihood function(considerando il logaritmodella PDF come una funzione di A), otteniamo

∂ ln p(x; A)

∂A= − N

2A+

1

A

N−1∑n=0

(x[n]− A) +1

2A2

N−1∑n=0

(x[n]− A)2

?= I(A)(A− A). (6.2)

Gia da una prima occhiata alla 6.2 non appare affatto facile esprimere laderivata nella forma richiesta dal teorema di Cramer-Rao; in effetti si puodimostrare che cio non e possibile, e di conseguenza che non esiste uno sti-matore efficiente. Si puo inoltre dimostrare che il CRLB per questo tipo diproblema e

var(A) ≥ A2

N(A + 12)

(6.3)

e che procedendo nel calcolo di uno stimatore efficiente mediante il teoremadella fattorizzazione di Neyman-Fisher, non si puo ne ottenere una statisticasufficiente che sia non polarizzata, ne ottenere una media condizionata ma-tematicamente maneggevole.Esaurite tutte le possibilita di utilizzare un approccio che porti ad uno sti-matore ottimo o subottimo, dobbiamo affidarci ad un metodo da cui si pos-sa derivare uno stimatore approssimativamente ottimo, ovvero che diventiefficiente per N che tende ad infinito. Questa condizione comporta che

se N →∞ allora=⇒

{E(A) → A

var(A) → CRLB, (6.4)

dove il CRLB e dato dalla 6.3. Uno stimatore A che soddisfi la condizione sul-la media nella 6.4 e detto asintoticamente non polarizzato. Se invece soddisfaentrambe le condizioni della 6.4, viene detto asintoticamente efficiente. Per


un insieme di dati finiti tuttavia non possiamo utilizzare dei metodi efficacial fine della valutazione delle caratteristiche di ottimalita.

Ritornando all’esempio, possiamo proporre lo stimatore

A = −1

2+

√√√√ 1

N

N−1∑n=0

x2[n] +1

4. (6.5)

che risulta polarizzato in quanto

E(A) = E

−1

2+

√√√√ 1

N

N−1∑n=0

x2[n] +1

4

6= −1

2+

√√√√E

(1

N

N−1∑n=0

x2[n]

)+

1

4∀A

= −1

2+

√A + A2 +

1

4= A

Tuttavia lo stimatore e ragionevole, in quanto, per la legge dei grandi numeri,quando N tende ad infinito

1

N

N−1∑n=0

x2[n] → E(x2[n]) = A + A2 (6.6)

pertanto dalla 6.5A → A.

Lo stimatore nella 6.5 e detto stimatore consistente. Si puo inoltre dimostrareche asintoticamente questo stimatore rispetta le condizioni espresse dalla 6.4.

Lo stimatore 6.5 proposto e proprio l’MLE; nei prossimi paragrafi pas-siamo alla descrizione delle metodologie analitiche utilizzate per ottenerel’espressione dello stimatore MLE.

6.2 Trovare lo stimatore MLE

Come gia anticipato nel paragrafo precedente l’MLE per un parametro θ eil valore di θ che massimizza la funzione di verosimiglianza p(x, θ) di tuttii valori possibili di θ. Essendo p(x, θ) una funzione di x, l’operazione di


Figura 6.1: MLE:p(x, θ) per x = x0

massimizzazione produce un θ che e ancora una funzione di x.In figura 6.1, la PDF e valutata per x = x0; il valore di p(x = x0, θ)dx per

ogni valore di θ indica la probabilita di osservare x in una regione di volumedx in RN , centrata attorno ad x0 per un dato valore di θ. Osservando lafigura, si puo capire che sarebbe sbagliato scegliere θ = θ1: infatti, per θ = θ1,la probabilita di osservare x = x0 e davvero piccola. E’ invece altamente piuprobabile che θ = θ2 sia il valore che produce la piu alta probabilita diosservare x = x0.

Riprendiamo l’esempio precedente. Riscriviamo la PDF della 6.1

p(x; A) =1

(2πA)N2

exp

[− 1

2A

N−1∑n=0

(x[n]− A)2

]. (6.7)

Considerando questa espressione come una funzione di A, possiamo vederlacome la log-likelihood function. Derivando rispetto ad A otteniamo

∂ ln p(x; A)

∂A= − N

2A+

1

A

N−1∑n=0

(x[n]− A) +1

2A2

N−1∑n=0

(x[n]− A)2

che uguagliata a zero produce

A2 + A− 1

N

N−1∑n=0

x2[n] = 0.

Le due soluzioni sono

A = −1

2±

√√√√ 1

N

N−1∑n=0

x2[n] +1

4,

e dovendo essere A > 0, scegliamo quella positiva.


Teorema 9 (Proprieta asintotiche dell’MLE) Se la PDF p(x; θ) dei da-ti x soddisfa la condizione di regolarita1, allora l’MLE del parametro scono-sciuto e distribuita asintoticamente (per N molto grande), secondo la

θa∼ N (θ, I−1(θ)) (6.8)

dove I(θ) e la Fisher Information valutata per il valore vero del parametrosconosciuto.

Il simboloa∼ significa proprio ‘distribuita secondo’. Questo risultato e molto

importante ed e uno strumento che offre la possibilita di valutare ed affermarel’ottimalita dell’MLE.

Esempio 6.2 - MLE per la stima di una fase di unasinusoide

Riprendiamo l’esempio 4.3, nel quale avevamo cercato di stimare la fase diuna sinusoide immersa in un rumore w[n] con caratteristiche WGN e varianzaσ2:

x[n] = A cos(2πf0n + φ) + w[n] n = 0, 1, ...N − 1

e siano note l’ampiezza A e la frequenza f0.Avevamo visto nel capitolo 4 che esiste una statistica sufficiente congiunta

T1(x) =N−1∑n=0

x[n] cos 2πf0n

e

T2(x) =N−1∑n=0

x[n] sin 2πf0n.

L’MLE dev’essere trovato massimizzando la

p(x; φ) =1

(2πσ2)N2

exp

{− 1

2σ2

N−1∑n=0

[x[n]− A cos(2πf0n + φ)]2

}

o, equivalentemente, minimizzando

J(φ) =N−1∑n=0

(x[n]− A cos(2πf0n + φ))2.

1La condizione di regolarita e E

(∂ ln p(x;A)

∂A

)= 0


Differenziando rispetto a φ

∂J(φ)

∂φ= 2

N−1∑n=0

(x[n]− A cos(2πf0n + φ))A sin(2πf0n + φ)

che posto uguale a zero, restituisce

N−1∑n=0

x[n] sin(2πf0n + φ) = A

N−1∑n=0

cos(2πf0n + φ) sin(2πf0n + φ).

Si puo dimostrare che quando f0 non sia troppo vicino a 0 e 1/2, il rapportofra la parte a destra dell’espressione ed N e nulla:

1

N

N−1∑n=0

cos(2πf0n + φ) sin(2πf0n + φ) ≈ 0.

Pertanto cio che rimane e:

N−1∑n=0

x[n] sin(2πf0n + φ) = 0,

che, sfruttando la propieta

sin(α + β) = sin α cos β + cos α sin β,

ci fornisce l’MLE della fase

φ = − arctan

N−1∑n=0

x[n] sin(2πf0n)

N−1∑n=0

x[n] cos(2πf0n)

(6.9)


Figura 6.2: Valori numerici per l’esempio 6.2; η = (A2/2)/σ2 e il rapporto segnale rumore

Figura 6.3: Media in funzione dell’SNR per l’esempio 6.2;N = 80

Figura 6.4: Varianza in funzione dell’SNR per l’esempio 6.2;N = 80


Teorema 10 (Proprieta di invarianza dell’MLE) L’MLE di un para-metro α = g(θ), dove la PDF p(x; θ) e parametrizzata su θ, e dato da

α = g(θ),

in cui θ e l’MLE per θ. L’MLE di θ e ottenuto massimizzando p(x; θ). Se gnon e biettiva, allora α massimizza la funzione modificata

p(x; α) = max{θ:α=g(θ)}

p(x; θ).

6.3 Valutazione numerica dell’MLE

Il metodo MLE permette di determinare lo stimatore anche per via numeri-ca: cio significa che nella pratica possiamo determinare lo stimatore in ognicaso per un certo set di dati, in quanto possiamo ottenere l’MLE medianteuna ricerca numerica del massimo della funzione di verosimiglianza. Ovvia-mente se la funzione ammette solo valori di θ compresi in un intervallo [a, b],allora possiamo massimizzare p(x; θ) solo per tale intervallo, semplificando dimolto il calcolo. Il modo piu sicuro per fare questo e ovviamente una ricercaesaustiva, in cui venga testata una griglia di valori generati: se i valori di θsono abbastanza vicini, possiamo essere certi di trovare l’MLE per un datoinsieme di dati. Se tuttavia il range di θ non e un intervallo finito, comepuo accadere nella stima della varianza di un processo di rumore, allora unaricerca a griglia non e realizzabile: in casi simili, siamo obbligati a svilupparedelle procedure iterative.Fra le tecniche iterative piu utilizzate, possiamo elencare

• il metodo di Newton-Raphson;

• il metodo di scoring ;

• il metodo di minimizzazione del valore atteso (Caso Vettoriale),

che in generale convergono ad una soluzione se le condizioni iniziali dell’al-goritmo sono abbastanza vicine al valore vero.In ogni caso, la funzione da massimizzare non e nota a priori: la funzionedi verosimiglianza cambia per ogni set di dati, pertanto e richiesto che lamassimizzazione avvenga per una funzione casuale.Procediamo con una descrizione pratica dei metodi, applicandoli ad un realeproblema di stima.


Esempio 6.3 - Metodo di Newton-Raphson - Metodo diScoring

Consideriamo l’insieme di dati

x[n] = rn + w[n] n = 0, 1, . . . N − 1

dove w[n] e un rumore WGN con varianza σ2.Vogliamo stimare r, e sia r > 0.L’MLE per r e il valore che massimizza la funzione di verosimiglianza

p(x; r) =1

(2πσ2)N2

exp

[− 1

2σ2

N−1∑n=0

(x[n]− rn)2

]

o equivalentemente che minimizza

J(r) =N−1∑n=0

(x[n]− rn)2.

Differenziando J(r) e ponendo la derivata uguale a zero, si ottiene

N−1∑n=0

(x[n]− rn)nrn−1 = 0 (6.10)

che e una equazione non lineare di r e non puo essere risolta direttamente.Andiamo dunque ad applicare il metodo di Newton-Raphson.

Il metodo iterativo mira a massimizzare la log-likelihood function trovandogli zeri della derivata

∂ ln p(x; θ)

∂θ= 0. (6.11)

Sia

g(θ) =∂ ln p(x; θ)

∂θ,

e assumiamo di avere una ipotesi iniziale per la 6.10, detta θ0. Allora, se g(θ)e approssimativamente lineare nell’intorno di θ0, possiamo approssimare lafunzione con

g(θ) ≈ g(θ0) +dg(θ)

dθ

∣∣∣∣θ=θ0

(θ − θ0) (6.12)

come mostrato in figura 6.5.Risolviamo l’equazione 6.12 rispetto a θ, e chiamiamo θ1 lo zero trovato:

θ1 = θ0 − g(θ0)

dg(θ)

dθ

∣∣∣∣θ=θ0

.


Figura 6.5: Metodo di Newton-Raphson per determinare gli zeri di una funzione

Linearizziamo ancora la funzione g, ma stavolta nell’intorno di θ1, e ripetiamol’operazione eseguita al primo passo per trovare ancora lo zero della funzionelinearizzata.In generale possiamo dire che l’iterazione al passo k+1 produce il valore θk+1

dalla conoscenza del valore θk, usando l’espressione

θk+1 = θk − g(θk)

dg(θ)

dθ

∣∣∣∣θ=θk

.

Come illustrato nella figura 6.5, il metodo, se ben condizionato, condurraallo zero della funzione g(θ), convergendo alla condizione θk = θk+1.Possiamo anche riscrivere la condizione come

θk+1 = θk −[∂2 ln p(x; θ)

∂θ2

]−1∂ ln p(x; θ)

∂θ

∣∣∣∣θ=θk

. (6.13)

Nell’utilizzo di questa procedura si devono pero tener conto due aspetti:

• l’iterazione puo anche non convergere: questo e particolarmente eviden-te se si osserva che al denominatore della 6.13, e presente la derivataseconda della funzione di verosimiglianza che se presenta valori par-ticolarmente piccoli puo provocare ampie fluttuazioni da iterazione aiterazione;

• anche quando la procedura converga, il punto trovato potrebbe esserenon un massimo assoluto, ma un massimo locale o addirittura un mini-mo locale. Per evitare questo tipo di problema tipicamente si scelgonopiu punti da cui far partire l’algoritmo.


Applicando la 6.13 all’esempio, si ottiene un metodo iterativo descritto da

rk+1 = rk −

N−1∑n=0

(x[n]− rnk )nrn−1

k

N−1∑n=0

nrn−2k [(n− 1)x[n]− (2n− 1)rn

k ]

(6.14)

Una seconda comune tecnica iterativa, e detta metodo di scoring, cheutilizza la Fisher Information (si veda il capitolo 2) come

∂2 ln p(x; θ)

∂θ2

∣∣∣∣θ=θk

≈ −I(θk) (6.15)

valida per campioni IID e alti valori di N . Si ottiene con questa approssima-zione

∂2 ln p(x; θ)

∂θ2=

N−1∑n=0

∂2 ln p(x[n]; θ)

∂θ2

= N1

N

N−1∑n=0

∂2 ln p(x[n]; θ)

∂θ2

≈ NE

[∂2 ln p(x[n]; θ)

∂θ2

]

= −Ni(θ)

= −I(θ)

per la legge dei grandi numeri. La sostituzione della derivata seconda con ilsuo valore atteso migliora la stabilita dell’iterazione. Il metodo diventa

θk+1 = θk + I−1(θ)∂ ln p(x; θ)

∂θ

∣∣∣∣θ=θk

(6.16)

che e proprio la formulazione del metodo di scoring, che ha gli stessi problemidi convergenza del metodo di Newton-Raphson. Applicato al problema chestavamo considerando, abbiamo

I(θ) =1

σ2

N−1∑n=0

n2r2n−2

rk+1 = rk −

N−1∑n=0

(x[n]− rnk )nrn−1

k

N−1∑n=0

n2r2n−2k

(6.17)


L’ultimo algoritmo, detto EM, si basa sul fatto che esistono insiemi didati che permettono una stima dell’MLE molto piu semplice rispetto all’in-sieme di dati originale: piu precisamente, si suppone che esista un insiemedi dati completo, e che esista una trasformazione invertibile g da un insiemeincompleto di dati ad uno completo:

x = g(y1, y2, . . . , yM) = g(y).

L’obiettivo in questo caso e trovare l’MLE di θ non massimizzando la log-likelihood function [ln px(x,θ)], ma la [ln py(y,θ)] o equivalentemente, la

Ex|y[ln py(y, θ)] =

∫[ln py(y, θ)][ln py(y|x, θ)]dy (6.18)

utilizzando nell’integrale l’informazione iniziale θk.

Capitolo 7

Metodo dei minimi quadrati(Least Squares)

Nei capitoli precedenti abbiamo cercato di trovare degli stimatori ottimi oquasi-ottimi (per grandi moli di dati), scegliendo all’interno di un insieme distimatori non polarizzati quello che presenta la varianza minima. Andiamoadesso ad investigare una classe di stimatori che in generale non possiedeproprieta di ottimalita, ma risulta efficace in molti casi di interesse pratico:e il metodo dei minimi quadrati, che risale al 1795, anno in cui Gauss applicoquesto metodo nello studio dei moti planetari.La caratteristica fondamentale di questo metodo e il fatto che utilizza soloun modello per il segnale in esame senza la necessita di fare ipotesi probabi-listiche sui dati: per questa caratteristica il metodo risulta applicabile in unalarghissima classe di applicazioni. Tuttavia non ci sono criteri per stabilirel’ottimalita del metodo, ed inoltre le performance statistiche non possonoessere valutate per la mancanza di specifiche assunzioni sulle proprieta sta-tistiche dei dati.Il metodo risulta concettualmente di semplice implementazione, in quantoconsiste nel minimizzare una quantita di errore ai minimi quadrati.

7.1 Generalita sul problema dei minimi qua-

drati

I metodi descritti in precedenza utilizzano la varianza dello stimatore comeuna misura della qualita della stima effettuata, con l’obiettivo comune diminimizzare la differenza fra la stima e il valore vero (ovviamente in media).Nel metodo dei minimi quadrati, LS (Least Squares), si cerca di mini-

67


mizzare la differenza quadratica fra i dati x[n] acquisiti e il modello di segnalesenza rumore. Si osservi la figura 7.1. Il segnale s[n], deterministico, viene

(a) Modello di segnale

(b) Errore ai minimi quadrati

Figura 7.1: Il problema dei minimi quadrati

generato da modelli che a turno dipendono dal parametro θ sconosciuto. Acausa del rumore di osservazione o inaccuratezze del modello, possiamo osser-vare solo una versione perturbata di s[n], che chiameremo x[n]. Lo stimatorea minimi quadrati (LSE ) di θ, calcola il valore che rende s[n] il piu vicinopossibile ai dati osservati x[n]. Il criterio con cui si valuta la vicinanza fra idati e l’errore LS

J(θ) =N−1∑n=0

(x[n]− s[n])2 (7.1)

che dipende da θ e da s[n]. Si noti che non e stata fatta nessuna assunzioneprobabilistica sui dati x[n]. Pertanto le performance del metodo dipendonodalle proprieta del rumore e dalle caratteristiche del modello utilizzato.Ovviamente, per via della stessa definizione del metodo, possiamo distingue-re il problema in due classi: se il segnale ha un comportamento pressochelineare in θ, otteniamo per J(θ) un comportamento quadratico (LLS, LinearLeast Square Problem); altrimenti abbiamo un problema non lineare, che


per la soluzione richiede l’utilizzo di metodi esaustivi, in cui vengono testatiuna serie di valori presenti in una griglia (grid searches), oppure di metodiiterativi.

Esempio 7.1 - Segnale sinusoidale

Sia s[n] = A cos 2πf0n, in cui sia nota la frequenza f0 e l’ampiezza A sia ilparametro da stimare. Si ha

J(A) =N−1∑n=0

(x[n]− A cos 2πf0n)2,

che e abbastanza semplice da differenziare e risolvere, in quanto J(A) e qua-dratico in A. Tuttavia, se dovessimo stimare f0 invece di A, il problemasarebbe di soluzione meno semplice, in quanto abbiamo a che fare con unproblema non lineare.Infine, un’ultima possibilita e, nel caso di stima di parametri vettoriali, ilproblema di valutare mediante l’LSE entrambi i parametri, nel qual caso

J(A, f0) =N−1∑n=0

(x[n]− A cos 2πf0n)2

e quadratico in A ma non quadratico in f0. Quando il problema, come inquest’ultimo caso, e lineare in alcuni parametri e non lineare in altri, vienedenominato problema ai minimi quadrati separabili.

7.2 Minimi quadrati lineari - caso scalare

Per questo tipo di problema, dobbiamo assumere

s[n] = θh[n] (7.2)

in cui h[n] e una sequenza nota. Il criterio LS diventa

J(θ) =N−1∑n=0

(x[n]− θh[n])2. (7.3)

La minimizzazione produce l’LSE

θ =

N−1∑n=0

x[n]h[n]

N−1∑n=0

h2[n]

. (7.4)


Sostituendo la 7.4 nella 7.3,

Jmin = J(θ) =N−1∑n=0

(x[n]− θh[n])(x[n]− θh[n]))

=N−1∑n=0

x[n](x[n]− θh[n])− θ

N−1∑n=0

h[n](x[n]− θh[n])

︸︷︷︸S

=N−1∑n=0

x2[n]− θ

N−1∑n=0

x[n]h[n]. (7.5)

L’ultimo passaggio segue dal fatto che la sommatoria S e nulla (si sostituiscaθ per verificarlo). Utilizzando la 7.4, possiamo riscrivere la 7.5 come

Jmin =N−1∑n=0

x2[n]−

(N−1∑n=0

x[n]h[n]

)2

N−1∑n=0

h2[n]

. (7.6)

che in definitiva puo essere riscritta come

Jmin = Ex −

(N−1∑n=0

x[n]h[n]

)2

Eh

. (7.7)

Dunque in generale possiamo dire che

0 ≤ Jmin ≤ Ex.

7.3 Minimi quadrati lineari - caso vettoriale

L’estensione di questa tecnica al caso vettoriale segue abbastanza semplice-mente. Perche il segnale s = [s[0]s[1] . . . s[N − 1]]T sia lineare nel parametroθ, vettore di dimensioni p× 1, dev’essere

s = Hθ (7.8)

dove H e una matrice N × p nota (N > p) di rango p, detta matrice diosservazione. L’LSE viene trovato minimizzando

J(θ) =N−1∑n=0

(x[n]− s[n])2

= (x−Hθ)T (x−Hθ). (7.9)


Riscriviamo la 7.9 come

J(θ) = xTx− xTHθ − θTHTx + θTHTHθ

= xTx− 2xTHθ − θTHTHθ.

Il gradiente e∂J(θ)

∂θ= −2HTx + 2HTHθ,

che posto uguale a zero, restituisce come soluzione

θ = (HTH)−1HTx. (7.10)

L’equazione HTHθ = HTx risolta in θ e denominata equazione normale;l’ipotesi che la matrice H abbia rango pari a p, assicura l’invertibilita di(HTH).Abbiamo trovato uno stimatore che ha la stessa espressione trovata per unostimatore ottimo mediante la tecnica BLUE; tuttavia, affinche il θ descrittonella 7.10 sia il BLUE, dev’essere E(x) = Hθ e Cx = σ2I, ed inoltre e richiestoche x sia Gaussiano.Si puo dimostrare che manipolando le equazioni 7.9 e 7.10, si ottiene lasoluzione generale per il metodo dei minimi quadrati lineari

Jmin = xT (x−Hθ) (7.11)

Il problema dei minimi quadrati lineari si estende immediatamente adun problema a minimi quadrati pesati, (weighted LS). Invece di minimizzaredirettamente l’equazione 7.9, possiamo introdurre una matrice W di pesi,N ×N e simmetrica:

J(θ) = (x−Hθ)TW(x−Hθ). (7.12)

La matrice W viene introdotta con l’obiettivo di dare piu importanza ai datiritenuti piu attendibili.In questo caso l’LSE e

θ = (HTWH)−1HTWx (7.13)

e sostituendo, si ottiene l’errore minimo

Jmin = xT (W−WH(HTWH)−1HTW)x. (7.14)


7.4 Interpretazione Geometrica

Riesaminiamo l’LS da un punto di vista geometrico, dal quale possono risul-tare piu chiare le proprieta dello stimatore.Ricordiamo che possiamo scrivere il modello di segnale come s = Hθ. Sianohi le colonne della matrice H, per cui abbiamo

s = [ h1 h2 . . .hp ]

θ1

θ2...θp

=

p∑i=1

θihi (7.15)

in modo tale che il modello del segnale possa essere visto come una combi-nazione lineare delle hi.

Esempio 7.2 - Sovrapposizione di sinusoidi

Supponiamo che il modello del segnale sia

s[n] = a cos 2πf0n + b sin 2πf0n, n = 0, 1, . . . , N − 1

in cui f0 identifica la frequenza conosciuta, mentre θ = [a b]T e il vettore chedev’essere stimato. In forma matriciale abbiamo

s[0]s[1]...

s[N − 1]

=

1 0cos 2πf0 sin 2πf0

......

cos 2πf0(N − 1) sin 2πf0(N − 1)

[ab

]. (7.16)

Le colonne h1 e h2 di H sono composte da campioni di sequenze sinusoidalie cosinusoidali. Possiamo inoltre riscrivere la 7.9 come

J(θ) = ‖ x−Hθ ‖2

= ‖ x−p∑

i=1

θihi ‖2 . (7.17)

Si puo vedere dalla 7.17 come l’LLS sia un procedimento che minimizza il

quadrato della distanza fra i dati contenuti nel vettore x e il vettore

p∑i=1

θih,


che dev’essere una combinazione lineare delle colonne di H (si ricordi cheH tale matrice ha rango p, pertanto genera uno spazio p−dimensionale, Sp

sottospazio di Rn). La figura 7.2 mostra un esempio per N = 3 e p = 2.Si noti che tutte le possibili scelte per θ1 e θ2 producono dei segnali cheappartengono al sottospazio S2, mentre in generale x non appartiene adesso. Dovrebbe essere abbastanza chiaro a questo punto che il vettore s nelsottospazio S2 e il piu vicino - in senso Euclideo - alla componente di x nelsottospazio (si osservi la figura 7.2(b)). In altri termini, s e la proiezioneortogonale di x sul sottospazio S2. Per determinare s per questo esempio

(a) Sottospazio del segnale (b) Proiezione ortogonale del segnale

Figura 7.2: Interpretazione geometrica dell’LLS in R3

utilizziamo proprio questa condizione di ortogonalita. Questa condizioneprevede che il vettore dell’errore sia ortogonale al sottospazio del segnale

(x− s) ⊥ S2.

Affinche questo sia vero dobbiamo avere

(x− s) ⊥ h1

(x− s) ⊥ h2,

in quanto il vettore errore dev’essere ortogonale a qualsiasi combinazionelineare di h1 e h2. Dalla definizione di ortogonalita,

(x− s)Th1 = 0

(x− s)Th2 = 0.

Sia s = θ1h1 + θ2h2, abbiamo

(x− θ1h1 − θ2h2)Th1 = 0


(x− θ1h1 − θ2h2)Th2 = 0.

In forma matriciale(x−Hθ)Th1 = 0

(x−Hθ)Th2 = 0.

Combinando le due equazioni si ottiene

(x−Hθ)TH = 0T

e in definitiva otteniamo l’LSE come

θ = (HTH)−1HTx.

Si noti che se ε = x−Hθ indica il vettore dell’errore, allora l’LSE puo esseretrovato dalla condizione

εTH = 0T .

Questa condizione, meglio nota come principio di ortogonalita, significa chel’errore dev’essere ortogonale alle colonne di H. In pratica, il principio cisuggerisce che l’errore rappresenta la parte di x che non puo essere descrittadal modello assunto per il segnale.Il minimo errore LS e dato da ‖ x− s ‖2,

‖ x− s ‖2 = ‖ x−Hθ ‖2

= (x−Hθ)T (x−Hθ).

In definitiva, possiamo dire che una soluzione di tipo LS puo essere interpre-tata come un problema di approssimazione di un vettore x in Rn con un altrovettore s, che e una combinazione lineare dei vettori {h1,h2, . . . ,hp}, in unsottospazio p−dimensionale di Rn. La soluzione e appunto un vettore s chesia la proiezione ortogonale di x nel sottospazio, cioe il vettore piu vicino insenso Euclideo.


7.5 Principali tecniche iterative LLS

Passiamo ad una breve panoramica sulle principali soluzioni ricorsive appar-tenenti alla categoria delle tecniche a minimi quadrati; in particolare fra essepossiamo elencare

• order recursive least squares;

• sequential least squares

• constrained least squares.

Come abbiamo spesso sottolineato, in molti casi il modello per il segna-le e sconosciuto e dev’essere ipotizzato. Ad esempio, consideriamo i datisperimentali mostrati in figura 7.3(a). Possiamo assumere come modello

s1(t) = A

s2(t) = A + Bt,

con 0 ≤ t ≤ T . Se campioniamo il modello ad intervalli regolari, otteniamoil corrispondente segnale tempo discreto

s1(n) = A

s2(n) = A + Bn.

Usando l’LSE con

H1 =

11...1

H2 =

1 01 1...

...1 N − 1

si ottiene lo stimatoreA = x

e

A2 =2(2N − 1)

N(N + 1)

N−1∑n=0

x[n]− 6

N(N + 1)

N−1∑n=0

nx[n]

B2 = − 6

N(N + 1)

N−1∑n=0

x[n] +12

N(N2 − 1)

N−1∑n=0

nx[n].

Nelle figure 7.3(b) e 7.3(c) si possono vedere s1(t) = A1 ed s2(t) = A2 + B2t,


Figura 7.3: Dati sperimentali approssimati a minimi quadrati


in cui T = 100. Ovviamente, utilizzando due parametri, l’approssimazionemigliora; l’errore LS minimo infatti deve decrescere all’aumentare dei para-metri, fino a saturare, nel senso che ad un certo punto l’aggiunta di parametrinon potra piu migliorare l’approssimazione, andando semplicemente ad ap-prossimare l’errore dovuto al rumore come se fosse appartenente ai dati.Il primo metodo, order recursive least squares, riduce il calcolo dell’LSE ag-giornandolo ad ogni passo. In particolare, il metodo permette di calcolarel’LSE basandosi su una matrice H di dimensioni N × (k + 1) da una soluzio-ne precedente basata su una matrice H di dimensioni N × k. Ad esempio,nell’esempio precedente, un approccio del genere risulta estremamente sem-plice. Supponiamo infatti di alterare l’intervallo di osservazione da [0, N ] a[−M,M ], facendo in modo che sia simmetrico. Adesso la matrice H2 presentale colonne ortogonalizzate,

H2 =

1 −M1 −(M − 1)...

...1 M

.

Pertanto l’LSE puo essere facilmente trovato per il fatto che

H2TH2 =

2M + 1 0

0M∑

n=−M

n2

.

e una matrice diagonale. Le soluzioni sono

A1 =1

2M + 1

M∑n=−M

x[n]

A2 =1

2M + 1

M∑n=−M

x[n]

B2 =

M∑n=−M

nx[n]

M∑n=−M

n2

.

In questo caso l’LSE di A non cambia se aggiungiamo dei parametri al mo-dello. In termini geometrici possiamo dire che l’ortogonalita delle colonne


della matrice ci permette di proiettare x lungo le due direzioni identificatedalle colonne separatamente e successivamente di sovrapporre i risultati (sinoti la somiglianza con il principio di sovrapposizione degli effetti). Ogniproiezione e indipendente dalle altre. In generale, le colonne della matriceH non sono ortogonali, ma possono essere sempre sostituite da un set di pvettori ortogonali (si ricordi infatti che la matrice ha rango p), tramite ilprocesso di ortogonalizzazione di Gram-Schmidt. Nel caso p = 2, la colonnah2 viene proiettata nello spazio ortogonale ad h1; successivamente, essendoh1 ed h′2 ortogonali, la stima LS diventa

s = h1θ1 + h′2θ′2,

in cui h1θ1 e anche l’LSE per il segnale, basato su H = h1. questa proceduraaggiunge ricorsivamente delle componenti al modello; si noti che il numeromassimo di componenti che e possibile aggiungere corrisponde proprio a p,rango della matrice, oltre il quale non e piu possibile aggiungere nessunainformazione aggiuntiva (o equivalentemente, non e possibile aggiungere unacomponente linearmente indipendente).

Il secondo metodo, sequential least squares, nasce invece dalla necessitadi analizzare ed elaborare dati in tempo reale, effettuando delle operazionidi stima senza attendere che il vettore di dati sia stato aquisito completa-mente. Questo approccio produce una sequenza di LSE: piu precisamente,assumendo di aver determinato l’LSE θ basandoci su una sequenza di N − 1campioni, il metodo permette di aggiornare il θ dopo l’assunzione del cam-pione N−esimo. Il risultato puo essere raggiunto semplicemente aggiungendo

Figura 7.4: Stimatore sequential least square


allo stimatore al passo (N − 1)−esimo un termine correttivo, che in generaledeve decrescere al crescere di N in quanto lo stimatore ‘al passo precedente’contiene sempre piu campioni e deve pertanto avere piu ‘peso’. Uno schemagenerale della procedura e riportato in figura 7.4.

Infine, prendiamo il caso in cui la stima sia vincolata ad una certa co-noscenza priori sui dati; ad esempio, potremmo essere nel caso della stimadell’ampiezza di diverse sinusoidi che sappiamo a priori essere la stessa. Inquesto caso, il numero di parametri puo essere ridotto e il problema affron-tato con una tecnica iterativa nota come constrained least squares. Nel casodi vincolo lineare - molto usuale nella pratica - questo approccio porta aduna soluzione in cui lo stimatore vincolato θc risulta una funzione linearedello stimatore θ tramite le matrici H e la matrice A che esprime il vincolo.La figura 7.5 mostra una interpretazione geometrica di questo tipo di solu-zione: l’effetto delle combinazioni lineari prodotte dalle matrici H e A sullostimatore θ sono in effetti delle rotazioni del vettore s, che in tal modo vienevincolato in un certo sottospazio identificato dal vincolo.

(a) LS non vincolato (b) LS vincolato

Figura 7.5: Interpretazione geometrica fra i metodi a minimi quadrati vincolato e nonvincolato

7.6 Problema LS non lineare

Andiamo ad esaminare le caratteristiche del problema dei minimi quadratinel caso non lineare. Ricordiamo che la procedura LS utilizza, per la stimadel parametro θ, il criterio

J = (x− s(θ))T (x− s(θ))

in cui s(θ) e il modello di segnale per x, con la sua dipendenza esplicita dalparametro θ (si noti che se x − s(θ) ∼ N (0, σ2I), l’LSE e anche l’MLE).


In generale non e possibile esprimere s(θ) come una combinazione lineares = Hθ, ma e una funzione N−dimensionale non lineare di θ. In questocaso, la minimizzazione di J diventa molto piu complicata, e molte volteaddirittura impossibile. Questo tipo di problema viene spesso denominatoproblema di regressione non lineare. Nella pratica la minimizzazione deveessere basata su un approccio iterativo: se le dimensioni di θ non superanop = 5, si preferisce spesso utilizzare dei criteri a griglia.Molto brevemente, possiamo dire che in generale esistono 2 strade per ridurrela complessita del problema:

1. trasformazioni di parametri

2. separazione dei parametri

Nel primo caso, si cerca una trasformazione biettiva di θ che produca unmodello lineare nel nuovo spazio. Sia quindi

α = g(θ)

in cui g e una funzione p−dimensionale di θ invertibile. Se g e tale che

s(θ(α)) = s(g−1(α)) = Hα

allora il modello sara lineare in α.Possiamo trovare l’LSE lineare di α e dunque l’LSE non lineare di θ da

θ = g−1(α)

in cuiα = (HTH)−1gTx.

Passiamo adesso alla descrizione di alcuni tipici problemi di signal pro-cessing in cui non sono utilizzabili tecniche MVU, e che pertanto vengonorisolti con tecniche LSE: la caratterizzazione statistica del rumore non puoessere conosciuta, o comunque lo stimatore MVU non puo essere trovato orisulta troppo complicato da implementare.

Esempio 7.3 - Digital Filter Design

Uno dei problemi piu comuni nell’elaborazione digitale del segnale riguardala progettazione di filtri digitali la cui risposta in frequenza sia rispondentea certe specifiche di progetto. Una strada alternativa, illustrata nel prossimoesempio, puo essere la specifica della risposta impulsiva.


In generale, una risposta impulsiva di tipo IIR (Infinite Impulse Response),puo essere descritta nel dominio trasformato z dal rapporto

H(z) =B(z)

A(z)

=b[0] + b[1]z−1 + . . . + b[q]z−q

a[0] + a[1]z−1 + . . . + a[p]z−p.

Se la risposta in frequenza desiderata e Hd(f) = Hd(exp[j2πf ]), allora larisposta impulsiva e

hd[n] = F−1{Hd(f)}.La risposta all’impulso e data da una equazione alle differenze ricorsiva,ottenuta dalla trasformata z inversa di H(z), nell’ipotesi a[0] = 1,

h[n] =

−

p∑

k=1

a[k]h[n− k] +

q∑

k=0

b[k]δ[n− k] n ≥ 0

0 n < 0.

La soluzione LS prevede che si scelgano degli {a[k], b[k]} in modo tale daminimizzare

J =N−1∑n=0

(hd[n]− h[n])2

in cui N sia sufficientemente grande in modo da poter considerare h[n] so-stanzialmente nulla. Questo approccio tuttavia produce un problema di tipononlineare (si noti che in questo caso hd[n] rappresenta il dato, mentre h[n]rappresenta il modello per il segnale). Ad esempio si consideri

H(z) =b[0]

1 + a[1]z−1,

allora

h[n] =

{b[0](−a[1])n n ≥ 00 n < 0.

e

J =N−1∑n=0

(hd[n]− b[0](−a[1])n)2,

che e non lineare in a[1]; infatti e proprio la presenza di A che fa in modoche l’errore LS sia non-quadratico (condizione richiesta per la linearita). Percercare di risolvere questo tipo di problema, possiamo filtrare hd[n] e h[n]utilizzando il polinomio A, come mostrato in figura 7.6. A questo punto


(a) Errore LS

(b) Errore LS filtrato

Figura 7.6: Conversione di un problema LS da non lineare a lineare

quindi possiamo minimizzare l’errore LS filtrato

Jf =N−1∑n=0

(hdf[n]− b[n])2

in cui hdf[n] e data da

hdf[n] =

p∑

k=0

a[k]hd[n− k]

e a[0] = 1. Allora l’errore LS filtrato diventa

Jf =N−1∑n=0

(p∑

k=0

a[k]hd[n− k]− b[n]

)2

,

che e una funzione quadratica degli a[k] e b[k].Estraendo dalla sommatoria in k il termine per k = 0 si ottiene

Jf =N−1∑n=0

[hd[n] +

(p∑

k=1

a[k]hd[n− k]− b[n]

)]2

.


Nel minimizzare la precedente espressione rispetto ai coefficienti del filtro, sinoti che b[n] = 0 per n > q. Come risultato abbiamo

Jf (a,b) =

q∑n=0

[hd[n]−

(−

p∑

k=1

a[k]hd[n− k] + b[n]

)]2

+N−1∑

n=q+1

[hd[n]−

(−

p∑

k=1

a[k]hd[n− k]

)]2

. (7.18)

La prima sommatoria puo essere minimizzata (annullata) ponendo

hd[n] +

p∑

k=1

a[k]hd[n− k] = b[n].

Dunque in forma matriciale l’LSE dei coefficienti del numeratore e

b = h + H0a

in cuih = [ hd[0] hd[0] . . . hd[q] ]T

H0 =

0 0 · · · 0hd[0] 0 · · · 0hd[1] hd[0] · · · 0

......

......

hd[q − 1] hd[q − 2] · · · hd[q − p]

.

Le dimensioni del vettore H sono (q + 1) × 1, mentre la matrice H0 ha ledimensioni (q + 1) × p. Per trovare l’LSE dei coefficienti del denominatoredobbiamo minimizzare

Jf (a, b) =N−1∑

n=q+1

[hd[n]−

(p∑

k=1

a[k]hd[n− k]

)]2

(7.19)

= (x−Hθ)T (x−Hθ)

in cui a = θ,x = [ hd[q + 1] hd[q + 2] . . . hd[N − 1] ]T

e

H =

hd[q] hd[q − 1] · · · hd[q − p + 1]hd[q + 1] hd[q] · · · hd[q − p + 2]

......

......

hd[N − 2] hd[N − 3] · · · hd[N − 1− p]

.


L’LSE per questo metodo, comunemente chiamato metodo di Prony, si trovacome

a = (HTH)−1HTx.

Esempio 7.4 - Stima AR di parametri per un modelloARMA

Andiamo adesso a descrivere un metodo per la stima dei parametri AR diun modello ARMA, autoregressive moving average. Un modello ARMA perun processo casuale WSS (Wide Sense Stationary) assume una distribuzionedi potenza

Pxx(f) =σ2

u|B(f)|2|A(f)|2

in cui

B(f) = 1 +

q∑

k=1

b[k] exp(−j2πfk)

A(f) = 1 +

p∑

k=1

a[k] exp(−j2πfk).

I b[k] vengono chiamati coefficienti MA, mentre gli a[k] vengono chiamaticoefficienti AR. Il processo e ottenuto eccitando un filtro causale la cui ri-sposta in frequenza sia B(f)/A(f), con un rumore bianco di varianza σ2

u.Poniamo l’attenzione sulla stima dei parametri a[k], lasciando la stima deiparametri b[k] e σ2

u ad una valutazione successiva. Cio e possibile analizzandoil problema come segue. Estendiamo la PSD al piano z e calcoliamo l’inversadella trasformata z:

Pxx(z) =σ2B(z)B(z−1)

A(z)A(z−1)

La PSD ovviamente puo essere ottenuta valutando P sul circolo unitario delpiano z.Effettuiamo la trasformata inversa della A(z)Pxx(z):

Z−1{A(z)Pxx(z)} = Z−1

{σ2B(z)

B(z−1)

A(z−1)

}.

Essendo il filtro causale, per n < 0

h[n] = Z−1

{B(z)

A(z)

}= 0


e per n > 0

h[−n] = Z−1

{B(z−1)

A(z−1)

}= 0.

Allora

Z−1

{σ2B(z)

B(z−1)

A(z−1)

}= σ2b[n] ? h[−n]

= 0 per n > q.

Continuando abbiamo

Z−1 {A(z)Pxx(z)} = Z−1

{σ2B(z)

B(z−1)

A(z−1)

}

= 0 per n > q.

In definitiva, possiamo scrivere una equazione alle differenze come

p∑

k=0

a[k]rxx[n− k] = 0 per n > q, (7.20)

in cui a[0] = 1. Questa equazione e denominata equazione di Yule-Walkermodificata, ed e identica all’equazione di Yule-Walker ad eccezione del fattoche in quel caso q = 0. In pratica, viene stimato

rxx[k] =1

N

N−1−|k|∑n=0

x[n]x[n + |k|],

assumendo di avere i valori reali di x[n] per n = 0, 1, . . . , N − 1. Sostituendol’ultima espressione nella 7.20,

p∑

k=0

a[k]rxx[n− k] = ε[n], n > q

in cui ε[n] indica l’errore dovuto all’effetto dell’errore di stima. Il modellodiventa

rxx[n] = −p∑

k=1

a[k]rxx[n− k] + ε[n] n > q

che puo essere interpretata come una funzione lineare dei coefficienti AR.L’LSE per la stima degli a[k] minimizza

J =M∑

n=q+1

[rxx[n]−

(−

p∑

k=1

a[k]rxx[n− k]

)]2

= (x−Hθ)T (x−Hθ)


con

x =

rxx[q + 1]rxx[q + 2]

...rxx[M ]

θ =

a[1]a[2]...

a[p]

H =

rxx[q] rxx[q − 1] . . . rxx[q − p + 1]rxx[q + 1] rxx[q] . . . rxx[q − p + 2]

......

......

rxx[M − 1] rxx[M − 2] . . . rxx[M − p]

.

L’LSE si ottiene come θ = (HTH)−1HTx, detta equazione di Yule-Walkermodificata ai minimi quadrati.

Appunti di teoria della stima Parte prima: approccio...

Documents

Transcript of Appunti di teoria della stima Parte prima: approccio...