7. Covariance and Semi-variogram Models
Covariance and Semivariogram Models
For an isotropic process with variance \(\sigma^2\) and no nugget effect:
\[ \gamma(h) =
\left\{
\begin{eqnarray*}
0 &\textrm{ if }& h=0 \\
\sigma^2(1-C(h)) &\textrm{ if }& h >0
\end{eqnarray*}
\right.
\]
Matern Class of Covariance Functions (Matern (1986))
\[C(h) = \sigma^2 \frac{1}{\Gamma(\nu)}
\left(\frac{\theta h}{ 2} \right) 2 K_{\nu}(\theta h), \nu > 0,
\theta > 0\] where \(K_{\nu}\) is the modified Bessel function
of the second kind of order \(\nu\)
\[K_{\nu}(t) \approx \frac{\Gamma(\nu)}{2}
\left(\frac{t}{2} \right)^{-\nu} \textrm{ as } t \to 0.\]
Matern Class of Covariance Functions (Matern (1986))
sigma <- 2
nu <- 8 #test what happens as nu is increased
theta <- 10 # test what happens as you change theta
curve((sigma^2)*(1/gamma(nu))*(theta*x/2)*2*besselY(theta*x,nu), from = 0, to = 250)
Matern Class of Covariance Functions (Matern (1986))
As \(\nu \to \infty\):
Gaussian model \[C(h) = \sigma^2
e^{-\theta h^2} = \sigma^2 e^{-3 \frac{h^2}{\alpha^2}}\] where
\(\alpha\) is the practical range (the
distance at which the correlations have decreased to \(\approx 0.05\) or less (\(e^{-3} = 0.04978\)))
Matern Class of Covariance Functions (Matern (1986))
sigma <- 10
alpha <- 100
curve((sigma^2)*exp(-3*(x^2)/(alpha^2)), from = 0, to = 250)
Matern Class of Covariance Functions (Matern (1986))
For \(\nu = 0.5\) :
exponential model \[C(h) = \sigma^2
e^{-\theta h} = \sigma^2 e^{-3 \frac{h}{\alpha}}\]
sigma <- 10
alpha <- 100
curve((sigma^2)*exp(-3*(x)/(alpha)), from = 0, to = 250)
Matern Class of Covariance Functions (Matern (1986))
Whittle believes these two models are only valid for \(\mathbb{R}\), and thus suggests the
following model for \(\mathbb{R}^2\)
(which is however more complicated to model):
For \(\nu =1\): Whittle
(1954) \[C(h) = \sigma^2 \theta h
K_1(\theta h)\]
Matern Class of Covariance Functions (Matern (1986))
sigma <- 2
theta <- 2 # test what happens as you change theta
curve((sigma^2)*(theta)*(x)*besselY(theta*x,1), from = 0, to = 250)
Basic Models when not Second-Order Stationary (but still
isotropic)
Power model: \(\gamma(h) = \theta
h^{\lambda}, \theta \ge 0, 0 \le \lambda < 2\). When \(\lambda = 1\), a linear semivariogram is
produced.
theta <- 10
lambda <- 1.5
curve(theta*x^(lambda), from = 0, to = 200)
Estimating the Semivariogram (the empirical semivariogram)
Matheron (1962): An unbiased estimator of \(\gamma(h)\) \[\hat{\gamma}(h) = \frac{1}{2|N(h)|} \sum_{N(h)}
(Z(s_i) - Z(s_j))^2\]
- At least 30 pairs of locations should be available at each \(h\), thus group into \(h \pm \epsilon\).
Estimating the Semivariogram (the empirical semivariogram)
## [1] 181072 181025 181165 181298 181307 181390 181165 181027 181060 181232
## [11] 181191 181032 180874 180969 181011 180830 180763 180694 180625 180555
## [21] 180642 180704 180704 181153 181147 181167 181008 180973 180916 181352
## [31] 181133 180878 180829 180954 180956 180710 180632 180530 180478 180383
## [41] 180494 180561 180451 180410 180355 180292 180283 180282 180270 180199
## [51] 180135 180237 180103 179973 179826 179687 179792 179902 180100 179604
## [61] 179526 179495 179489 179414 179334 179255 179470 179692 179852 179140
## [71] 179128 179065 179007 179110 179032 179095 179058 178810 178912 178981
## [81] 179076 180151 179211 181118 179474 179559 179022 178953 178875 178803
## [91] 179029 178605 178701 179547 179301 179405 179462 179293 179180 179206
## [101] 179618 179782 179980 180067 180162 180451 180328 180276 180114 179881
## [111] 179774 179657 179731 179717 179446 179524 179644 180321 180162 180029
## [121] 179797 179642 179849 180265 180107 180462 180478 180347 180862 180700
## [131] 180201 180173 180923 180467 179917 179822 179991 179120 179034 179085
## [141] 179236 179456 179550 179445 179337 179245 179024 178786 179135 179030
## [151] 179184 179085 178875 179466 180627
## [1] 333611 333558 333537 333484 333330 333260 333370 333363 333231 333168
## [11] 333115 333031 333339 333252 333161 333246 333104 332972 332847 332707
## [21] 332708 332717 332664 332925 332823 332778 332777 332687 332753 332946
## [31] 332570 332489 332450 332399 332318 332330 332445 332538 332578 332476
## [41] 332330 332193 332175 332031 332299 332157 332014 331861 331707 331591
## [51] 331552 332351 332297 332255 332217 332161 332035 332113 332213 332059
## [61] 331936 331770 331633 331494 331366 331264 331125 330933 330801 330955
## [71] 330867 330864 330727 330758 330645 330636 330510 330666 330779 330924
## [81] 331005 330353 331175 333214 331304 331423 330873 330742 330516 330349
## [91] 330394 330406 330557 330245 330179 330567 330766 330797 330710 330398
## [101] 330458 330540 330773 331185 331387 331473 331158 330963 330803 330912
## [111] 330921 331150 331245 331441 331422 331565 331730 330366 331837 331720
## [121] 331919 331955 332142 332297 332101 331947 331822 331700 333116 332882
## [131] 331160 331923 332874 331694 331325 331242 331069 330578 330561 330433
## [141] 330046 330072 329940 329807 329870 329714 329733 329822 329890 330082
## [151] 330182 330292 330311 330381 330190
plot(x = meuse$x, y = meuse$y)
Estimating the Semivariogram (the empirical semivariogram)
# no trend:
coordinates(meuse) = ~x+y
vv1 <- gstat::variogram(log(zinc)~1, meuse)
plot(vv1)
vv11 <- gstat::variogram(zinc~1, meuse)
plot(vv11)
Estimating the Semivariogram (the empirical semivariogram)
# residual variogram w.r.t. a linear trend:
vv2 <- variogram(log(zinc)~x+y, meuse)
plot(vv2)
Estimating the Semivariogram (the empirical semivariogram)
# directional variogram:
vv3 <- variogram(log(zinc)~x+y, meuse, alpha=c(0,45,90,135))
plot(vv3)
Estimating the Semivariogram (the empirical semivariogram)
vv4 <- variogram(log(zinc)~1, meuse, width=50, cutoff=1100)
plot(vv4)
Estimating the Semivariogram (the empirical semivariogram)
# GLS residual variogram:
v.fit <- fit.variogram(vv2, vgm(1, "Sph", range = 700, nugget = 1)) #Try "Mat", "Exp"
v.fit
## model psill range
## 1 Nug 0.08234213 0.000
## 2 Sph 0.38866509 1098.571
Estimating the Semivariogram (the empirical semivariogram)
## np dist gamma dir.hor dir.ver id
## 1 57 79.29244 0.1060834 0 0 var1
## 2 299 163.97367 0.1829983 0 0 var1
## 3 419 267.36483 0.2264256 0 0 var1
## 4 457 372.73542 0.2847192 0 0 var1
## 5 547 478.47670 0.3162418 0 0 var1
## 6 533 585.34058 0.3571578 0 0 var1
## 7 574 693.14526 0.3701742 0 0 var1
## 8 564 796.18365 0.4201289 0 0 var1
## 9 589 903.14650 0.4216983 0 0 var1
## 10 543 1011.29177 0.4772549 0 0 var1
## 11 500 1117.86235 0.5075874 0 0 var1
## 12 477 1221.32810 0.4617632 0 0 var1
## 13 452 1329.16407 0.5512305 0 0 var1
## 14 457 1437.25620 0.4352155 0 0 var1
## 15 415 1543.20248 0.4556815 0 0 var1
Estimating the Semivariogram (the empirical semivariogram)
g <- gstat(NULL, "log-zinc", log(zinc)~x+y, meuse, model=v.fit, set = set)
vg <- variogram(g)
plot(vg)
Estimating the Semivariogram (the empirical semivariogram)
