Appunti per il corso di
Laboratorio di Calcolo 2
anno 2002-03
Paolo Baldi
1
Analisi multivariata
1.1 Statistica Descrittiva
In questo primo capitolo vedremo dei metodi di analisi di insiemi di dati. Vedremo che queste
tecniche si riconducono essenzialmente a delle questioni di algebra e di geometria, che si
acquisiscono tipicamente nei primi due anni di corso. Cominciamo con un insieme di dati
semplice: vengono osservati n individui e per ciascuno di essi viene osservata una quantità
numerica. Il ricercatore si trova quindi di fronte a n numeri x1 , . . . , xn .
Si chiama media del campione la quantità
n
1
1X
x̄ = (x1 + . . . + xn ) =
xi
n
n
(1.1)
i=1
La media di un campione è uno di quelli che si chiamano indici di centralità. In altre parole la
media indica il centro del campione.
Si vede subito che se poniamo yi = axi + b con a, b ∈ R, allora la media ȳ del campione
y1 , . . . , yn è uguale a a x̄ + b. In altre parole la media è invariante per un cambio di unità di
misura: fare un cambio di scala e poi la media è lo stesso che calcolare la media e poi fare il
cambio di scala. Infatti
ȳ =
n
n
n
n
i=1
i=1
i=1
i=1
1 X 1 X
1X
1X
yi =
(axi + b) = a
xi +
b=
n
n
n
n
= a x̄ + b
Si chiama varianza del campione la quantità
n
(1.2)
1X
(xi − x̄)2
n
i=1
2
Capitolo 1. Analisi multivariata
La varianza si indica σ 2 , oppure σx2 per indicare che si tratta della varianza del campione
x = {x1 , . . . , xn }, nel caso si stia lavorando con più di un campione e ci sia pericolo di confusione.
La radice quadrata σ (oppure σx ) della varianza è lo scarto quadratico o deviazione standard.
La varianza è un indice di dispersione: essa misura se i valori del campione siano o no lontani
dal loro baricentro x̄. Valori grandi della varianza indicano che vi sono valori del campione
lontani da x̄. Viceversa se σ 2 = 0 allora tutti gli xi sono necessariamente uguali a x̄, perché ciò
vorrebbe dire che nella somma della formula (1.2) tutti i termini sono nulli.
Se b ∈ R e yi = xi + b, allora σx2 = σy2 . Infatti sappiamo già che ȳ = x̄ + b e dunque σy2 è
uguale a
n
n
n
1X
1X
1X
2
2
(yi − ȳ) =
(xi + b − x̄ − b) =
(xi − x̄)2 = σx2
n
n
n
i=1
i=1
i=1
Del resto, intuitivamente, sommare b a tutti gli elementi del campione significa spostare di b
anche il baricentro x̄, ma ciò non cambia la dispersione del campione rispetto al baricentro. È
anche facile verificare che se yi = a xi allora σy2 = a 2 σx2 . Un’altra espressione per la varianza
si ottiene sviluppando il quadrato nella (1.2):
n
n
n
1 X 2 2x̄ X
1X
2
(xi − x̄) =
xi −
xi + x̄ 2 =
n
n
n
i=1
=
1
n
i=1
n
X
i=1
i=1
xi2 − x̄ 2
Media e varianza danno quindi delle indicazioni utili sui dati: ne indicano il valore centrale e
una misura della dispersione dei valori rispetto al valore centrale. Vediamo ora cosa si può dire
quando su ogni individuo vengono fatte due misurazioni. Nel prossimo paragrafo affronteremo
invece il caso, più interessante in cui, per ogni individuo, vengono invece misurate n quantità
numeriche.
In questo caso quindi i dati si presentano sotto forma di vettori (x1 , y1 ), (x2 , y2 ), . . . , (xn , yn ).
Chiediamoci se vi sia una relazione di tipo lineare tra i due caratteri, ovvero se esistano due
numeri a e b tali che, grosso modo, si abbia yi = axi + b, i = 1, . . . , n. Dobbiamo cioè cercare
la retta ax + b = y che approssima meglio y come funzione di x (e la chiameremo la retta di
regressione). Più precisamente: in generale non esisteranno dei numeri a, b ∈ R tali che si
abbia esattamente yi = axi + b per ogni i = 1, . . . , n. Esisterà però una retta che rende minimi
gli scarti tra i valori yi e axi + b, in un certo senso la retta che meglio delle altre spiega i valori
yi come funzione lineare-affine degli xi .
Cerchiamo dunque due numeri a, b ∈ R in maniera che la quantità
n
X
(yi − axi − b)2
S(a, b) =
i=1
sia minima. Conviene sommare e sottrarre le medie all’interno del quadrato, cioè
S(a, b) =
n
X
i=1
((yi − ȳ) − a(xi − x̄) + (ȳ − b − a x̄))2 .
1.1 Statistica Descrittiva
3
Sviluppando il quadrato si trova
=
n
X
i=1
S(a, b) =
n
X
2
2
(yi − ȳ) + a
(xi − x̄)2 +
+n(ȳ − b − a x̄)2 − 2a
+(ȳ − b − a x̄)
n
X
i=1
i=1
n
X
i=1
(yi − ȳ)(xi − x̄)+
(yi − ȳ) −2a(ȳ − b − a x̄)
|
{z
}
=0
Le due quantità indicate sono nulle perché, ad esempio,
n
X
i=1
(yi − ȳ) =
n
X
i=1
n
X
i=1
|
(xi − x̄) .
{z
=0
}
yi − nȳ = nȳ − nȳ = 0.
Se poniamo
n
(1.3)
σxy
1X
=
(xi − x̄)(yi − ȳ),
n
i=1
allora si può scrivere
1
S(a, b) = σy2 + a 2 σx2 − 2aσxy + (ȳ − b − a x̄)2 .
n
Poiché b compare unicamente nell’ultimo termine, è chiaro che, per ottenere il minimo, deve
essere b = ȳ − a x̄. Basta quindi trovare il valore di a che rende minimo il trinomio
a → a 2 σx2 − 2aσxy + σy2 .
Derivando e uguagliando a 0 si trova 2aσx2 − 2σxy = 0. Dunque,
(1.4)
a=
σxy
,
σx2
b = ȳ −
σxy
x̄.
σx2
La retta y = ax + b si chiama la retta di regressione del carattere y sul carattere x. In realtà
i valori di a e b appena calcolati sono solo dei valori nei quali le derivate parziali di S si
annullano (cioè quello che si chiama un punto critico) e ciò in generale non garantisce che
si tratti di un punto di minimo (potrebbe anche essere di massimo o di sella o altro ancora).
Non è difficile però vedere che se la variabile (a, b) tende all’infinito, allora anche il valore
4
Capitolo 1. Analisi multivariata
di S tende all’infinito. Inoltre sappiamo che una tale funzione, continua e che vale infinito
all’infinito, possiede necessariamente un minimo ed inoltre che tale punto è necessariamente
un punto critico. Poiché il punto (a, b) calcolato in (1.4) è l’unico punto critico di S, se ne
deduce che esso è anche l’unico punto di minimo della funzione S. La quantità σxy si chiama
la covarianza di x e y. Sviluppando i prodotti nella (1.3), si trova l’espressione equivalente
n
σxy
1X
=
xi yi − x̄ ȳ.
n
i=1
σ
Osserviamo che il coefficiente angolare a = xy2 della retta di regressione ha lo stesso segno
σx
della covarianza σxy . Se dunque σxy > 0, ciò indica che y tende a crescere al crescere di x.
Ci si può aspettare qualcosa del genere, ad esempio, se x e y fossero rispettivamente il peso e
l’altezza di individui. Se invece σxy < 0, l’effetto di x su y sarà antagonista: al crescere dei
valori di x i valori di y tenderanno a decrescere.
•
•
•
•
•
•
•
•
•
•
• • •
•
• ••
•
• •
•
••
• •
•
•
•
•
•
•
•
•
••
••
•
•
•
•
•
•
• ••
•• • •
•
•
••
•• •• • • • •
•
• •••• •• • • •• •
•
••
••
••
•
•
•
•
••
•
• •••
•
•
•
•
• • •
•
• •
• •
• •
•
•
•• •
•
•
•
•
•• • • •
•
•
•
• • •••• •
•
•
•
•
• ••
••
•• • •
•• • • •
•
• • • • •
• ••
•
• •• • •• •
•••
••
••
•
• •
• ••
• • •
•• • • •
•
••
•
•
•
•
•
•
•
•
•
•
Figura 1.1 Esempio di grafico di due caratteri con una correlazione positiva, con la relativa retta di
regressione.
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
• • •••
•
•
••
.....................
••
.............•
•
•• ••
•
....•
••
•
.....•
..........•...•
•
••
• •.......•................•.....•.........
•
•
•
•
.
.
•
.
.
.
.
.
.
.
.
.
•
.
•
• •• ............... •
• •
••
•
•
•••
• •• • •• • • •....•.......•.................•................. • ••••
•
•
• •
•• • •
•
•
• •..................•..•...•....•.....•....... •
••
•
•
•
.
.
.
•
•
.
.
•
•
.
............•
•
•
••
•• •
............•.•........
•
•
• • • •• ••
•
• .............•.•...........
••
• • ••
•
•
• .............•.....•..•... •
•
• • •
...........•..........
•
•• •
•
• • •
•
•
.....................
•
•
•
• ••
••
..........
•
•• •
• •
•
••
•
•
• ••
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
Figura 1.2 Qui, invece, la correlazione è negativa.
Si chiama coefficiente di correlazione dei caratteri x e y la quantità
σxy
·
ρx,y =
σx σy
1.2 Analisi in componenti principali
5
Si può dimostrare che si ha sempre −1 ≤ ρx,y ≤ 1. Un valore di ρx,y vicino a 0 indica scarsa
correlazione tra i due caratteri. Viceversa valori di ρx,y vicini a −1 indicheranno una forte
correlazione antagonista e valori vicini a 1, al contrario, una forte correlazione all’unisono.
Il pregio del coefficiente ρx,y , rispetto alla covarianza, è che esso è insensibile ai cambiamenti
di scala: se tutti i valori di xi vengono moltiplicati per uno stesso numero a > 0 ed i valori di yi
per uno stesso numero b > 0, il coefficiente di correlazione non cambia. Infatti σxy risulterebbe
moltiplicato per ab, σx per a e σy per b.
1.2 Analisi in componenti principali
Si può dire che lo scopo della statistica descrittiva è di estrarre dai dati raccolti le informazioni
più rilevanti che essi contengono.
Nella situazione più semplice, in cui i dati consistono nel rilevamento di un singolo carattere
per ogni individuo di un campione, media e varianza forniscono delle informazioni semplici sul
baricentro e sulla dispersione dei valori osservati. In presenza di 2 caratteri un’informazione
supplementare è data dalla covarianza (o dal coefficiente di correlazione), che misura la dipendenza tra i caratteri.
Sono però possibili situazioni più complesse in cui occorre considerare più di 2 caratteri per
individuo. In questo caso infatti è più difficile tenere conto delle dipendenze tra i caratteri e
non è più possibile fare un grafico dei dati come nelle figure del paragrafo precedente.
In questo capitolo svilupperemo dei metodi di analisi di più caratteri. A differenza dei casi
semplici del capitolo precedente, in cui calcoli e grafici si sarebbero potuti fare con un po’
di pazienza anche a mano, per i metodi di questo capitolo sarà inevitabile l’uso di calcolatori
e di software specializzato. Lo studio della teoria dovrà dunque avanzare di pari passo con
l’acquisizione della capacità di utilizzare questi strumenti di calcolo.
In generale supporremo che i dati siano della forma seguente. Vengono osservati n individui, per
ciascuno dei quali sono rilevati p caratteri (o variabili) che supporremo tutti numerici. Siamo
quindi in presenza di p variabili che indicheremo X(1) , . . . , X(p) .
Rappresenteremo i dati nella forma della matrice
X = (xij ) i=1,...,n
j =1,...,p
in cui le righe corrispondono agli individui e le colonne ai caratteri: xij sarà dunque il valore
del carattere j -esimo misurato sull’individuo i-esimo.
I dati si possono anche vedere come una nuvola di n punti in Rp (sono i punti le cui coordinate
sono date dalle righe della matrice X dei dati). Indicheremo con x1 , . . . , xn questi vettori.
Viceversa ognuna delle colonne della matrice X contiene i valori assunti da una delle variabili
X(1) , . . . , X(p) .
Possiamo iniziare l’analisi calcolando per ciascuna delle variabili le medie, che indicheremo con
x̄·1 , . . . , x̄·p , e le varianze e poi le covarianze di ognuna delle variabili con le altre. Chiameremo
baricentro dei dati considerati il vettore x̄ di cui i numeri x̄·1 , . . . , x̄·p sono le coordinate.
Si chiama matrice di covarianza delle variabili X (1) , . . . , X(p) la matrice p × p che ha come
elemento di posto ij la covarianza tra la variabile X (i) e la variabile X(j ) . Si tratta quindi di
6
Capitolo 1. Analisi multivariata
1a var.
↓

