Previsione spaziale
L’obbiettivo dell’analisi spaziale è prevedere valori non osservati di una grandezza misurabile in predefinite locazioni spaziali a partire dalle sue misurazioni campionarie rilevate su un numero finito di siti \(s_1,...,s_n\).
Si suppone che il fenomeno sia definito in un regione di studio \(D\) e che possa essere rappresentato da delle coordinate spaziali \(s\).
\[ Y(s):s \in D \]
Quindi \(Y(s_1),...,Y(s_n)\) sono i valori osservati in un insieme di di siti noti (\(s_1,...s_n\)), solitamente \(s_i\) corrisponde una coppia ordinata (ascissa, ordinata) per esempio (longitudine, latitudine).
La previsione può essere puntuale (o intervallare), binaria rispetto a soglie di riferimento e su aree nella regione di studio.
Esistono due approcci per la preisione spaziale:- Previsione stocastica
- Previsione non stocastica
Nella previsione di tipo stocastico si assume l’esistenza di un modello probabilistico, mentre nella previsione di tipo non stocastico non si suppone l’esistenza di un tale modello sottostante. In questo articolo ci occuperemo solo dell’approcio stocastico.
Kriging
Il Kriging è un metodo utilizzato la previsione spaziale introdotto da Matheron e Gandin nei primi anni ’60 formalizzando studi precedenti condotti da Krige (1951).
L’obbiettivo princiaple di questo metodo è quello di ottenere una mappa del processo.
Nella classe dei previsori lineari, il kriging è quello che minimizza l’errore quadratico medio di previsione e, per tal motivo, è considerato un metodo di previsione ottimo (lineare). Il kriging può essere classificato in tre tipologie a seconda delle assunzioni sulla funzione media del processo spaziale che definisce la variabile regionalizzata. In particolare, se il trend è noto il kriging utilizzato è detto semplice, se è ignoto ma costante il metodo di previsione è detto kriging ordinario mentre se è una funzione delle locazioni (quindi ignoto e non costante) si parla di kriging universale.
- Kriging semplice (SK): \(\mu(s)\) nota
- Kriging ordinario (OK): \(\mu(s)=\mu\) ignota
- Kriging universale (UK): \(\mu(s)\) funzione non nota di \(s\)
Il kriging è un metodo di regressione usato nell’ambito dell’analisi spaziale (geostatistica) che permette di interpolare una grandezza nello spazio, minimizzando l’errore quadratico medio.
Un punto di forza di questo tipo di interpolazione è che non solo genera un modello spaziale interpolato, ma genera anche una stima dell’incertezza di ogni punto in quel modello. A differenza della regressione lineare o dell’interpolazione ponderata per la distanza inversa,l’interpolazione di kriging si basa principalmente su osservazioni empiriche, i punti dati del campione osservato, piuttosto che su un modello preassunto. L’interpolazione attribuisce maggiore peso ai punti campione nelle vicinanze di una posizione rispetto a quelli più importanti e, al fine di ridurre la distorsionedi campionamento, pesa i cluster meno pesantemente dei singoli punti. Il valore di ogni punto viene calcolato in modo tale da ridurre al minimo l’errore previsto per quel particolare punto.
Ipotizioamo
\[Y(s)=\mu(s)+Z(s)\] dove:
- \(s \in D\)
- \(E[Y(s)]=\mu(s)\)
- \(Z(s)=Y(s)-\mu(s)\), Un PPS (gaussiano) a media nulla
Si ipotizzi che il processo abbia un variogramma pari a \(2\gamma(h)\) (e di conseguenza semivariograma \(\gamma(h)\)).
Per semplicità di trascrizione in questo articolo ipotizzeremo che \(H\) sia la mtrice di distanze tra tutte le coordinate \(s\) e che \(\Gamma\) sia il semivariogramma calcolato su tutte le distanze di \(H\).
Nel si vuole prevedere il valore della variabile regionalizzata in \(s_0\) date alcune osservazioni \(Y(s_1),..., Y(s_n)\).
Nel Kriging ordinario (OK) per ottenere i coefficienti per le coordinate \(s_0\) calcolo:
\[\Gamma_0 = \left(\begin{array}{cc} \Gamma & 1\\ 1 & 0 \end{array}\right) \]
\[\gamma_0 = \left(\begin{array}{cc} \gamma(||s_1-s_0||) \\ \gamma(||s_1-s_0||)\\ ...\\ \gamma(||s_n-s_n||)\\ 1 \end{array}\right)\]
Infine stimo i coefficienti:
\[ \alpha_0=\Gamma_0^{^{-1}} \gamma_0 \]
e trovo \[\hat Y(s_0)=\sum_{i=1}^n \alpha_{0 i}\space\space s_{0 i}\]
Nel Kriging semplice (SK) si pone come assunzione che il processo stocastico sia un PPS con media uguale a \(\mu(s)\) nota quindi si applica lo stesso ragionamento fatto per il kriging ordinario ma su \(Z(s)=Y(s)-\mu(s)\) e quando poi occorre fare previsioni si aggiunge alla stima della componente di piccola scala quella di larga scala nota.
Nel kriging universale (UK) occoore prima stimare \(\mu(s)\) con una regressione polinomiale di grado \(p\)
\[Y(s)=f_1(s)+...f_p(s)+Z(s)\]
creiamo la matrice \(X\) calcolando per ogni i-esima coordinata e per ogni j-esimo grado del polinomio, la funzione \(f_{j-1-n}(s_i)\)
\[\Gamma_0 = \left(\begin{array}{cc} \Gamma & X\\ X^T & 0 \end{array}\right) \]
\[\gamma_0 = \left(\begin{array}{cc} \gamma(||s_1-s_0||) \\ \gamma(||s_1-s_0||)\\ ...\\ \gamma(||s_n-s_n||)\\ 1\\ f_1(s_0)\\ ...\\ f_p(s_0) \end{array}\right)\]
Infine stimo i coefficienti:
\[ \alpha_0=\Gamma_0^{^{-1}} \gamma_0 \]
e trovo \[\hat Y(s_0)=\sum_{i=1}^n (\alpha_{0 i}\space\space \gamma(s_{0 i}-s_i))+\sum_{i=n+1}^{n+p+1} (\alpha_{0 i}\space\space f_{i-n-1}( s_{0 i}))\]
# Importo i dati
require(geoR)## Loading required package: geoR## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.7-5.2.1 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------dat.om<-soja98; str(dat.om)## 'data.frame': 256 obs. of 10 variables:
## $ X : num 5.6 15.2 24.8 34.4 44 53.6 63.2 72.8 82.4 92 ...
## $ Y : num 3.6 1.6 4.6 5.6 1.6 1.6 2.6 4.6 3.6 2.6 ...
## $ P : num 4.4 3.7 3.5 3.9 3.6 3.7 2.9 2.2 3 2.7 ...
## $ PH : num 6.3 6.3 5.4 5.2 5.6 5.2 5.2 5.6 5.7 5.4 ...
## $ K : num 0.53 0.38 0.47 0.48 0.52 0.45 0.47 0.42 0.36 0.35 ...
## $ MO : num 47.5 48.7 51.8 57.9 59.7 ...
## $ SB : num 80.6 82.2 61.5 59.1 71.6 ...
## $ iCone: num 31 29.5 23 24.5 18.5 21.5 13 23 20 20 ...
## $ Rend : num 6.4 7.05 5.88 8.52 7.47 ...
## $ PROD : num 2.56 2.82 2.35 3.41 2.99 3.73 2.91 3.47 3 3.58 ...obj <- dat.om[,c("X","Y","MO")]
summary(obj)## X Y MO
## Min. : 1.60 Min. : 1.60 Min. :36.55
## 1st Qu.: 38.60 1st Qu.: 29.60 1st Qu.:48.07
## Median : 75.60 Median : 57.60 Median :53.61
## Mean : 75.85 Mean : 57.29 Mean :52.55
## 3rd Qu.:112.60 3rd Qu.: 85.60 3rd Qu.:56.97
## Max. :149.60 Max. :113.60 Max. :68.35attach(obj)
# Stimo vaoigramma con WLS
soja.geo <- as.geodata(obj,coords.col=1:2,data.col=3)
soja.var <- variog(soja.geo, estimator.type="classical")## variog: computing omnidirectional variogramplot(soja.var, main="variogramma empirico")
soja.var.fit <-variofit(soja.var,ini.cov.pars=c(25,35),weights="cressie",
cov.model="exponential", fix.nugget=FALSE,nugget=20)## variofit: covariance model used is exponential
## variofit: weights used: cressie
## variofit: minimisation function used: optimlines(soja.var.fit); soja.var.fit
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
## tausq sigmasq phi
## 14.8772 37.1710 56.4090
## Practical Range with cor=0.05 for asymptotic range: 168.9862
##
## variofit: minimised weighted sum of squares = 118.3311a<-soja.var.fit$cov.pars
size=30
YY<-round(seq(ceiling(min(Y)),floor(max(Y)), length=size),2)
XX<-round(seq(ceiling(min(X)),floor(max(X)), length=size),2)
griglia<-expand.grid(X=XX,Y=YY)
dim(griglia)## [1] 900 2Kriging ordinario
krg.or <- krige.conv(geodata=soja.geo, loc=griglia,
krige=krige.control(cov.pars=a,
cov.model="exponential",
nugget=soja.var.fit$nugget))## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhoodpar(cex=0.5,mfrow=c(2,2))
# 'altosinistra'
image(krg.or,xlim=c(-15,floor(max(X))),ylim=c(0,120))
legend.krige(x.leg=c(-13,-8), y.leg=c(30,100), val=krg.or$predict,
vert=TRUE, off=0.7, cex=0.9, scale.vals=pretty(krg.or$predict))
contour(krg.or,add=T,coords.data=soja.geo$coords,cex=0.3)
# 'altodestra'
persp(krg.or,xlab="longitudine",ylab="latitudine",zlab="log-catch",ticktype="detailed")
# 'bassosinistra'
image(krg.or, val=sqrt(krg.or$krige.var), main="kriging std. errors", xlim=c(min(X),floor(max(X))))
legend.krige(x.leg=c(-13,-8), y.leg=c(30,100), val=krg.or$krige.var,vert=TRUE, off=0.7, cex=0.9, scale.vals=pretty(krg.or$krige.var))
# 'bassodestra'
points.geodata(soja.geo,pt.divide="data.proportional")
par(cex=1,mfrow=c(1,1))Kriging universale
obj.var.u <- variog(soja.geo, estimator.type="classical", trend="2nd") ## variog: computing omnidirectional variogramobj.var.fit.u <- variofit(obj.var.u,ini.cov.pars=c(10,40), cov.model="exponential", fix.nugget=FALSE,nugget=15)## variofit: covariance model used is exponential
## variofit: weights used: npairs
## variofit: minimisation function used: optimplot(obj.var.u)
lines(obj.var.fit.u)
krg.uk <- krige.conv(geodata=soja.geo,
loc=griglia,
krige=krige.control(type.krige = "ok",
trend.d="2nd",
trend.l="2nd",
cov.pars=obj.var.fit.u$cov.pars,
nugget=obj.var.fit.u$nugget,
cov.model="exponential"))## krige.conv: model with mean given by a 2nd order polynomial on the coordinates
## krige.conv: Kriging performed using global neighbourhoodpar(cex=0.5,mfrow=c(2,2))
# 'altosinistra'
image(krg.uk, main="kriging universale",xlim=c(-15,floor(max(X))),ylim=c(0,120))
legend.krige(x.leg=c(-13,-8), y.leg=c(30,100), val=krg.uk$predict,
vert=TRUE, off=0.7, cex=0.9, scale.vals=pretty(krg.uk$predict))
contour(krg.uk,add=T)
# 'altodestra'
persp(krg.uk,xlab="longitudine",ylab="latitudine",zlab="log-catch",ticktype="detailed")
# 'bassosinistra'
image(krg.uk, val=sqrt(krg.uk$krige.var), main="kriging std. errors", xlim=c(min(X),floor(max(X))))
legend.krige(x.leg=c(-13,-8), y.leg=c(30,100), val=krg.uk$krige.var,vert=TRUE, off=0.7, cex=0.9, scale.vals=pretty(krg.uk$krige.var))
# 'bassodestra'
points.geodata(soja.geo,pt.divide="data.proportional")
par(cex=1,mfrow=c(1,1))Risultati
krg.uk$beta ## beta0 beta1 beta2 beta3 beta4
## 4.998546e+01 -3.453553e-02 2.682163e-02 -9.281561e-06 7.409176e-05
## beta5
## 7.958683e-04krg.or$beta## beta
## 52.76726