Statistica spaziale: dati e prime analisi

www.dst.unive.it\ ∼gaetan
Statistica spaziale:
dati e prime analisi
Carlo Gaetan
Dipartimento di Statistica
Università Ca’ Foscari Venezia
Venezia, 2 dicembre 2009
IUAV
1/ 73
www.dst.unive.it\ ∼gaetan
Sommario
Introduzione: tipi di dati spaziali
Dati puntuali
Dati geostatistici
Dati di area
IUAV
2/ 73
www.dst.unive.it\ ∼gaetan
Esempi di dati spaziali (dataset in R)
1. Altezze piezometriche in una regione detta Wolfcamp Aquifer
(Texas, USA)
2. Biomassa di blue grama rilevata in una zona di 200 x 200
metri vicino a Elgin (Arizona, USA)
3. Numero di morti dovuti a SIDS nelle contee dello stato del
North Carolina (USA)
4. Valori dei pixel di una immagine satellitare
5. Posizione e diametri di alberi in una regione di 200 x 200
metri nel southern Georgia (USA)
6. Altri ...
IUAV
3/ 73
www.dst.unive.it\ ∼gaetan
Ingredienti di base I
• Le osservazioni provengono da ben identificati siti in una parte
di uno spazio. Consideremo lo spazio geografico.
• La localizzazione di questi siti sono note e etichettano le
osservazioni.
• Le osservazioni e/o le localizzazioni sono modellate come
variabili casuali.
• Variabili risposta:
I
I
I
univariate, multivariate
continue, categoriali
a valori reali oppure no
localizzazione delle osservazioni sulla variabile risposta:
I
I
I
IUAV
punti: regioni, segmenti, curve
irregolare, su una griglia
di forma regolare oppure no
4/ 73
www.dst.unive.it\ ∼gaetan
Ingredienti di base II
I
spazio euclideo oppure no
• Meccanismo di generazione della localizzazione
I
I
noto, non noto
casuale, non casuale
• Notazioni
I
I
I
I
I
IUAV
spazio delle posizioni possibili: S (consideremo S come un
piano)
punto arbitrario in S: s
regione di studio: D sottoinsieme di S
localizzazione di n osservazioni: s1 , . . . , sn
osservazioni: Z (s1 ), . . . , Z (sn )
5/ 73
www.dst.unive.it\ ∼gaetan
Dati geostatistici: altezze piezometriche
La variabile d’interesse esiste in ogni punto della regione ma si
osserva solo la risposta in un insieme finito di localizzazioni.
Piezometric head measurements
at the Wolfcamp Aquifer (Texas, USA)
●
●
●
●●
●●
100
●
●
●
●
●
●
●
●
●
50
1200
●
●●
●
●
●
●
●
600
●
●
●
●
●
● ●
● ●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
● ●
●
●●
●
● ●
●
●
150
●
100
●
●
0
−50
0
●
x
100
200
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
●
●
●
●
●
●
●
●
−200
●
●
●
●
●
●
−150
200
0
●
●
●
−100
−150
−100
●
●
50
●
●
●
●
●
●
●
●
−50
●
●
●
● ●
● ●
●●
●
●
●● ●
●
● ●
y
●
●
●
●
−100
800
●
●
●●
●
●
●
●
●
●
●
y
1000
●
●
●
400
Piezometric head
●
●
●
●
−200
●
●
●
−300
●
●
●●
●
●
●
●
●
●
●
−100
0
100
x
IUAV
6/ 73
www.dst.unive.it\ ∼gaetan
Dati geostatistici: biomassa del blue grama
Si noti come la localizzazione possa essere regolare.
1500
Bluegrama
●
●
y
30
●
1000
40
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
10
●
●
●
●
●
●
●
●
●
●
●
●
0
●
●
●
●
●
●
●
0
●
●
●
●
500
●
●
●
1000
x
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
1500
●
●
1000
●
●
●
500
●
●
●
●
1500
0
●
0
●●
●
●
2000
●
●
500
●
●
●
y
●
20
biomass
●
●
2000
0
500
1000
1500
x
IUAV
7/ 73
www.dst.unive.it\ ∼gaetan
Dati su reticolo
La variabile d’interesse esiste ed è osservata solo in un insieme
finito di localizzazioni.
0.0
32
0.2
34
0.4
36
0.6
38
0.8
40
1.0
Sudden infant deaths
in North Carolina for 1974−78
−84
IUAV
−82
−80
−78
−76
0.0
0.2
0.4
0.6
0.8
1.0
8/ 73
www.dst.unive.it\ ∼gaetan
Dati puntuali
I dati sono le localizzazioni (supposte casuali). Qualche volta nelle
localizzazioni si osserva una variabile risposta (mark), in questo
caso il diametro.
Locations seedlings and saplings
●
●
●
●
Locations and diameters of Longleaf pine trees
● ● ●●
●
●
●
● ●● ●●● ● ● ●
●
●
●
●
●
● ●
● ● ●●
●
●●
●
●
●●
●● ● ●
●●●●
●● ● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
●
●
●
●●●
●
●
●
●
●
●
●●
● ●●● ●
●● ●
●● ●
●
●●
●
●
●
●
●● ●●●●● ●
● ● ●●●● ●●● ● ●●
●
●●●● ●
●
●
●
●● ● ●
●
● ●
●
●
●
●
●
●
●
●
● ● ● ● ●●●
●●
●●●
●
●
●●
● ● ●
●●●
●
●
●
●
●
●
●
● ●● ● ●
● ●
●
● ●●
●
●●
●
●
●
●
●
●
●
●
●
●● ● ●
● ●
●● ●
● ● ●●
●●●●●● ●●●●●● ● ●
●
●
●● ● ●● ●
●
●● ●
●
●
●
●●●
●
●
● ●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●● ●
●●●
●
● ●●● ●●●●
●● ●
●
● ●●
●
●
● ● ●● ● ●● ●
●●
●
●
●
●
●
●
●
●
●
● ●●
●●
●
●
●● ●
●●
●
●● ● ●
●
●
● ●●
● ●●
●
●
● ● ●● ●
●
●● ●
●
● ● ●●
●
● ● ● ● ●●●
● ●
●
●● ● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
● ●
●
●
●
●●●
●
● ●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●●
●
●
●
●
●●●●
●
●
●
● ●●
●
●
●●
●
●
●
● ●●
●●●
●
●
●
●●●●
● ●
●
●
●● ●
●
●
●
●
IUAV
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●●
●
●
●
●
●
●●
● ●
●
●
●
9/ 73
www.dst.unive.it\ ∼gaetan
Tipi di strutture spaziali I
Legge di Tobler
“Osservazioni prese da siti vicini tendono ad essere più simili di
osservazioni prese a siti distanti.”
• Strutture di larga scala
I
I
I
Funzione media per dati geostatistici
Vettore delle medie per dati areali
Intensità per un processo di punto
• Strutture di piccola scala
I
I
I
IUAV
Variogramma, funzione di covarianza
Matrice di vicinanza per dati areali
Funzione K di Ripley
10/ 73
www.dst.unive.it\ ∼gaetan
Tipi di strutture spaziali II
• Stazionarietà (il comportamento del processo generatore dei
dati è simile in tutti i sottoinsiemi di S)
I
I
struttura di larga scala costante
struttura di piccola scala che dipende dalla posizione relativa
• Isotropia (il processo è stazionario e la struttura di piccola
scala dipende solo dalla distanza tra siti).
IUAV
11/ 73
www.dst.unive.it\ ∼gaetan
Alcuni obiettivi della statistica spaziale I
• Inferenza sulla struttura spaziale (verifica sull’esistenza di una
struttura spaziale, stima di questa). Importante per la validità
delle procedure statistiche !
• Inferenza sulla struttura non spaziale (stima degli effetti di un
trattamento, effetti delle covariate, stima del numero di punti)
Residential burglaries and vehicle
thefts per thousand households, 1980
6
7
8
Zi
9
10
11
12
13
14
under 9.84
9.84 − 11.64
11.64 − 14.81
14.81 − 18.76
over 18.76
11
12
13
14
under 23.08
23.08 − 30.48
30.48 − 39.1
39.1 − 45.83
over 45.83
11
11
12
13
14
under 19.02
19.02 − 29.33
29.33 − 39.03
39.03 − 53.16
over 53.16
15
Household income (in $1,000)
15
15
Housing value (in $1,000)
6
7
8
9
10
11
HOVALi
6
7
8
9
10
11
INCi
Zi = β0 + β1 INCi + β2 HOVALi + ηi
IUAV
12/ 73
www.dst.unive.it\ ∼gaetan
Alcuni obiettivi della statistica spaziale II
• Previsione di variabili non osservate (kriging, disegno spaziale
ovvero come scegliere nuovi siti per nuove osservazioni ovvero
come riposizionare i siti esistenti)
Average daily ozone values over 1987 summer in Chicago
Average daily ozone values: prediction
30
●
●
20
●
●
●
●
●
●
●
●
10
50
●
●
●
●
●
●
●
●
30
●
20
●
●
●
0
35
●
−10
−10
10
●
●
●
0
●
●
40
●
y
●
40
Ozone
45
●
●
●
−20
●
−30
−20
−10
0
10
x
IUAV
20
30
40
−20
30
●
−30
●
●
●
●
●
−20
−10
0
10
20
30
13/ 73
www.dst.unive.it\ ∼gaetan
Alcuni obiettivi della statistica spaziale III
• estrazione di segnale, ricostruzione di immagini
IUAV
x
z
ẑ(0)
ẑ
14/ 73
www.dst.unive.it\ ∼gaetan
Dati puntuali
IUAV
15/ 73
www.dst.unive.it\ ∼gaetan
Processi di punto
• Un processo di punto è un meccanismo generatore che
determina la posizione di un insieme di n punti s1 , . . . , sn nello
spazio.
• Da un punto di vista probabilistico un processo di punto può
essere caratterizzato dal numero di punti che cadono in una
regione prefissata A, N(A). Questo numero è una variabile
casuale. L’area di A sarà indicata con |A|.
Centri di 69 città in un altopiano della Spagna
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
IUAV
●
●
●
●
●●
●
●
●
●
●
16/ 73
www.dst.unive.it\ ∼gaetan
Esempi di configurazioni spaziali I
Casualità completa
Centri di 69 città in un altopiano della Spagna
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
IUAV
●
●
●
●
●●
●
●
●
●
●
●
17/ 73
www.dst.unive.it\ ∼gaetan
Esempi di configurazioni spaziali II
Regolarità
Cellule di una sezione istologica
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
IUAV
●
●
18/ 73
www.dst.unive.it\ ∼gaetan
Esempi di configurazioni spaziali III
Aggregazione
Posizione di alberi di sequoia
●
●
●
●
●
●
●
●
●
●
●●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
● ●
●
●
●
●
●
●
IUAV
19/ 73
www.dst.unive.it\ ∼gaetan
Caratteristiche di un processo di punto I
• N(A) è una variabile casuale con una certa media E(N(A))
• Questa media varia con A. Se consideriamo una regione di
area infinitesima ds possiamo definire l’intensità del processo
come
E (N(ds))
λ(s) = lim
|ds|
|ds|→0
Tale funzione è detta funzione d’intensità del primo ordine.
• Un processo stazionario è tale per cui λ(s) assume un valore
costante λ su tutta la regione ; quindi E (N(A)) = λ|A|
IUAV
20/ 73
www.dst.unive.it\ ∼gaetan
Caratteristiche di un processo di punto II
• La funzione di intensità di secondo ordine caratterizza la
dipendenza spaziale, ed è definita come
λ2 (si , sj ) =
E (N(dsi )E (N(dsj )
|dsi ||dsj |
|dsi |,|dsj |→0
lim
• Per un processo stazionario,
λ2 (si , sj ) = λ2 (si − sj ) = λ2 (h)
con h il vettore differenza tra si ed sj per un processo
stazionario.
IUAV
21/ 73
www.dst.unive.it\ ∼gaetan
Casualità spaziale completa e Processo di Poisson
omogeneo I
• Si consideri una regione prefissata A.
• La casualità spaziale completa (Complete Spatial
Randomness) nella regione A corrisponde alla generazione di
un processo di Poisson omogeneo.
• Un processo di Poisson omogeneo d’intensità λ ha due
caratteristiche:
1. per ogni n punti nella regione “i punti si distribuiscono
uniformemente in A”
2. il numero di punti nella regione A ha una distribuzione di
Poisson con media λ|A| (|A| è l’area di A)
IUAV
22/ 73
www.dst.unive.it\ ∼gaetan
Casualità spaziale completa e Processo di Poisson
omogeneo II
λ=100
λ=200
●
●
●
●
●
●
● ●
●
●●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●●
●
● ●●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●● ● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●●
●
●
●
●
●
● ●
● ●●
●
●● ●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
● ●
● ●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
IUAV
●
●
●
●
●
●
● ●
●
●
●
●
● ●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
23/ 73
www.dst.unive.it\ ∼gaetan
Casualità spaziale completa e Processo di Poisson
omogeneo III
Il numero atteso di linee che si intersecano una regione convessa
del piano A è pari a λp(A), p(A) è il perimetro di A. La lunghezza
attesa totale delle linee che attraversano una regione del piano è
pari a λπ|A|
λ=10
IUAV
24/ 73
www.dst.unive.it\ ∼gaetan
Stima non parametrica della funzione di intensità I
n
X
1
s − si
λ̂(s) =
K
τ2
τ
i=1
dove K (·) è una densità di probabilità bivariata simmetrica rispetto
l’origine, detta kernel. Il parametro τ > 0 è noto come ampiezza di
banda (bandwidth) e determina il grado di “lisciamento”
dell’intensità stimata (rappresenta il raggio del disco centrato su s
entro cui i punti si contribuiscono alla stima di λ(s)).
IUAV
25/ 73
www.dst.unive.it\ ∼gaetan
Stima non parametrica della funzione di intensità II
●
●
●
●
0.8
●
●
100
0.4
0.6
x
0.8
1.0
300
800
400
0.6
900
y
500
0.4
400
100
200
300
400
800
900
500
600
0.2
700
900
800
1000
0
0.2
200
700
●
●
0.0
0.6
y
0.4
●
●
●
300
0.2
●
60
0.0
●
●
●
●
●
●
●
●●
●
●
● ●
●
0.0
●
●
●
● ●
● ●
●●
● ●
●
● ●
●
●
●
●
● ●
●
●●
●
●
● ●●
●
● ● ●
● ●● ● ●
●
● ●
●●
●
●
●●
●●
●
● ●●
●
●●
●
●●
●● ●●●
●
●
● ●
● ●
●
●●
●●
●
●
●
●
● ● ●
● ● ●●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ● ●●
●
●
●
●
● ●● ● ●
●
●● ●
●
●
●
●
●●
●●
●
●
●
●
●●
●
●
●●
●
● ●●
●
● ●●
●
●
●
●
● ●
●
●
●
●● ●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
● ●
●
●
●
●
●●
●
● ●
●
●
●
●
●
●
●
●
●
●
● ●
●
● ● ●
●●
●
●
● ●
●●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●
● ● ●
●●● ●
● ●
● ●
●
●
●●
●
●● ● ●
●
●
●
●
●
●●●
●
●
●
●
●
●●
● ●
●
●
●●
●
●
●●
●
● ●
●●
●
●● ●
●
● ●
●●●●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●●
●
●
●
● ●● ●●
●
●
● ●
●
●
● ●●
●●
●●
●
●
●
●
●
●
●
● ●●
●
●●
●
●● ● ●
●
●
●
●
●●
●
●
●
●
●
● ●●
●
●
● ●
●
●●
●
●
●
●
●
●
●
● ●● ●
●
● ● ●
● ●
●
●
●
●
●
●
●
● ●
●
●
●
●
●● ● ●
●
●
●
● ● ●●
● ●
●
●
●
●
●
●
● ●
●
● ● ●
●
●
●
● ● ●
●
●
●●
●
●
●
●
● ●
● ●
●
●
● ● ●
●
●
●
●
●
1.0
●
●
●
0.8
1.0
τ=0.01
0.0
0.2
0.4
0.6
0.8
1.0
x
Posizione di aceri in un terreno di 19.6 acri (Lansing Woods,
Clinton County, Michigan USA)
IUAV
26/ 73
www.dst.unive.it\ ∼gaetan
Funzione K di Ripley I
• La funzione K di Ripley
K (h) =
E(N0 (h))
λ
dove
I
I
IUAV
N0 (h) è il numero degli eventi che stanno dentro un cerchio di
raggio h con centro un evento arbitrario
λ è l’intensità del processo.
27/ 73
www.dst.unive.it\ ∼gaetan
Funzione K di Ripley II
• Nel caso di un processo di Poisson omogeneo (CSR)
K (h) = πh2
infatti in un cerchio di raggio h in media abbiamo un numero
atteso di punti
E(N0 (h)) = λ × ’Area cerchio di raggio h’ = λπh2
• Nel caso di aggregazione K (h) > πh2
• Nel caso di regolarità K (h) < πh2
IUAV
28/ 73
www.dst.unive.it\ ∼gaetan
Funzione K di Ripley III
LA stima di K è data da
b (h) =
K
1
λ2 |A|
una stima per λ è
IUAV
{#di coppie (si , sj ) con distanza ≤ h}
b= n
λ
|A|
29/ 73
www.dst.unive.it\ ∼gaetan
Funzione K di Ripley IV
Residenze di giovani criminali, i cui crimini sono stati registrati nel
1971 a Cardiff (Galles, UK).
K
X
●
●
● ●
●●
●
●
●
●
●
●●●
●●
●
●
●● ●
●
● ●
●
● ●●
●
●
● ● ●●
●● ● ●
● ●
● ● ● ●● ●●●
●
●
● ●●
●
●
●●
●
●●
●
●
● ●
●●
●
●
●
●
●●
●
●
●
●●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●●
●
●●
●
●●
●
● ●
●
●
●
●
●
●
●
●
● ●
1500
1000
● ●
●
●
●
K(r)
●
●●●
●
●
●
●
●
●
●
●
500
●
●
●●
● ●
●
●
●
●
●
●
●
●
0
●
● ●
●
●
●● ●
●
● ●
●●
●
●
iso
theo
●
0
5
10
15
20
r
Presenza di aggregazione !
IUAV
30/ 73
www.dst.unive.it\ ∼gaetan
Funzione K di Ripley V
0.20
K
Cellule di una sezione istologica
●
iso
theo
●
●
0.15
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
0.05
●
●
0.10
●
●
K(r)
●
●
●
●
●
●
0.00
●
●
0.00
0.05
0.10
0.15
0.20
0.25
r
Presenza di regolarità
IUAV
31/ 73
www.dst.unive.it\ ∼gaetan
Funzione K di Ripley VI
350
K
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
50
●
●
●
●
●
●
250
●
●
●
●
iso
theo
●
●
●
●
●
●
200
●
●
●
K(r)
●
●
●
●
●
●
●●
●
●
●
●
●
● ●
●
●
●
●
0
●
●
●
150
●
100
●
300
Centri di 69 città in un altopiano della Spagna
0
2
4
6
8
10
r
IUAV
32/ 73
www.dst.unive.it\ ∼gaetan
Dati geostatistici
IUAV
33/ 73
www.dst.unive.it\ ∼gaetan
Notazioni
Considereremo una regione S che è un sottoinsieme di R2 (per
semplicità), un punto s avrà coordinate s = (x, y )0 . In n punti
distinti s1 , . . . , sn
●
100
●
●
●
●●
●●
●
●
●
●● ●
●
50
●
●
●●
●
●
●
●
●
●
●
0
●
y
●
●
●
●
●
●
●
−50
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●●
●●
●
●● ●
●
●
−150
−100
●
●
●
●
●●
●
−200
● ●●
●
●
●
●
−100
0
●
●
●
●●
●
●
●
●
●
100
x
vengono considerate n osservazioni Z (s1 ), . . . , Zs (n)
IUAV
Piezometric head measurements
at the Wolfcamp Aquifer (Texas, USA)
34/ 73
www.dst.unive.it\ ∼gaetan
La stima della componente di larga scala I
Si supponga che
Z (s) = µ(s) + ε(s)
con ε(s) componente d’errore. Si vuole trovare una stima µ
b(s) di
µ(s), mediante le n osservazioni Z = (Z( s1 ), . . . , Z (sn ))0
disponibili. Per questo consideremo vari metodi.
Media aritmetica
• Ipotesi: Si suppone che la superficie Z (·) si muova attorno ad
un valore ’centrale’ (µ ’costante’).
• Stima
µ
b(s) =
n
X
Z (si )/n
i=1
o
µ
b(s) = med{Z( s1 ), . . . , Z (sn )}0
IUAV
35/ 73
www.dst.unive.it\ ∼gaetan
La stima della componente di larga scala II
• Proprietà: semplicità, velocità di calcolo.
Media mobile
• Ipotesi: si suppone che la superficie µ(·) sia “sufficientemente
liscia”.
• Stima
µ
b(s) =
n
X
Z (si )I (ds,si ≤ r )/
i=1
n
X
I (ds,si ≤ r )
i=1
dove ds,si = ks0 − si k, oppure (k vicini) se
d(s,s1 ) ≤ · · · ≤ d(s,sn )
µ
b(s) =
n
X
Z (si )I (ds,si ≤ d(s,sk ) )/k
i=1
IUAV
36/ 73
www.dst.unive.it\ ∼gaetan
La stima della componente di larga scala III
• Proprietà: è necessario scegliere r o k (grado di lisciamento),
richiede algoritmi di ordinamento.
Medie pesate seconde le distanze Le singole osservazioni non
contribuiscono più in maniera uguale.
µ
b(s) =
n
X
i=1
−p
ds,s
Z (si )/
i
n
X
−p
ds,s
,
i
p≥0
i=1
Proprietà: semplicità, velocità di calcolo
Dati originari
Media pesata
z
z
x
x
y
IUAV
y
37/ 73
www.dst.unive.it\ ∼gaetan
Covarianza e Covarianza campionaria
Date due variabili casuali (non spaziali) X e Y
• Covarianza:
Cov(X , Y ) = E[(X − E(X ))(Y − E(Y ))]
Cov(X , X ) = E[(X − E(X ))2 ] = Var(X )
(varianza)
• Covarianza campionaria: date n osservazioni {x1 , . . . , xn },
{y1 , . . . , yn }
Pn
(xi − x)(yi − y )
sxy = i=1
n−1
sxx = sx2
IUAV
(varianza campionaria).
38/ 73
www.dst.unive.it\ ∼gaetan
Correlazione e Correlazione campionaria
Date due variabili casuali (non spaziali) X e Y
• Correlazione:
Cov(X , Y )
Corr(X , Y ) = p
Var(X )Var(Y )
misura standardizzata −1 ≤ Corr(X , Y ) ≤ 1
• Correlazione campionaria:
rxy =
sxy
sx sy
misura standardizzata −1 ≤ rxy ≤ 1
IUAV
39/ 73
www.dst.unive.it\ ∼gaetan
Auto covarianza
• Covarianza di una variabile con se stessa ?
• Covarianza di una variabile osservata in un sito (indichiamolo
con s) con una variabile osservata in un altro (indichiamolo
con s 0 )
Cov(Z (s), Z (s 0 )) = E{[Z (s) − E(Z (s))][Z (s 0 ) − E(Z (s 0 ))]}
IUAV
40/ 73
www.dst.unive.it\ ∼gaetan
Auto covarianza campionaria
Consideriamo un campione {z(s1 ), . . . , z(sn )}
• Covarianza di una variabile osservata in un sito (indichiamolo
con si ) con una variabile osservata in un altro (indichiamolo
con sj ). Ammettiamo che abbiano la stessa media
campionaria z.
ssi ,sj = (z(si ) − z)(z(sj ) − z)
• n(n − 1)/2 coppie (tante!)
• non molto utile (ovvero molto variabile) in quanto si riferisce
ad una sola coppia !
• potremmo pensare di raccogliere insieme le coppie con le
stesse caratteristiche ad esempio le coppie i cui siti sono
separati dalla medesima quantità si − sj = h
P
si −sj =h (z(si ) − z)(z(sj ) − z)
sh =
N(h)
IUAV
41/ 73
www.dst.unive.it\ ∼gaetan
(Semi) varianza campionaria delle differenze
E’ calcolata non sui dati ma su tutte le possibili differenze
(z(si ) − z(sj ))2
2
• è una misura di dissimilarità
• coppie vicine nello spazio dovrebbero avere valori simili e
quindi differenze vicine allo zero
IUAV
42/ 73
www.dst.unive.it\ ∼gaetan
Nuvola del variogramma I
Diagramma di dispersione delle coppie
si − sj ,
(z(si ) − z(sj ))2
2
150000
0
50000
semivariance
250000
nuvola del variogramma
0
100
200
300
400
distance
IUAV
43/ 73
www.dst.unive.it\ ∼gaetan
Nuvola del variogramma II
Utile nel rilevare
• valori anomali globali;
1e+06
200
dato errato
●
●●
−100
● ●
●
●
●
●●
●
●
●
●
●
● ●
●
●
●
● ●●
●
●
●● ●●●
●
●
●
●
●●
●
●● ●
●
●●
●
●
●
●●●
●
●
●●
●●
●
●
●
●●
0e+00
−200
●
−200
−100
0
X Coord
IUAV
●
●
semivariance
●
0
●
●
8e+05
●
6e+05
●
●
●
Y Coord
●●
●●
●
●
●●
4e+05
●
●
2e+05
100
●
●
100
0
100
200
300
400
distance
44/ 73
www.dst.unive.it\ ∼gaetan
Nuvola del variogramma III
• valori anomali locali;
200
4e+05
dato errato
●
●
●
●
●
●●
●●
●
●
●
●
● ●
●
●
●
●
●●
●
●
●
●●
●
●
●●
● ●
●
● ●
●
●● ●
●
●
●
●
0e+00
−200
−100
0
X Coord
IUAV
●
●
●
−200
●
●●
● ●
●●
3e+05
● ●●
●
2e+05
0
●●●
●
●
● ●●
−100
Y Coord
●●
●●●
●
●
●
●
●
semivariance
100
●●
●
1e+05
●
●
100
0
100
200
300
400
distance
45/ 73
www.dst.unive.it\ ∼gaetan
Variogramma campionario
Il semi-variogramma campionario riassume la nuvola del
variogramma
X
1
γ
b(h) =
(z(si ) − z(sj ))2
2N(h)
si −sj =h
Nel caso in cui i siti non siano disposti su di un reticolo dobbiamo
creare dei contenitori (bins) per contenere le coppie che sono
separate approssimativamente da h
Piezometric head measurements
at the Wolfcamp Aquifer (Texas, USA)
250000
Variogramma campionario
●
●
●
●
●
●●
800
●
●
●
●
●
●
●
●
●
●
● ●
● ●
●●
●
●
●● ●
●
●
●●
●
●
●
●
●
600
●
●
●
●
● ●
●
● ●
● ●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
● ●
●
●●
●
● ●
●
●
150
●
●
●
●
100
●
●
50
●
●
●
0
−50
50000
●
●
y
1000
●
400
Piezometric head
●
150000
semivariance
1200
●
●
IUAV
−150
−200
−100
0
x
100
●
● ● ●
0
200
−100
−300
●
● ●
●
●
200
0
100
200
300
400
46/ 73
www.dst.unive.it\ ∼gaetan
Il variogramma nelle varie direzioni
Anisotropia: la variabilità dipende dalla direzione
Variogramma campionario nelle quattro direzioni :
N−S, NE−SW, E−W e SE−NW
0°
45°
90°
135°
200000
Piezometric head measurements
at the Wolfcamp Aquifer (Texas, USA)
●
●
●
●
●
●
● ●
● ●
●
●●
●
●
●● ●
●
● ●
●
●●
●
●
●
●
●
600
●
●
●
●
●
● ●
● ●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
● ●
●
●●
●
● ●
●
●
150
●
150000
semivariance
800
●
●
●●
100
●
●
50
●
●
●
0
50000
●
●
●
●
●
y
1000
●
●
●
400
Piezometric head
●
100000
1200
●
●
−50
−300
−150
−200
−100
0
x
100
0
200
−100
200
0
100
200
300
400
distance
IUAV
47/ 73
www.dst.unive.it\ ∼gaetan
Esempio: precipitazioni in Svizzera
250
20000
I dati si riferiscono alla pioggia caduta sulla Svizzera l’8 maggio
1986, giorno in cui la nube di Chernobyl si è trovata al di sopra
dell’Europa centrale. I dati sono stati oggetti di una competizione
per trovare il miglior previsore di questi.
100
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●● ● ●
●
●●
●
●
●
●
●
●
●●
● ●
●
●
●
●
●
●
●
●
●
●
15000
150
●●
●
●
●
50
Y Coord
●
●
●
●
●
●●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
10000
●
●
●●
● ●
●
semi−variogram
200
●●
●
●
●
●
● ●
●
●
●
●
●
●
●
−50
5000
0
●
●
0
IUAV
●
50
100
150
200
X Coord
250
300
350
50
100
150
h
200
48/ 73
www.dst.unive.it\ ∼gaetan
Modelli per la covarianza I
• Consideriamo la covarianza tra due siti Cov(Z (s), Z (s 0 )) in
generale questa dipenderà da s e s 0 .
• Supponiamo ora che Z (s) e Z (s 0 ) abbiano la stessa media e
che la covarianza dipenda solo da h = s − s 0
C (s − s 0 ) = Cov(Z (s), Z (s 0 )).
Tale ipotesi è detta stazionarietà,
Campo di colza in Lombardia
IUAV
Campo aleatorio stazionario
49/ 73
www.dst.unive.it\ ∼gaetan
Modelli per la covarianza II
• Introduciamo un’ulteriore restrizione la covarianza dipende
solo dalla distanza khk = ks − s 0 k
C (ks − s 0 k) = Cov(Z (s), Z (s 0 )).
’Campi’ lunari
IUAV
Campo aleatorio isotropico
50/ 73
www.dst.unive.it\ ∼gaetan
Variogramma (teorico)
• Definizione
Si supponga che tutti i siti s1 e s2
E(Z (s1 )) = E(Z (s2 ))
2γ(s1 − s2 ) = Var(Z (s1 ) − Z (s2 )),
La quantità 2γ è detta variogramma e γ è detto
semivariogramma.
• Relazione tra variogramma e covariogramma
Se il processo è stazionario
γ(h) = (σ 2 − C (h)).
IUAV
51/ 73
www.dst.unive.it\ ∼gaetan
Interpretazione del (semi)variogramma
• sill: asintoto del semivariogramma. Solamente i processi
stazionari in senso debole presentano un sill
• range: è la distanza h alla quale γ(h) diviene costante
(C (h) ' 0). Lo definiamo in maniera approssimata perché per
alcuni processi la correlazione è zero solo asintoticamente.
• nugget: discontinuità del variogramma
1.0
0.0
0.5
γ(h)
1.5
2.0
Variogram
0.0
0.5
1.0
1.5
2.0
2.5
3.0
distance
IUAV
52/ 73
www.dst.unive.it\ ∼gaetan
Esempi di (semi)variogramma (processi stazionari)
Modello esponenziale
γ(h) = σ 2 {1 − exp {−h/φ}}
σ 2 = 1, φ = 0.3
γ(h)
0.0
0.0
0.2
0.2
0.4
0.4
γ(h)
0.6
0.6
0.8
0.8
1.0
1.0
Modello pepita
(white noise)
0, h = 0
γ(h) =
σ 2 h 6= 0
σ2 = 1
0.0
0.5
1.0
1.5
distance
IUAV
2.0
2.5
3.0
0.0
0.5
1.0
1.5
2.0
2.5
3.0
distance
53/ 73
www.dst.unive.it\ ∼gaetan
10
8
6
y
4
2
0
0
2
4
y
6
8
10
Regolarità
0
2
4
6
8
10
0
2
4
x
8
10
8
6
y
4
2
0
0
2
4
y
6
8
10
esponenziale (φ = 0.2)
10
pepita
0
2
4
6
8
10
x
esponenziale (φ = 1)
IUAV
6
x
0
2
4
6
8
10
x
esponenziale (φ = 4)
54/ 73
www.dst.unive.it\ ∼gaetan
Esempi di (semi)variogramma (processi stazionari)
Modello wave
γ(h) = σ 2 {1 − φh sin( φh )}
σ 2 = 1, φ = 0.3
γ(h)
0.6
0.0
0.0
0.2
0.2
0.4
0.4
γ(h)
0.6
0.8
0.8
1.0
1.2
1.0
Modello gaussiano
2
γ(h) = σ 2 {1 − exp{− hφ }}
σ 2 = 1, φ = 0.3
0.0
0.5
1.0
1.5
distance
IUAV
2.0
2.5
3.0
0.0
0.5
1.0
1.5
2.0
2.5
3.0
distance
55/ 73
www.dst.unive.it\ ∼gaetan
10
8
6
y
4
2
0
0
2
4
y
6
8
10
Regolarità
0
2
4
6
8
10
0
2
4
x
10
10
8
0
2
4
y
6
8
6
y
4
2
0
0
2
4
6
8
x
wave (φ = .1)
IUAV
8
gaussiano (φ = 2)
10
√
gaussiano (φ = 0.6/ 3)
6
x
10
0
2
4
6
8
10
x
wave (φ = 2)
56/ 73
www.dst.unive.it\ ∼gaetan
Modelli di regressione:
dati geostatistici
variabilità=larga scala+ piccola scala
Z (s) = µ(s) + ε(s)
dove µ(s) componente di grande scala funzione deterministica
sufficientemente ’liscia’ e ε(s) componente di piccola scala
aleatoria.
IUAV
57/ 73
www.dst.unive.it\ ∼gaetan
Superfici di trend e covariate
1. Superfici: la componente di larga scala dipende dalle
coordinate
1.1 µ(s) = β0 + β1 x + β2 y ,
s = (x, y )0
1.2 µ(s) = β0 + β1 x + β2 y + β3 x 2 + β4 xy + β5 y 2
2. Covariate
µ(s) = β0 + β1 u(s).
IUAV
58/ 73
www.dst.unive.it\ ∼gaetan
Esempio: Acquifero del Texas
Perspective plot
Perspective plot
als
residu
terp$z
trend
wolf.in
Z (s)
Z (s)
x
x
x
y
IUAV
Perspective plot
y
=
=
µ(s)
β0 + β1 x + β2 y
y
+
+
ε(s)
ε(s)
59/ 73
www.dst.unive.it\ ∼gaetan
Esempio: Potenziali di ritenzione
Si considerano dei campioni prelevati dal suolo. Successivamente i
campioni vengono saturati con acqua e si misura il potenziale di
ritenzione ad una determinata pressione. Per ogni campione è nota
la frazione granulometrica d’argilla.
Potential
Clay
Residuals
residu
u
z
als
Z (s)
Z (s)
IUAV
x
x
x
y
y
=
=
µ(s)
β0 + β1 u(s)
y
+
+
ε(s)
ε(s)
60/ 73
www.dst.unive.it\ ∼gaetan
Previsione spaziale
• Problema: prevedere il valore Z (s) quando s è un punto non
campionato basandosi sui valori osservati z(s1 ), . . . , z(sn ).
• La soluzione è detta anche interpolazione
• Obbiettivo: costruzione di mappe
Si suppone che il modello sia dato da
Z (s) = µ(s) + ε(s).
Due approcci
b
1. Proporre un modello per µ(s), stimarlo (m(s)),
sottrarlo dai
b
dati εb(s) = Z (s) − m(s),
quindi modellare e prevedere i
residui εb(s)
2. Modellare e stimare simultaneamente µ(s) e ε(s)
IUAV
61/ 73
www.dst.unive.it\ ∼gaetan
Kriging
b (s) non distorto ottimo secondo il criterio
• Previsore lineare Z
dell’errore quadratico medio
b (s)]2
E[Z (s) − Z
• E’ un previsore ottimo solo rispetto al criterio e al modello
scelto
• Basato sulla teoria dei processi stocastici con funzioni di
covarianza
IUAV
62/ 73
www.dst.unive.it\ ∼gaetan
Kriging: caratteristiche
• Il previsore in ogni punto è una media pesata dei valori
osservati z(s1 ), . . . , z(sn )
b (s) = λ1 z(s1 ) + · · · + λn z(sn )
Z
• L’errore quadratico medio di previsione è calcolato
automaticamente nella procedura di calcolo.
• La superficie prevista è liscia il grado di lisciamento dipende
dal variogramma.
• I punti campionati sono previsti esattamente.
IUAV
63/ 73
www.dst.unive.it\ ∼gaetan
Kriging: esempio precipitazione in Svizzera
120
0
●
●
●35
●
60
60
70
80
●
●
70
●●
80
70
●
● ●
●
●
90
80
●
●
80
●
70
●
100
●
7●0
110
80
●
100
10
80
●
●
70
120
0
●
120
●
●
90
0
150
0
80
●
●
●
●
●
●
●
●
90
8●0
80
●
●
●
●
●
●
100
150
300
●
70
●
●
50
●
●
80
90
●
●
●
100
●
0
50
100
150
200
250
300
0
xgrid
Kriging
IUAV
70
●
●
●
●
●
●
●
100
●
90
●
●
0
300
●
●
50
●
● ●
80
●
90
70●
80
80
● 15
80
●
50
●
80
●
70
●
●
●
350
●
●
●
●
●
70
●
●
●●
●
●
●
80
●
●
80
●
●
●
●
●
●
70
●
●
●
●
200
●
●
150
● 250
●
●
90
●
●
●
400
80
●
●
●
●
●
●
●
80
100
●
90
●
●
●
●
●
●
ygrid
●
●
●
●
●
●
100
350
●
●
80
●
80
●
●
●
70
●
●
●
150
●
60
●
●
●
●
●
60
300
●
●
●
●
●
70
●
70
●
●
80
200
●
●
●
110
●
●
●
60●
●
●
●
●
●60●
●
150
●
80
●
80
●
●
●
90
●
●
●
●
●
250
200
●
6●0
●
120
100
●
●
●
150
●
●
●
●
●
ygrid
●
●
●
● ●
●
50
●
●
●
50
●
●
70
200
●
●
●
●
60
●
●
150
50
100
150
200
250
300
xgrid
√
Errore quadratico medio
64/ 73
www.dst.unive.it\ ∼gaetan
Kriging: esempio acquifero del Texas
Errore quadratico medio
di previsione
Kriging Universale
●
●
60
●
●
●
70
50
50
●
●
●
0
y
0
●
40
y
●
60●
●
60
●
●
●
50
●
●
●
●
●
●
●●
50
●
●
●
●
60
●
●
●●
●●
●
●
●
●
●
●
0
●
60
●
30
●
●●
●
50
●
●
●
50
100
●
50
●●
●●
100
●
●
●
●
●
0
●
0
70
●
60●
0
x
0
100
60
60
●
●
●
●
●
●
●
●
●
●
●
50
●
●
●
●
●
6●0
●
●
● ●●
60●
●
−200
●
●
●
●
50
50 ●
●
●
●
●
●
●
●
●
●
●
60●
60
● ●●
●
●
50
0
−150
●
●
●
●
●
●
●
●
−100
−50
−50
90
0●0
●
●
50
0
●
●
●
●
●
60
●
●
●
●
●
−100
●
80
−100
●
●
10
●
●
●
●
●
500●
●
●
●
●
●
●
●
60
●
●
−200
IUAV
●
●
●
●
60
●
●
●
●
●
●
●
●
−150
●
●
●
50●
●
●
●
●
●
−100
0
100
x
65/ 73
www.dst.unive.it\ ∼gaetan
Dati di area
IUAV
66/ 73
www.dst.unive.it\ ∼gaetan
Dati di area: notazioni
Una regione D è suddivisa in n aree Ai , i = 1, . . . , n e viene
osservato il valore su ogni area Zi = Z (Ai ), i ∈ S = {1, 2, · · · , n}.
0.0
0.2
0.4
0.6
under 27.91
27.91 − 29.26
29.26 − 31.02
over 31.02
0.8
1.0
Percentage with blood group A in Eire
0.0
IUAV
0.2
0.4
0.6
0.8
1.0
67/ 73
www.dst.unive.it\ ∼gaetan
Misure di vicinanza I
• Distanze tra i centroidi (dati continui)
• Matrice di prossimità (W = (wij ) di dimensione n × n.
Esempi:
Percentage with blood group A in Eire
Graph
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
IUAV
68/ 73
www.dst.unive.it\ ∼gaetan
Misure di vicinanza II
I
Matrice di contiguità (induce un grafo)
1 se Ai ha un confine in comune con Aj
wij =
0 altrimenti
I
W può essere non simmetrica
wij =
I
IUAV
lij
li
doveP
lij è la lunghezza della frontiera in comune tra Ai e Aj e
li = j lij è il perimetro di Ai
W(k) matrice di distanza al ritardo spaziale k ad esempio
(2)
wij = 1 se Ai e Aj non sono confinanti ma confinano tutte e
due con una regione Ak .
69/ 73
www.dst.unive.it\ ∼gaetan
Lisciamento
Medie mobili spaziali
µ̂i =
zi +
Pn
1+
Pj=1
n
j=1 wij
6100
6050
5900
5950
6000
under 27.91
27.91 − 29.26
29.26 − 31.02
over 31.02
5850
0
IUAV
Percentage with blood group A in Eire
5800
5800
5850
5900
5950
6000
6050
6100
Percentage with blood group A in Eire
under 27.91
27.91 − 29.26
29.26 − 31.02
over 31.02
wij zj
100
200
300
0
100
200
300
70/ 73
www.dst.unive.it\ ∼gaetan
Autocorrelazione e sue misure
Una W -misura dell’ autocovarianza spaziale è data da
Pn Pn
i=1
j=1 wi,j (zi − z)(zj − z)
P
C=
n Pn
i=1
j=1 wi,j
Questa autocovarianza deve essere normalizzata per ottenere una
autocorrelazione. Ad esempio con la varianza campionaria
aggiustata
Pn
(zi − z)2
(n − 1)s 2
2
=
s∗ = i=1
n
n
Ciò ci conduce alla definizione dell’indice di Moran
IUAV
71/ 73
www.dst.unive.it\ ∼gaetan
Indice di Moran
P P
n ni=1 nj=1 wi,j (zi − z)(zj − z)
C
P
M = 2 = Pn Pn
s∗
{ i=1 j=1 wi,j } ni=1 (zi − z)2
• L’indice di Moran M può uscire dall’intervallo [−1, +1] (non è
proprio un’autocorrelazione !)
• È tanto più grande quanto valori in siti vicini sono simili (è
tanto più piccolo quanto valori in siti vicini non sono simili) →
aggregazione o autocorrelazione positiva (repulsione o
autocorrelazione negativa)
• Il caso in cui le osservazioni sono indipendenti corrisponde ad
un valore di InM vicino a 0.
Nel nostro caso M = 0.554.
IUAV
72/ 73
www.dst.unive.it\ ∼gaetan
Indice di Geary
Riprendiamo la misura di dissimilarità (zi − zj )2 .
Una misura media
Pn Pn
2
i=1
j=1 wi,j (zi − zj )
P
P
D=
2 ni=1 nj=1 wij
Una misura media standardizzata è l’indice di Geary
P P
(n − 1) ni=1 nj=1 wi,j (zi − zj )2
D
P
G = 2 = Pn Pn
2
s
2 i=1 j=1 wij
i=1 (zi − z)
IUAV
• assume valori nell’intervallo [0, 2].
• valori prossimi a 1 indicano assenza di autocorrelazione. Valori
inferiori (superiori) a 1 indicano autocorrelazione positiva
(negativa)
• G assume valori vicino allo zero quando a siti vicini
corrispondono valori simili (legge di Tobler)
Nel nostro caso G = 0.38.
73/ 73