1◦ individuo →
x11
2◦ individuo →  x21






...





◦
n individuo →
xn1
...
...
...
...
...
p a var.
↓

x1p
x2p 











xnp
Figura 1.3 Il tipico aspetto di una matrice dei dati: con molte righe (quanti sono gli individui recensiti)
rispetto alle colonne. Ognuna delle righe può essere pensata come un vettore di Rp (p è il numero delle
variabili).
una matrice che ha sulla diagonale le varianze delle singole variabili (la covarianza di X (i) con
se stessa non è altro che la varianza di X (i) ) mentre riporta fuori della diagonale le covarianze
delle variabili tra di loro. Si tratta chiaramente di una matrice simmetrica, poiché la covarianza
di X(i) con X(j ) è la stessa di quella di X(j ) con X(i) .
Numericamente gli elementi della matrice di covarianza, che indicheremo con 0, sono espressi
in termini della matrice X dei dati da
n
n
h=1
h=1
1X
1X
0ij =
(xhi − x̄·i) (xhj − x̄·j ) =
xhi xhj − x̄·i x̄·j
n
n
Analogamente si definisce la matrice di correlazione: sarà la matrice che ha come elemento di
posto ij il coefficiente di correlazione tra la variabile X(i) e la variabile X(j ) . Naturalmente gli
elementi sulla diagonale della matrice di correlazione sono tutti = 1, perché il coefficiente di
correlazione tra una variabile e se stessa vale sempre 1.
L’idea dell’analisi multivariata è di proiettare i punti su un piano ed osservare il grafico dei
punti proiettati. Naturalmente facendo così si perde dell’informazione: punti che appaiono
vicini nella proiezione possono in realtà essere molto distanti tra di loro. Per questo sarà
importante
a) innanzitutto cercare il piano migliore su cui effettuare la proiezione;
b) essere in grado di stimare quanto si è perso facendo la proiezione, in modo da valutare
quanto i punti della proiezione siano una rappresentazione fedele dei punti originali.
Discutiamo ora il punto a), partendo dal caso più semplice della proiezione degli n punti su una
retta per l’origine. Cerchiamo cioè la retta ‘‘migliore’’ su cui proiettare i punti, ricordando che
la proiezione di un punto x su una retta è il punto della retta che si trova più vicino a x.
Osserviamo intanto che una retta per l’origine può essere individuata da un vettore v ∈ Rp di
lunghezza 1: dato un tale vettore l’insieme dei punti di Rp della forma tv, al variare di t ∈ R
descrive una retta, e d’altra parte ogni retta passante per l’origine è di questa forma. In questa
situazione la proiezione di un punto x ∈ Rp sulla retta individuata da v è il punto tv (un punto
della retta, quindi) dove t è pari al prodotto scalare hx, v i.
1.2 Analisi in componenti principali
7
Ricordiamo infatti che il prodotto scalare hx, y i di due vettori di Rp è il numero |x ||y | cos θ ,
dove θ è l’angolo compreso tra i due vettori x e y. Poiché nel nostro caso supponiamo |v | = 1,
dunque hx, v i = |x | cos θ. La Figura 1.4 illustra questo rapporto tra prodotto scalare e proiezione
di un punto su una retta.
•x
..
......
... .....
...
...
.
.
...
.
...
...
...
...
.
.
...
.
.
.
...
.
.
...
.
.
.
...
.
.
.
...
.
.
.
.
...
.
. ...
.
.
.
......... .... ...
.
.
.
......
.
.
.
.
.
.
.
.
..
....
.. ..
.. ..
... .....
........
.
... .......
.
.
.
.
.
.
.
.
..
..
.
.. ..............
...
. .....
...
........
...
.
.
.
.
.
.
.
.
.
.
.......
...
........
...
........
.................
.
.
..........
.......
•
-
•v
θ
0
..
.. ..
.. ..
Figura 1.4 Se |v| = 1, il valore del prodotto scalare dei vettori x e v è dato dalla distanza tra l’origine ed
il punto segnato con una freccia. Quest’ultimo è anche la proiezione di x sulla retta individuata da v. Si
tratta chiaramente del punto della retta che si trova più vicino a x.
Il prodotto scalare di due vettori x e v si esprime facilmente anche in termini delle loro coordinate.
Si può infatti mostrare che vale la formula
hx, y i =
(1.5)
p
X
xi yi
i=1
In particolare il prodotto scalare è lineare, cioè si ha
hx, αy + βzi = α hx, y i + β hx, zi
per ogni α, β ∈ R e x, y, z ∈ Rp .
Tornando agli n punti x1 , . . . , xn in Rp , quale sarà il vettore unitario v tale che la loro proiezione
sulla retta individuata da v sia la più fedele possibile ? Ricordando il significato intuitivo
della nozione di varianza, una proiezione sarà tanto più fedele quanto più grande sarà la
dispersione dei punti ottenuti dalla proiezione sulla retta, ovvero quanto più grande sarà la
varianza dell’insieme di punti hv, x1 i, . . . , hv, xn i, che ora costituiscono un campione di numeri
reali. Osserviamo comunque che, a causa della formula (1.5), possiamo vedere questi punti
come i valori di una nuova variabile (o carattere). Si tratterà del carattere che nell’osservazione
i-esima prende il valore
v1 xi 1 + . . . + vp xip
che è combinazione lineare delle variabili X (1) , . . . , X(p) con i coefficienti v1 , . . . , vp .
Cominciamo col calcolare la media µ(v) di questi valori: poiché il prodotto scalare è lineare
n
n
k=1
k=1
1X 1X
µ(v) =
hv, xk i = v,
xk = hv, x̄ i
n
n
8
Capitolo 1. Analisi multivariata
Dunque la media µ(v) dei punti proiettati non è altro che la proiezione sulla retta individuata da
v del baricentro x̄ dei punti x1 , . . . , xn . Calcoliamone la varianza. Anzi cerchiamo una formula
per trovare la covarianza dei punti proiettati su due direzioni v e w, che indicheremo σ (v, w).
Poiché hxk , v i − hx̄, v i = hxk − x̄, v i, si ha
n
σ 2 (v, w) =
1X
(hxk , v i − hx̄, v i)(hxk , wi − hx̄, w i) =
n
k=1
n
1X
=
hxk − x̄, v ihxk − x̄, w i
n
k=1
Sviluppando il prodotto
hxk − x̄, v ihxk − x̄, w i =
=
p
X
i,j =1
p
X
i=1
(xki − x̄·i )vi
p
X
j =1
(xkj − x̄·j )wj =
(xki − x̄·i )(xkj − x̄·j )vi wj
e riprendendo il calcolo e scambiando le due somme vediamo apparire dei termini legati alla
matrice di covarianza 0;
p
n
1XX
σ (v, w) =
(xki − x̄i )(xkj − x̄j )vi wj
n
2
k=1 i,j =1
p
1
=
vi wj
n
i,j =1
|
X
n
X
k=1
(xki − x̄i )(xkj − x̄j ) =
{z
}
=0ij
p
X
0ij vi wj
i,j =1
Quest’ultima espressione si può scrivere in maniera compatta
p
X
i,j =1
0ij vi wj =
p
X
i=1
vi
p
X
j =1
0ij vj = h0v, w i
| {z }
(0v)i
In conclusione (la matrice 0 è simmetrica)
σ 2 (v, w) = h0v, wi = h0w, v i
Scegliendo v = w si ha dunque che la varianza dei punti proiettati nella direzione v è
σ 2 (v) = h0v, v i
1.2 Analisi in componenti principali
9
Poiché una varianza è sempre ≥ 0, si ottiene dunque che, per ogni vettore v (unitario o no),
h0v, v i ≥ 0 .
Dunque la matrice di covarianza 0 è sempre semidefinita positiva e tutti i suoi autovalori sono
≥ 0. Quindi il problema di trovare il vettore unitario v lungo il quale la nuvola di punti
considerata abbia la massima dispersione equivale a trovare il vettore unitario v per il quale la
quantità h0v, v i sia massima.
È abbastanza facile vedere che il vettore v di modulo 1 per cui la quantità h0v, v i è massima è
il vettore v1 , autovettore relativo all’autovalore più grande λ1 di 0. Proveremo questo fatto in
due modi diversi e che fanno uso di due tipi di ragionamento completamente diversi. 1◦ modo.
Indichiamo con v1 , . . . , vn gli autovettori di 0, che possiamo scegliere di modulo unitario. Essi
costituiscono una base e, se v è un vettore tale che |v | = 1, allora si può scrivere
v=
Poiché supponiamo |v | = 1, deve essere
p
X
αi vi
i=1
Pp
2
i=1 αi
2
= 1. Infatti
p
1 = |v | = hv, v i =
X
αi vi ,
j =1
i=1
p
=
p
X
αj vj =
p
X
i,j =1
αi αj hvi , vj i =
X
αi2
i=1
Infatti i vettori v1 , . . . , vp sono a due a due ortogonali e hvi , vj i = 0 se i 6= j , mentre hvi , vj i =
|vi |2 = 1 se i = j ). D’altra parte poiché 0vi = λi vi
p
h0v, v i =
=
p
X
i,j =1
X
αi 0vi ,
i=1
p
X
j =1
λi αi αj hvi , vj i =
αj vj =
p
X
λi αi2
i=1
Se supponiamo α1 = 1, α2 = . . . = αp = 0, e cioè v = v1 , allora
h0v, v i = λ1
Per ogni altro vettore v di modulo 1 d’altra parte si ha, poiché λ1 è il più grande degli autovalori,
h0v, v i =
p
X
i=1
λi αi2
≤ λ1
p
X
i=1
αi2 = λ1
10
Capitolo 1. Analisi multivariata
In conclusione la quantità h0v, v i è sempre ≤ λ1 ed è = λ1 solo se v = v1 . Dunque v = v1 è il
vettore di modulo 1 per cui la quantità σ 2 (v) = h0v, v i è massima.
In realtà siamo stati un po’ inesatti nel ragionamento appena fatto. Perché ? 2◦ modo. Il
problema considerato non è altro che la determinazione del massimo della funzione F (u) =
h0u, ui sulla sfera di raggio 1. Si tratta quindi di un problema di massimo vincolato, che si può
affrontare col metodo dei moltiplicatori di Lagrange. La sfera si può indicare come {G = 1},
dove G(u) = |u|2 . Poiché sappiamo che la sfera è compatta e F è continua, esiste sicuramente
un punto di massimo di F su {G = 1} ed il metodo dei moltiplicatori di Lagrange assicura che
questo massimo deve essere soluzione di
F 0 (u) + λG0 (u) = 0
(1.6)
Calcoliamo il gradiente F 0 . Scrivendo in coordinate si ha
k
k
X
∂F
∂ X
0ij ui uj =
(0`i ui + 0i` ui ) = (20u)`
=
∂u` ∂u`
i,j =1
i=1
Ovvero, in forma vettoriale,
F 0 (u) = 20u
Se ne deduce che anche G0 (u) = 2u, dunque l’equazione dei moltiplicatori (1.6) diventa
0u + λu = 0
ed ha soluzione se e solo se −λ è un autovalore di 0 ed u è l’autovettore corrispondente. Se ui
è l’autovettore associato a λi , allora
F (ui ) = λi
Il massimo è quindi raggiunto in corrispondenza dell’autovettore associato al più grande degli
autovalori ed il massimo è appunto il più grande degli autovalori.
Abbiamo dunque visto una caratterizzazione del più grande autovalore di una matrice simmetrica 0 come il massimo della forma quadratica u → h0u, ui sulla sfera unitaria. Questa
caratterizzazione (che è valida anche se la matrice non è definita positiva come le matrici di
covarianza) è un risultato d’interesse generale.
Esercizio Sapreste esprimere in termini di autovalori la quantittà
sup hAu, ui
|u|=1
quando la matrice A non è simmetrica ?
Una volta stabilito che la direzione individuata da v1 è quella per cui la proiezione dei punti ha
dispersione massima, si può dimostrare allo stesso modo che v2 è il vettore per cui h0v, v i è
1.2 Analisi in componenti principali
11
massima, tra i vettori che sono ortogonali a v1 , che v3 è la direzione ortogonale sia a v1 che a
v2 per cui h0v, v i è massima e così via.
I vettori v1 , . . . , vp si chiamano anche le direzioni principali. Anzi v1 è la prima direzione
principale, v2 la seconda e così via. In conclusione il piano individuato dalle direzioni v1 e
v2 è quello cercato: la proiezione della nuvola di punti su questo piano è quella nella quale i
punti proiettati hanno la massima dispersione. Questo piano si chiama anche il primo piano
principale; anche gli altri piani individuati da coppie di direzioni principali vi , vj vanno tenuti
in considerazione e le proiezioni relative possono fornire informazioni utili; le proiezioni su
questi (che chiameremo piani principali) possono infatti dare informazioni complementari a
quelle fornite dalla proiezione sul primo piano principale.
Prima di vedere degli esempi concreti poniamoci la questione di valutare quanto i valori proiettati
siano una immagine fedele dei dati: come abbiamo già osservato è possibile che punti in realtà
distanti tra di loro appaiano vicini nella proiezione.
Chiameremo dispersione totale dei dati la quantità
n
1X
|xi − x̄ |2
n
i=1
Osserviamo che, se (vi )i=1,...,k è una base formata da vettori ortogonali e di modulo uguale a
1, allora si ha, per ogni z ∈ Rk ,
2
(1.7)
|z| =
p
X
j =1
hzj , vj i2
Infatti (vi )i=1,...,k costituisce una base, per cui si ha
z=
p
X
α i zi
i=1
e per di più si ha αi = hz, vi i. Dunque
p
2
|z| = hz, zi =
=
p
X
i,j =1
X
i=1
hz, vi ivi ,
p
X
j =1
hz, vj ivj =
hz, vi ihz, vj ihvi , vj i =
p
X
αi2
i=1
da cui segue la (1.7) osservando che i prodotti scalari hvi , vj i sono uguali a 0 per i 6= j e uguali
a 1 per i = j . Torniamo al calcolo della dispersione:
n
n
p
1 XX
1X
|xi − x̄ |2 =
hxi − x̄, vj i2 =
n
n
i=1
i=1 j =1
12
Capitolo 1. Analisi multivariata
=
p
n
X
1X
n
j =1
|
i=1
2
hxi − x̄, vj i =
{z
}
p
X
λj
j =1
=dispersione del carattere vj =λj
Dunque la dispersione totale è sempre uguale alla somma degli autovalori della matrice di
covarianza. Quanto vale la dispersione dei punti proiettati sul piano principale ? Basta osservare
che la loro matrice di covarianza è
λ1 0
h0v1 , v1 i h0v1 , v2 i
0̃ =
=
0 λ2
h0v1 , v2 i h0v2 , v2 i
Per il calcolo appena fatto dunque la dispersione dei punti proiettati sul primo piano principale
è uguale alla traccia di 0̃, e dunque a λ1 + λ2 . Più in generale, la dispersione della proiezione
sul piano principale formato dai caratteri vh e vk è uguale a λh + λk . Definiremo la fedeltà della
proiezione su un piano come il rapporto tra la dispersione dei punti proiettati e la dispersione
totale. Per la proiezione sul primo piano principale la fedeltà vale dunque
λ1 + λ2
λ1 + . . . + λp
Si tratta di un numero sempre compreso tra 0 e 1 (gli autovalori sono tutti ≥ 0 e quindi il
denominatore è sempre più grande del numeratore), e che si usa come indicatore della bontà
della proiezione: tanto più la fedeltà risulta prossima a 1, tanto più i punti proiettati costituiscono
una immagine fedele dei dati originali.
1.3 Analisi discriminante
Supponiamo di essere in presenza di misurazioni multidimensionali, come nel paragrafo precedente, ma provenienti da più popolazioni. Cioè le osservazioni sono suddivise in gruppi corrispondenti a popolazioni diverse (e lo studioso sa a quale popolazione appartenga ognuna delle
osservazioni).
Vi è una differenza tra le misurazioni di individui appartenenti a popolazioni diverse ? È possibile decidere, solo a partire dalle misurazioni, a quale popolazione appartenga un singolo
individuo ? Ovvero è possibile determinare una regola che ad ogni individuo, in funzione delle
sue misurazioni, associ la popolazione a cui esso appartiene ?
Rispondere a queste domande ha un indubbio interesse sia pratico che teorico. Ad esempio se
le misurazioni fossero le osservazioni cliniche su un gruppo di pazienti ed i gruppi fossero due,
composti da coloro che sono affetti da una certa malattia e da quelli sani rispettivamente, si
potrebbe determinare una procedura automatica che a partire dalle misurazioni stabilisca se un
nuovo paziente sia malato o no.
Da un punto di vista teorico, invece, è interessante determinare il carattere (cioè la funzione
delle misurazioni) che meglio permette di discriminare tra i vari gruppi.
Supporremo ancora che i dati siano della forma di matrice
X = (xij ) i=1,...,n
j =1,...,p
1.3 Analisi discriminante
13
dove xij indica la misurazione della variabile j -esima per lo i-esimo individuo; n è ancora il
numero totale di individui e p il numero di variabili.
Supporremo però che vi siano m popolazioni (o gruppi) presenti e che i primi n1 individui
appartengano alla popolazione 1, che quelli con indice compreso tra n1 +1 e n1 +n2 appartengano
alla seconda e così via fino a quelli con indice compreso tra n1 + . . . + nm−1 + 1 e n1 + . . . +
nm−1 + nm , che apparterranno alla m-esima popolazione. Quindi n1 è l’effettivo della prima
popolazione, n2 quello della seconda e così via. Per semplicità indicheremo con Ik l’insieme
degli indici compresi tra nk−1 + 1 e nk (e quindi corrispondenti agli individui della k-esima
popolazione).
Cominceremo col chiederci se non vi sia una direzione u che discrimini il meglio possibile le
diverse popolazioni. Ovvero tale che nella proiezione nella direzione u le varie popolazioni
appaiano il più possibile ben separate tra di loro. Si tratta, come si vede, di un problema non
molto diverso da quello del paragrafo precedente, anche se ancora dobbiamo chiarire come
si possa caratterizzare matematicamente una tale direzione u. Anche qui la soluzione verrà
dall’analisi della matrice di covarianza 0.
Ricordiamo che il baricentro dei dati x̄ è stato definito nel paragrafo precedente come il vettore
le cui coordinate x̄·j sono date da
n
1X
x̄·j =
xij
n
i=1
Possiamo però considerare anche i baricentri delle singole popolazioni, cioè i vettori ȳ1 , . . . , ȳm
le cui coordinate sono
ȳ1j =
1 X
1 X
1 X
xij
xij , ȳ2j =
xij , . . . , ȳmj =
n1
n2
nm
i∈I1
i∈I2
i∈Im
La matrice di covarianza 0 si può ora scrivere
n
0hk
(1.8)
1X
(xih − x̄·h )(xik − x̄·k ) =
=
n
1
=
n
i=1
m hX
X
j =1 i∈Ij
i
(xih − x̄·h )(xik − x̄·k )
dove abbiamo semplicemente spezzato la somma tra le varie popolazioni. Studiamo più da
vicino le somme all’interno delle parentesi quadre. Sommando e sottraendo la media della
popolazione j -esima e sviluppando abbiamo
X
i∈Ij
(xih − x̄·h )(xik − x̄·k ) =
X
=
(xih − ȳj h + ȳj h − x̄·h )(xik − ȳj k + ȳj k − x̄·k ) =
i∈Ij
14
Capitolo 1. Analisi multivariata
=
X
X
(xih − ȳj h )(xik − ȳj k ) +
(ȳj h − x̄·h )(ȳj k − x̄·k )+
+
X
i∈Ij
i∈Ij
i∈Ij
X
(ȳj h − x̄·h )(xik − ȳj k ) +
(xih − ȳj h )(ȳj k − x̄·k ) =
i∈Ij
Le ultime due somme nella relazione precedente sono nulle. Infatti, poiché il termine ȳj h − x̄·h
non dipende da i,
X
X
(ȳj h − x̄·h )(xik − ȳj k ) = (ȳj h − x̄·h )
(xik − ȳj k ) =
i∈Ij
i∈Ij
hX
i
xik − nj ȳj k = 0
= (ȳj h − x̄·h )
i∈Ij
(la quantità tra parentesi quadre è nulla, ricordando la definizione dei baricentri ȳj k ). Sostituendo
nella (1.8) abbiamo
0hk
m
i
1 XhX
(xih − x̄·h )(xik − x̄·k ) =
=
n
j =1 i∈Ij
=
m
i
1 XhX
(xih − ȳj h )(xik − ȳj k ) +
n
j =1 i∈Ij
{z
}
|
=Dhk
m
i
1 X hX
+
(ȳj h − x̄·h )(ȳj k − x̄·k )
n
j =1 i∈Ij
{z
}
|
=Thk
Uno sguardo alle somme che definiscono la matrice D mostra che quest’ultima non è altro che
la somma delle covarianze Dentro le singole popolazioni, mentre T descrive l covarianze Tra
le singole popolazioni rispetto al baricentro generale x̄.
Più precisamente D è la matrice di covarianza della nuvola di punti che si ottiene sottraendo ad
ognuno dei gruppi il suo baricentro, mentre T è la matrice di covarianza della nuvola di punti
che si ottiene sostituendo ad ogni osservazione il baricentro della popolazione alla quale essa
appartiene.
La decomposizione
(1.9)
0 =D+T
che abbiamo appena ottenuto si chiama formula di Huygens.
Se ora u è una direzione, abbiamo visto nel paragrafo precedente che la dispersione lungo la
direzione u è data da
h0u, ui = hDu, ui + hT u, ui
1.3 Analisi discriminante
15
La separazione delle popolazioni lungo la direzione u sarà dunque tanto più grande quanto più è
piccola la dispersione interna delle popolazioni hDu, ui e grande la dispersione esterna hT u, ui.
Tenendo conto che comunque deve essere
hDu, ui hT u, ui
=1
+
h0u, ui
h0u, ui
la direzione lungo la quale la separazione tra le popolazioni è la più grande sarà quella per cui
la quantità
hDu, ui
h0u, ui
è minima.
Mostriamo che un tale vettore u è l’autovettore relativo al più piccolo autovalore della matrice
0 −1 D.
La funzione F (u) = hDu,ui
h0u,ui è omogenea di grado 0, cioè tale che F (αu) = F (u) per ogni
α > 0. Dunque il problema di determinarne il minimo è equivalente a quello del minimo di
v → hDv, v i con il vincolo h0v, v i = 1; dopo di che basterà scegliere
u=
v
|v |
per avere un vettore di lunghezza 1.
Abbiamo già visto come si usano i moltiplicatori di Lagrange per questo tipo di problemi.
Poiché sappiamo già che i gradienti delle funzioni G(v) = h0v, v i e F (v) = hDv, v i sono uguali
rispettivamente a G0 (v) = 20v e F 0 (v) = 2Dv rispettivamente, l’equazione dei moltiplicatori
di Lagrange diviene
Dv + λ0v = 0
ovvero
0 −1 Dv = −λv
Dunque v deve essere autovettore della matrice 0 −1 D. Inoltre, poiché si ha
hDv, v i = λh0v, v i = λ
il vettore v per cui il minimo di hDv, v i è raggiunto tra quelli per cui h0v, v i = 1 è necessariamente l’autovettore relativo al più piccolo autovalore della matrice 0 −1 D.
Per completare la nostra analisi occorre studiare gli autovalori di 0 −1 D, che non è in generale
una matrice simmetrica, anche se sia 0 −1 che D lo sono. Il prodotto di due matrici simmetriche
non è una matrice simmetrica, a meno che le due matrici commutino.
L’analisi spettrale di 0 −1 D si può però ricondurre a quella di una matrice simmetrica. Sappiamo
infatti che esiste una matrice, che indicheremo con D 1/2 , che è simmetrica e tale che
D 1/2 D 1/2 = D
16
Capitolo 1. Analisi multivariata
nel senso del prodotto di matrici. Moltiplicando a sinistra ambo i membri nella relazione
0 −1 Dv = −λv
per la matrice D 1/2 otteniamo
D 1/2 0 −1 Dv = −λD 1/2 v
ovvero
(D 1/2 0 −1 D 1/2 )D 1/2 v = λD 1/2 v
In conclusione se v è autovettore per l’autovalore −λ della matrice 0 −1 D allora D 1/2 v è autovettore della matrice simmetrica D 1/2 0 −1 D 1/2 rispetto allo stesso autovalore e viceversa
(veramente questo ragionamento è incompleto: perché ?). Inoltre poiché
hD 1/2 0 −1 D 1/2 v, v i = h0 −1 D 1/2 v, D 1/2 v i ≥ 0
gli autovalori sono tutti ≥ 0.
Abbiamo dunque mostrato che la matrice 0 −1 D ha tutti i suoi autovalori reali e positivi. Altre
direzioni discriminanti si ottengono considerando gli autovettori associati agli autovalori via
via più grandi di 0 −1 D . Si tratta di direzioni lungo le quali la discriminazione tra le varie
popolazioni diminuisce man mano che si considerano autovalori più grandi, ma che possono
sempre dare informazioni utili.
L’analisi discriminante quindi procede come l’ACP: si calcolano i due autovettori, w1 e w2 ,
associati ai due più piccoli autovalori di 0 −1 D, e si proiettano i dati sul piano generato da w1 e
w2 . Sarà ragionevole all’inizio centrare e ridurre le variabili. Ora però c’è una complicazione
supplementare: a differenza di quanto accadeva per l’ACP, w1 e w2 non sono ortogonali. Come
fare allora per calcolare la proiezione ortogonale ?
Basterà cercare due vettori ortogonali che individuino lo stesso piano. Ad esempio v1 più un
altro vettore w, ortogonale a v1 e complanare a v1 e v2 . La scelta naturale è w = v2 − hv1 , v2 iv1 .
Essendo una combinazione lineare di v1 e v2 , w si trova nello stesso piano di v1 e v2 ed inoltre
hv1 , w i = hv1 , v2 i − hv1 , v2 i hv1 , v1 i = 0
| {z }
=1
Dunque w e v1 sono ortogonali. Il piano discriminante è quindi individuato da questi due
vettori, che sono per di più di modulo 1 se si sostituisce w con w/|w |.
2
Integrazione numerica
Scilab prevede delle routine che effettuano il calcolo numerico di un integrale. Ad esempio
-->deff("y=f(x)","y=1/x");intg(1,2,f)
produrrà un valore approssimato (molto bene. . . ) dell’integrale tra 0 e 1 di x1 (cioè log 2). Quali
sono i metodi numerici che i programmi come Scilab usano per calcolare gli integrali ?
2.1 Le somme di Riemann
Il primo metodo che viene in mente sono chiaramente. Nei corsi di Analisi Matematica si dimostra che, se una funzione f è continua, allora le somme di Riemann convergono all’integrale.
Più precisamente, se f : [a, b] → R è una funzione continua e suddividiamo l’intervallo [a, b]
in n intervallini di ampiezza n1 (b − a), allora una somma di Riemann associata a questa suddivisione è la quantità
n
X
1
Sn (f ) = (b − a)
f (ξk )
n
k=1
dove ξi è un punto qualunque nello k-esimo intervallino Ik = [a + ni (b − a), a + i+n 1 (b − a),
i = 0, . . . n − 1. Qual è l’errore che si commette se si approssima l’integrale con una somma di
Riemann ?
Vediamo di stimare, per cominciare, l’errore nel primo intervallino. Per semplicità supporremo
a = 0 e indicheremo h = n1 (b − a) l’ampiezza di questi intervallini. Usando il teorema di
Lagrange, lo scarto tra l’integrale su [0, h] ed il valore hf (ξ ) è
Z
0
h
Z h
f (x) dx − hf (ξ ) = f (x) − f (ξ ) dx =
0
Z h
0
=
f (cx )(x − ξ ) dx 0
18
Capitolo 2. Integrazione numerica
Ora, se poniamo
M1 = sup |f 0 (x)| ,
a≤t≤b
si ha
Z
Z h
f (cx )(x − ξ ) dx ≤ M1
|x − ξ | dx =
0
0
Z h
Z ξ
(x − ξ ) dx +
(ξ − x)| dx =
= M1
h
0
ξ
0
= M1
1
ξ 2 + (h − ξ )2
2
2
1
Si vede facilmente che, al variare di ξ ∈ [0, h], la quantità ξ 2 + (h − ξ )2 può valere al massimo
h2 , per cui abbiamo trovato che
Z
h
0
Ora, ricordando che h = n1 (b − a),
1
f (x) dx − hf (ξ ) ≤ M1 h2
2
n
X
1
f (x) dx − (b − a)
f (ξk ) ≤
n
a
k=1
Z
n
X
1
|f (x) − f (ξk )| dx ≤ M1 h2 · n =
≤
2
Ik
Z
b
k=1
=
1
M1 (b − a)2
2n
Quindi l’errore tende a 0 come n1 (n = numero di intervallini). Questa stima vale qualunque sia
il punto ξ scelto all’interno di ogni intervallino. Osserviamo però che se la funzione è crescente,
scegliendo come ξ l’estremo sinistro si avrà comunque una minorazione, mentre con l’estremo
destro si avrà una maggiorazione.
Abbiamo visto che l’errore tende a 0 come n1 . Non è un gran che. Si può fare di meglio ?
In effetti se si sceglie come ξ il punto di mezzo di ogni intervallino, succede qualcosa. Se ξ =
e supponiamo f due volte derivabile,
Z
h
f (x) dx
0
− hf ( h2 )
Z h
=
(f (x) − f ( h2 )) dx =
0
Z h
0 h
h
h 2
1 00
=
(f ( 2 )(x − 2 ) + 2 f (cx )(x − 2 ) ) dx 0
h
2
2.2 Polinomi d’interpolazione
h
2
dove abbiamo fatto uno sviluppo di f attorno a
Z
19
fino al second’ordine. Ora
h
f 0 ( h2 )(x − h2 ) dx = 0
0
per cui, scrivendo M2 = supa≤t≤b |f 00 (x)|,
Z
h
f (x) dx
0
− hf ( h2 )
1
≤ M2
2
Z
h
Z
=
0
h
1 00
2 f (cx )(x
(x − h2 )2 dx =
0
Quindi in ogni intervallino l’errore è maggiorato da
Z
b
a
f (x) dx −
≤
≤
n Z
X
k=1
Ik
1
24
−
h 2
2 ) dx ≤
1
M 2 h3
24
M 2 h3 e l’errore globale invece
n
X
1
(b − a)
f (ξk ) ≤
n
k=1
|f (x) − f (ξk )| dx ≤
1
1
M2 h3 · n =
M2 (b − a)3
2
24
24n
che è un risultato molto migliore.
Si potrebbe anche pensare di migliorare l’approssimazione approssimando la funzione f con
una funzione lineare a tratti. Ora la superficie corrispondente ad ogni intervallino [xi−1 , xi ]è un
trapezio, le cui basi sono alte f (xi−1 ) e f (xi rispettivamente. Dunque ognuno di questi trapezi
ha area n1 (f (xi−1 ) + f (xi )) e l’approssimazione dell’integrale è
n−1
n
X
1X
1 1
1
(f (xi−1 ) + f (xi )) =
f (xi ) + f (b)
f (a) +
I∼
n
n 2
2
i=1
i=1
2.2 Polinomi d’interpolazione
A pensarci bene le somme di Riemann si possono vedere come una approssimazione della
funzione integranda f mediante una funzione costante a tratti, facile da integrare. Ma ci sono
anche altre funzioni facili da integrare: i polinomi. L’idea quindi è: approssimiamo f con un
polinomio e approssimiamo l’integrale di f su ogni intervallino con l’integrale del polinomio
approssimante. Il metodo dei trapezi rientra in questo tipo di metodi, quando il polinomio
approssimante è di grado 1.
Siano f una funzione e x0 , x1 , . . . , xn numeri reali. Si chiama polinomio interpolante di f di
ordine n (associato a x0 , x1 , . . . , xn ) un polinomio P di grado n che coincida con f nei punti
x0 , x1 , . . . , xn . Due questioni
20
Capitolo 2. Integrazione numerica
a) Esiste sempre ?
b) È unico ?
L’unicità è immediata, dato che se P e Q fossero polinomi interpolanti di f di ordine n, allora
il polinomio P − Q si annullerebbe nei punti x0 , . . . , xn . Ma un polinomio di ordine n non si
può annullare in n + 1 punti distinti senza essere identicamente nullo. Quindi necessariamente
P = Q.
Anche l’esistenza è facile da stabilire, dato che è facile costruire esplicitamente il polinomio
interpolante. Per ogni i = 0, . . . , n, il polinomio
(2.1)
Qi (x) =
Y x − xj
xi − xj
j 6=i
è di grado n e si annulla nei punti xj , j 6= i, mentre è immediato che Qi (xi ) = 1. Dunque il
polinomio
n
n Y
X
X
x − xj
Qi (x)f (xi ) =
f (xi )
P (x) =
xi − xj
i=0
i=0 j 6=i
è di grado n ed è tale che P (xi ) = f (xi ).
Sia ora f una funzione definita su un intervallo [0, h]. Il suo polinomio d’interpolazione di
ordine 2 associato ai punti 0, h2 , h è
P (x) = Q0 (x)f (0) + Q1 (x)f ( h2 ) + Q2 (x)f (h)
dove Q0 , Q1 , Q2 sono i polinomi (2.1) associati ai punti x0 = 0, x1 =
esplicitamente
2
Q0 (x) = 2 (x − h2 )(x − h)
h
4
Q1 (x) = − 2 x(x − h)
h
2
Q0 (x) = 2 x(x − h2 )
h
1
0
h
2,
x2 = h. Più
............................................................
............
.............
..............
..............
.
.
.
.
.
.
.
.
....... ..
....
...........
..................
...... ..
.
..........
................
..... ..
.
.
.
..........
..............
.....
.
.
.
.
...........
..............
......
.
.
..... ...
.
.
.........
..............
..... .
.
.
.
.........
..............
.......
.
......
........
.
.
.
.
.........
.......
.....
..
.
.
.
.
.
.
..........
...........
.......
.
.
.
.
.....
....
.
0
π
Figura 2.2 Grafico della funzione seno con il suo polinomio approssimante di ordine 2 (trattini).
2.2 Polinomi d’interpolazione
1
.
....................
..
.....
.....
.
....
.....
..
.....
.
....
.
.
.
.
.
....
.
....
....
...
..
...
.
...
...
...
...
.
.
...
.
.
.
.
...
.
..
..
...
...
.
...
.
..
.
...
..
...
.
.
...
.
.
...
..
.
.
... ...
...
. ...
... .
... ...
. ....
....
..
.
..
.. ..
.. ....
.. ..
.
.
.
.
.
. ....
.
.
...
...
..
.
...
...
.
.
...
.
.
...
...
..
.
.
.
.
...
.
.
.
.
.
...
.
.
.
.
...
.
.
.
.
.
.
.
...
.
.
.
.
...
.
.
.
.
...
.
.
.
...
...
..
.
.
...
.
.
.
.
.
.
.
...
.
.
.
.
...
.
..
.
.
.
...
.
.
.
.
.
.
...
.
.
.
.
.
. ...
..
......
.
..
..
.... .
.
.
.
.
......
.
.
.
.
.
.
. .
.
....
.......
. . . . . . .
....... ....... .......
0
1
Figura 2.1 Grafico dei tre polinomi Q0 (puntini), Q1 e Q2 (trattini) per h = 1
L’integrale
Z
h
f (x) dx
0
verrebbe quindi approssimato con
Z h
Z h
Z h
h
f (0) Q0 (x) dx + f ( 2 ) Q1 (x) dx + f (h) Q2 (x) dx
0
0
0
Con pazienza calcoliamo
Z
h
Z
1
0
= 2h
0
Z
2
Q0 (x) dx = 2
h
(x −
1
2 )(x
=
h
(x − h2 )(x − h) dx =
0
− 1) dx = 2h
2h( 31
3
4
−
+
Z
1
2)
1
0
=
(x 2 − 23 x + 21 dx =
h
6
mentre
Z
Z h
4
Q1 (x) dx = − 2
x(x − h) dx =
h 0
0
Z 1
Z 1
x 2 − x dx =
= −4h
x(x − 1) dx = −4h
h
0
0
=
e
Z
−4h( 13
h
0
Q2 (x) dx =
−
Z
1
2)
=
2
3
h
h
0
Q0 (x) dx =
h
.
6
21
22
Capitolo 2. Integrazione numerica
Dunque l’approssimazione è
Z h
f (x) dx ' h
0
1
6
f (0) + 23 f ( h2 ) + 16 f (h)
Questa è l’approssimazione di Simpson.
Supponiamo ora di volere approssimare l’integrale di f su un intervallo [a, b]. Suddividiamo
questo intervallo in n sottointervalli di uguale ampiezza e indichiamo con a = x0 < x1 < . . . <
xn = b le estremità di questi sottointervalli. Indichiamo poi con y1 , . . . yn i punti di mezzo di
questi sottointervalli. Dunque yi = 21 (xi−1 + xi ). L’approssimazione di Simpson dell’integrale
di f sull’intervallino [xk−1 , xk ] è
Z xk
1
f (x) dx ' (b − a) 61 f (xk−1 ) + 23 f (yk ) + 16 f (xk )
n
xk−1
mentre, sommando i contributi dei singoli intervallini,
Z b
f (x) dx '
a
n−1
n
1
1
1X
2X
f (xk ) +
f (yk )
' (b − a) (f (a) + f (b)) +
n
6
3
3
k=1
k=1
Insomma è come se si mettesse un peso pari a 61 sugli estremi dell’intervallo, un peso pari a 13
sugli altri estremi degli intervallini ed infine un peso pari a 23 sui punti di mezzo degli intervallini.
Quanto è buona l’approssimazione di Simpson ?
Teorema. 2.1 L’errore dell’approssimazione di Simpson è maggiorato da
1
M4 (b − a)5
4
1780 n
dove M4 = supx∈[a,b] |f (4) (x)|.
Dunque l’approssimazione di Simpson è molto buona se f è almeno 4 volte derivabile e con
una derivata quarta non troppo grande.
Esempio 2.2 Cerchiamo un’approssimazione di
Z 2
1
log 2 =
dx
1 x
L’approssimazione di Simpson con n = 1 dà
1 22 1 1
+
+
= 0.6944
6 33 6 2
Il Teorema 2.1 dice che l’errore è minore di 24/1780 = 0.013 (qui è facile vedere che M4 = 24).
In realtà log 2 = 0.6931 e | log 2 − 0.6944| = 0.0012.
log 2 '
2.3 Un altro approccio: simulazione
23
2.3 Un altro approccio: simulazione
5 Sia f : [0, 1] → R una funzione integrabile (senza nessuna ipotesi di regolarità) e sia (Xn )n
una successione di v.a. indipendenti uniformi su [0, 1]. Allora la successione di v.a. (f (Xn ))n
è ancora formata da v.a. indipendenti, tutte di media E(f (X1 )); per la legge dei grandi numeri
quindi
n
1X
f (Xk )
n
(2.2)
P
→
n→∞
k=1
E(f (X1 )) .
D’altra parte sappiamo che
E(f (X1 )) =
Z
1
f (x) dx .
0
Queste osservazioni suggeriscono un metodo di calcolo numerico dell’integrale di f : basterà
disporre di un generatore aleatorio X1 , X2 , . . . di legge uniforme su [0, 1] e quindi calcolare
n
1X
f (Xk ) .
n
k=1
Questa quantità per n grande è un’approssimazione di
Z
1
f (x) dx .
0
Più in generale, se f è una funzione integrabile su un aperto limitato D ⊂ Rd , allora l’integrale
di f su D si può approssimare numericamente in modo simile: se X1 , X2 , . . . è un generatore
aleatorio uniforme su D, allora
n
1X
f (Xk )
n
k=1
P
→
n→∞
1
|D |
Z
f (x) dx .
D
Ad esempio, sappiamo che se (Xn )n una successione di v.a. indipendenti uniformi su [0, 1] e
a < b, allora Zn = a + (b − a)Xn è una successione di v.a. indipendenti e uniformi sull’intervallo
[a, b]. Dunque, se f è una funzione integrabile su [a, b], allora
n
1X
f (Zk )
n
k=1
to
n→∞
E[f (Z1 )] =
1
b−a
Z
b
f (x) dx
a
Questo algoritmo per il calcolo d’integrali o di medie è un tipico esempio di quelli che si
chiamano metodi Montecarlo. Vedremo ora di sperimentarlo su alcune situazioni concrete, per
vedere quali sono i suoi vantaggi ed i suoi difetti.
24
Capitolo 2. Integrazione numerica
Esempio 2.3 Come fare per calcolare con il metodo Monte Carlo l’integrale
Z
π
sin x dx ?
0
(che sappiamo beve essere uguale a 2) ? Prima di procedere occorre apprendere l’uso dei generatori aleatori in scilab. Questo sarà l’oggetto del prossimo paragrafo. per ora accontentiamoci
di sapere che il comando
-->xx=rand(10000,1,"u")
produce un vettore xx di numeri aleatori e uniformi su [0, 1] dunque con i comandi
-->yy=%pi*xx;app=%pi*mean(sin(yy))
vengono ottenuti prima i numeri yy che sono uniformi su [0, π] e poi il valore app che è
appunto l’approssimazione Monte Carlo dell’integrale. Provate ora a
• Fare questa operazione e controllare il risultato.
• Ripetere questa operazione 200 volte mettendo i valori ottenuti in un vettore e fare poi
l’istogramma dei risultati. Quanto vale la varianza dei numeri ottenuti ?
• Ripetere l’operazione, ma stavolta facendo la media su 40 000 numeri aleatori ogni volta.
Come si è comportata la varianza ?
Come si può valutare la velocità di convergenza ? Come si può paragonare questo metodo con
quelli numerici dei paragrafi precedenti. Vedremo che I metodi Monte Carlo sono molto lenti,
ma hanno anche dei vantaggi di altro tipo.
Lo strumento di base per valutare la velocità di convergenza è il Teorema Limite Centrale.
Infatti, se indichiamo Yi = (b − a)f (Xi ) e con Y n = n1 (Y1 + . . . + Yn la loro media, per il
Teorema Limite Centrale
√
(2.3)
Y1 + . . . + Yn − nµ
∼ N(0, σ 2 )
√
n
n (Y n − µ) =
Rb
dove abbiamo indicato con µ = E(Y1 ) = a f (x) dx l’integrale da approssimare e con σ 2 la
varianza delle v.a. Yi . Quest’ultima si può calcolare facilmente, almeno in teoria, dato che
E[Y12 ]
2
= E[(b − a) f (a + (b − a)Xi )] = (b − a)
Z
b
f (x)2 dx
a
per cui
2
σ = (b − a)
Z
a
b
2
f (x) dx −
Z
b
f (x) dx
a
2
2.3 Un altro approccio: simulazione
25
Per ottenere σ 2 occorre quindi calcolare degli integrali, tra cui proprio quello che speravamo di
ottenere con il metodo Monte Carlo. In realtà vedremo che si può fare a meno del valore esatto
di σ 2 , di cui basta avere una maggiorazione.
√
La (2.3) permette di scrivere che, per n grande, la v.a. σn (Y n − µ) segue approssimativamente
una legg N(0, 1). Dunque se indichiamo con Z una v.a. N(0, 1), per δ > 0 e n grande si ha
P | Y n − µ| ≥
δσ
√
n
' P(|Z | ≥ δ)
Ma
P(|Z | ≥ δ) = P(Z ≤ −δ) + P(Z ≥ δ) = 8(−δ) + 1 − 8(δ) =
= 2(1 − 8(δ))
In conclusione
P |Y n − µ| ≥
δσ
√
n
' 2(1 − 8(δ))
Rb
Se vogliamo che Y n approssimi µ = a f (x) dx con un errore che è ≤ ε con probabilità α,
possiamo procedere nel modo seguente. Prima scegliamo δ in maniera che
2(1 − 8(δ)) = α
ovvero δ deve essere tale che 8(δ) = 1 − α2 . Il che è lo stesso che richiedere δ = φ1−α/2 , quantile
δσ
di ordine 1 − α2 della legge N(0, 1). Dopo di che n dovrà essere scelto in modo che √
≤ ε.
n
Esempio 2.4 (continuazione dell’Esempio 2.3) Quanto deve essere grande n perché Y n approssimi l’integrale con un errore che è ≤ 10−3 con una probabilità del 99% ?
Qui α = 0.99 e quindi 1 − α2 = 0.995. Dalle tavole si ricava φ.995 = 2.576. n dovrà dunque
essere tale che
2.576σ
≤ 10−3
√
n
ovvero
n ≥ 106 · 2.5762 · σ 2
Resta da determinare il valore di σ 2 . Come indicato prima, in realtà basta una maggiorazione,
che in questo caso potrebbe essere
Z π
2
σ ≥π
f (x)2 dx ≥ π 2
0
si ottiene dunque che deve essere
n ≥ 106 · 2.5762 · π 2 ≥ 65 492 484
Il vero valore di σ 2 qui è
π2
2
− 4 = 0.935 e con questo valore si sarebbe ottenuto n ≥ 6 203 138.
26
Capitolo 2. Integrazione numerica
C’è un altro modo di stimare σ 2 , ancora a partire dalla Legge dei Grandi Numeri. In effetti
sappiamo che
n
1X
(Yi − Y n )2 → σ 2
n→∞
n
i=1
È quindi possibile, calcolando la varianza empirica dei risultati della simulazione ricavare
una stima di σ 2 , che può permettere di ridurre l’effetto della maggiorazione, come si vede
nell’Esempio 2.4.
In realtà occorre molta precauzione quando si sostituisce il valore stimato della varianza a quello
vero come in questo caso. Qui però il valore di n è in genere talmente elevato da giustificare
questa operazione.
Osserviamo però che questo metodo di stima della varianza non permette di stabilire a priori
quanto debba essere grande n, dato che il valore approssimato di σ 2 sarà noto solo alla fine.
Come si vede dall’Esempio 2.4 il numero di simulazioni richiesto per ottenere una approssimazione ragionevole è comunque molto elevato. I metodi Monte Carlo sono quindi inferiori
come prestazioni a quelli dei paragrafi precedenti. Hanno però due vantaggi
• Convergono sempre anche se la funzione da integrare non è regolare.
• Funzionano anche per gli integrali in dimensione più grande di 1 (la stima della velocità
di convergenza è esattamente la setssa), laddove gli altri metodi numerici diventano molto più
complicati.
E2.1 Sia
2
f (x) = ex .
Vogliamo calcolare
Z
1
f (x) dx
0
a) Determinare un’approssimazione di I con i metodi delle somme di Riemann con punto
a destra, a sinistra, di mezzo e trapezi, per n = 5.
b) Calcolare per n = 10, 20, 50, 100 le stesse approssimazioni di a).
b1) produrre e stampare una matrice avente come prima colonna il valore di n, come
seconda l’approssimazione ottenuta, come terza colonna la maggiorazione teorica dell’errore
(tranne che per i trapezi) e come quarta colonna l’errore numerico (=valore approssimato-valore
dell’integrale, ottenuto con il comando intg di scilab)
b2) Disegnare un grafico con n in ascisse e le approssimazioni ottenute in ordinate. Aggiungere al grafico la retta orizzontale di ordinata uguale al e le stime degli errori dei due metodi.
(Per il calcolo delle quantità M1 = sup0≤t≤1 |f 0 (t)| e M2 = sup0≤t≤1 |f 00 (t)| si potrà ammettere
che la funzione f ha le derivate di tutti gli ordini che sono delle funzioni crescenti
c) Vogliamo ottenere numericamente un’approssimazione di I con un errore in modulo
≤ 10−4 . Che valore di n bisogna scegliere se si usano i metodi
1. punto di sinistra
2. punto di destra
3. punto di mezzo.
Esercizi
27
d) Scrivere una funzione scilab che, in funzione di n effettua le operazioni precedenti.
E2.2 Ripetere le operazioni precedenti per le funzioni
a)
Z
π
sin x dx
I=
b)
I=
Z
0
π
e−x
2 /2
sin(12 x) dx
0
(per questa funzione si potranno ammettere le approssimazioni M1 = 1, M2 ' 134.)
c)
Z
π
I=
√
x sin x dx
0
E2.3 Per le funzioni degli Esercizi 2.1 e 2.2 calcolare le approssimazioni con il metodo di
Simpson e confrontare per n = 5 e n = 10 e confrontare con i risultati degli altri metodi.
E2.4 Vogliamo calcolare numericamente l’integrale
I=
Z
π2
√
sin x dx
0
a) In quanti sottointervalli occorrerebbe suddividere l’intervallo [0, π 2 ] perché si possa essere
sicuri che l’errore commesso sia ≤ 10−2 con il metodo del punto di sinistra ? Et del punto di
mezzo ? E dei trapezi ? E di Simpson ?
b) Decidiamo piuttosto di calcolare l’integrale per simulazione:
N
X
1
I ' π2
f (Xi ) := IN
N
i=1
dove le v.a. Xi sono uniformi su [0, π 2 ].
b1) A priori, come si deve scegliere N perché l’errore dello stimatore IN sia ≤ 10−2 con
probabilità 99% ?
b2) Calcolare IN per N = 25000 et determinare qual è l’errore, sempre ad un livello 99%, prima
con la maggiorazione delle varianza e poi stimando la varianza a partire dai dati. In questo
caso, confrontare il risultato con quello dei vicini.
c) Improvvisamente una illuminazione! Con un cambio di variabile si ha
I=
Z
π2
0
√
sin x dx = 2
Z
π
y sin y dy
0
Se si è disposti ad una maggiorazione grossolana della derivata quarta della funzione y →
y sin y, in quanti sottointervalli bisogna suddividere [0, π] perché l’errore sia ≤ 10−2 con il
metodo di Simpson ?
28
Capitolo 2. Integrazione numerica
d) Rispondere alle domnde di b) per il calcolo dell’integrale
J=
E2.5 Calcolare
Z
Z
1/π 2
sin √1x dx
0
e−(x
2 +y 2 )/2
dx dy
A
dove A è
a) il quadrato [−1, 1] × [−1, 1];
b) il cerchio di raggio 1 e di centro l’origine.
E2.6 Calcolare
Z
e−|x|
2 /2
dx dy
|x|≤1
2.4 Probabilità e simulazione in scilab
Scilab ha varie funzioni per gli aspetti numerici della Probabilità e della Statistica e le
simulazioni.
Prima di tutto ci sono i calcoli relativi alle funzioni di ripartizione. Ad esempio per le leggi
gaussiane il comando cdfnor produce la funzione di ripartizione ed i quantili delle leggi
gaussiane. L’uso di questi comandi non è proprio amichevole. Vediamo un po’. Se viene dato
il comando
-->[P,Q]=cdfnor("PQ",x,mu,sigma)
P sarà il valore della funzione di ripartizione di una legge gaussiana di media mu e deviazione
standard sigma calcolato in x . Il secondo valore restituito dal programma, Q è semplicemente uguale a 1-P . Attenzione agli aspetti seguenti:
• sigma è la deviazione standard e non la varianza.
• x può essere un vettore e in questo caso anche P e Q saranno dei vettori. Per approfittare
di questa possibilità però anche mu e sigma devono essere dei vettori, anche se le medie e le
varianze sono tutte uguali. Ad esempio per avere i valori della funzione di ripartizione di una
N(0, 1) nei valori 1, 1.5, 2 occorrerà scrivere
-->x=[1,1.5,2];
-->[P,Q]=cdfnor("PQ",x,zeros(x),ones(x));
Cambiando la sintassi il comando cdfnor produce i quantili: scrivendo
-->x=[1,1.5,2];
-->x=cdfnor("X",0,1,P,Q);
2.4 Probabilità e simulazione in scilab
29
il valore ottenuto x sarà il quantile corrispondente alla probabilità P (qui come prima Q deve
essere uguale a 1-P ). Come prima il comando è vettoriale, nel senso che P può essere un
vettore. In questo caso anche tutti gli altri argomenti devono essere dei vettori aventi la stessa
dimensione.
Per distribuzioni diverse dalla normale esistono i comandi appropriati:
cdfbet (leggi beta)
cdfbin (binomiali)
cdfchi (chiquadrato)
cdfchn (chiquadrato decentrato)
cdfgam (gamma)
cdff (F )
cdffnc (F decentrata)
cdfnbn (binomiale negativa)
cdfpoi (Poisson)
cdft (t di student)
La loro sintassi è simile a quella di cdfnor, sia per il calcolo delle funzioni di ripartizione che
per quello dei quantili.
Per la simulazione esistono due generatori aleatori in scilab
Esempio 2.5 (Il grafico dei quantili) Siano X1 , . . . , Xn delle v.a. indipendenti, aventi tutte la
stessa legge (continua o discreta, poco importa). Indichiamo con F la loro f.r. Se poniamo
At =] − ∞, t], allora si vede subito che
E(1At (X1 )) = P(X ≤ t) = F (t) .
Quindi per la (2.2) applicata alla funzione f = 1At si ha
n
(2.4)
1X
1At (Xk )
n
k=1
P
→
n→∞
E(1At (X1 )) = F (t) .
Indichiamo con F̄n (t) il termine a sinistra nella (2.4), cioè
n
1X
1
F̄n (t) =
1At (Xk ) = × numero di indici k tali che Xk ≤ t .
n
n
k=1
La funzione F̄n (t) è una funzione di t che si chiama funzione di ripartizione empirica . Un
altro modo, per certi versi più semplice, di vedere la f.r. empirica è il seguente. Indichiamo con
X(1) , . . . , X(n) i numeri X1 , . . . , Xn riordinati in senso crescente: in particolare X(1) è il valore
più piccolo e X(1) ≤ X(2) ≤ . . . ≤ X(n) . Allora è chiaro che si ha anche

0 se t < X(1)



k
F̄n (t) =
se X(k) ≤ t < X(k+1) , k = 1, . . . , n − 1

n


1 se X(n) ≤ t .
30
Capitolo 2. Integrazione numerica
Infatti se X(k) ≤ t < X(k+1) , ciò vuol dire che vi sono esattamente k indici i tali che Xi ≤ t e
dunque
n
1X
k
1At (Xk ) = ·
n
n
k=1
Osserviamo comunque che, per ogni valore di t, F̄n (t) è una quantità aleatoria, dato che
dipende dai valori X1 , . . . , Xn . Dunque la (2.4) non fa che affermare che, per ogni t fissato, la
f.r. empirica converge in probabilità alla f.r. F .
1
•
•
........
..
. . .•
.. . .
.. . . .
...
.
..
...
x(1) x(2)
...
..
..
.
..
..
x(3)
..
. ..
•. . .
..
.
•.. . .
..
...
...
.
• . . . . ........
.
. . . .
.•. .
.
.
..
.
..
x(4)
x(5)
x(6)
x(7)
Figura 2.3 Esempio di funzione di ripartizione empirica per un campione di rango 7 di una legge N (0, 1),
a confronto con la funzione di ripartizione teorica (a puntini).
Queste osservazioni hanno delle conseguenze pratiche interessanti. Talvolta in presenza di dati
sperimentali X1 , . . . , Xn è interessati a stabilire se essi seguano una distribuzione assegnata,
N(0, 1) ad esempio. Se n è abbastanza grande, per la Legge dei Grandi Numeri, la differenza
tra la f.r. empirica F̄n e quella teorica F è piccola. Poiché sappiamo che F̄n assume il valore ni
nel punto X(i) , ciò vuol dire che i numeri ni e 8(X(i) ) devono essere molto vicini. Applicando
la funzione 8−1 , che è continua, ai due numeri ni e 8(X(i) ), otteniamo che anche i due valori
8−1 ( ni )
e
X(i)
devono trovarsi vicini tra loro. Il punto 8−1 ( ni ) è il numero in cui la f.r. 8 prende il valore ni
e dunque non è altro che il quantile φi/n . Ciò significa che se si disegnano sul piano i punti di
coordinate X(i) e φi/n , essi si devono disporre lungo la bisettrice del primo e terzo quadrante.
Questa osservazione fornisce un metodo per valutare ‘‘ad occhio’’ (è il caso di dirlo) se le
osservazioni seguano o no una legge N(0, 1). Esso naturalmente può essere applicato anche ad
altre distribuzioni continue (esponenziali, χ 2 ,. . . ), semplicemente sostituendo ai quantili φi/n
quelli della f.r. appropriata.
L’interesse di questo metodo grafico è però ancora limitato perché esso permette di vedere se
le osservazioni seguano una prefissata distribuzione, mentre spesso si è interessati piuttosto a
vedere se esse seguono una distribuzione di una determinata famiglia, ad esempio se esse sono
normali N(µ, σ 2 ) per qualche valore ignoto di µ e σ 2 . Gli argomenti precedenti possono però
essere perfezionati per ricavare un metodo grafico per trattare questo problema più complesso
(anche se ci limiteremo, qui, al caso di leggi normali). Basta per questo ricordare la relazione tra
2.4 Probabilità e simulazione in scilab
31
i quantili di una legge N(µ, σ 2 ) e quelli di una N(0, 1): il quantile di ordine ni di una N (µ, σ 2 )
è uguale a σ φi/n + µ Riprendendo il ragionamento precedente, se le osservazioni seguono una
distribuzione N(µ, σ 2 ), allora i numeri
σ φi/n + µ
e
X(i)
devono essere vicini tra loro; in altre parole i punti di coordinate X(i) e φi/n devono allinearsi
lungo la retta di equazione x = σy + µ. Basta quindi disegnare questi n punti sul piano e poi
vedere se essi si dispongono lungo una retta o no (ora non più necessariamente la bisettrice del
primo e terzo quadrante).
.
... . . .
.
.
.
.
.
.
.
.
.....
................
.
.
.
.
.
.
.
.
.
.
.
.
.
........
............
.
.
.
.
.
.
.
.
.
.
.
. ..
. . .
.
Figura 2.4 Grafico dei quantili per un campione ottenuto per simulazione di una legge normale.
Il grafico dei punti (X(i) , 8−1 ( ni )) si chiama, nel gergo della statistica il grafico quantili contro
quantili, poiché si tratta di confrontare i quantili della distribuzione empirica F̄n contro quelli
della N(0, 1). La sua realizzazione manuale può risultare laboriosa, ma di solito tutti i software
statistici hanno prevedono questa operazione.
.
..
........
.
.
.
.
..........
.
.
.
.
.
.
.
.
.....
........
.
.
.
.
.
.
.
.
...
............
. . . . . ...........
.
.
Figura 2.5 Grafico dei quantili per un campione ottenuto per simulazione di una legge χ 2 (4). L’andamento
si discosta da una retta.