Cross Validated Asked by Zheyuan Li on January 7, 2022
Using grf
function from R package geoR
, I simulated 6 replicates (each with 1000 samples) of a Gaussian random field on [0, 1] x [0, 1]
, with zero mean, zero nugget, and an exponential spatial covariance $$gamma(h) = expleft(frac{h}{phi}right)sigma^2,$$ with $sigma^2 = 1$ and $phi = 0.5$. In the following figure, the blue solid curve gives the semi-variogram of this true model and the blue dashed line gives the sill (which is $sigma^2 = 1$). Then for each replicate, I
variog
(black solid dots);However, I find it difficult to understand the results / behavior as illustrated by the figure.
Is a spatial process expected to behave like this? If I examine an AR(1) process by simulating a time series and fitting it using Yule-Walker method, I can get back to the truth, i.e., true autocorrelation coefficient and variance. This is not possible for a spatial process?
Reproducible R code
set.seed(0)
phi <- 0.5
sigmasq <- 1
## unconditional simulation of Gaussian random field with exponential covariance
EXP_GRF <- grf(1000, cov.pars = c(sigmasq, phi), cov.model = "exponential",
nsim = 6)
## estimation of semi-variogram
semi_variog <- variog(EXP_GRF)
semi_variog_i <- semi_variog
semi_variog_i$v <- NULL
par(mfrow = c(2, 3), mar = c(2.5, 2.5, 1, 1))
for (i in 1:6) {
## the i-th empirocal semi-variogram
semi_variog_i$v <- semi_variog$v[, i]
## mc envelops
env <- variog.mc.env(coords = EXP_GRF$coords, data = EXP_GRF$data[, i],
obj.variog = semi_variog_i, nsim = 99, messages = FALSE)
## plotting range
MAX <- max(env$v.lower, env$v.upper, semi_variog_i$v, sigmasq)
with(semi_variog_i, plot(u, v, ylim = c(0, MAX), type = "n"))
## empirical semi-variogram
polygon(c(semi_variog_i$u, rev(semi_variog_i$u)),
c(env$v.lower, rev(env$v.upper)), col = 8, border = NA)
with(semi_variog_i, points(u, v, pch = 19))
## variance of the data
var_data <- semi_variog$var.mark[i]
abline(h = var_data, lty = 2)
## semi-variogram of the true model
lines.variomodel(EXP_GRF, col = 4, lwd = 2)
abline(h = sigmasq, lwd = 2, col = 4, lty = 2)
## fit an variogram for parameter estimation
model <- variofit(semi_variog_i, c(var_data, phi), "exponential",
fix.nugget = TRUE, nugget = 0, message = FALSE)
## fitted variogram
lines.variomodel(model)
## add estimated `sigmasq` and `phi` as title
title(sprintf("sigmasq = %.2f, phi = %.2f", model[[2]][1], model[[2]][2]))
}
I didn't have a chance to look carefully into the code, but I can already guess what you are experiencing. It may have to due with the fact that you are simulating random fields on finite domains (in a computer).
Try to increase the domain to a larger box [0,10]x[0,10]
or try to decrease the range of the variogram to something that is at least smaller than half of your domain size, say 0.1
. Then plot the variograms up to a lag that is at most half of your domain. Everything that is separated by a lag that is equal or larger to half of your domain is not useful in empirical variogram calculations (too few samples, and hence spurious oscillations).
Geostatistics theory is more general than AR models, and there are theorems related to the issue you are experiencing, but I have no time at the moment to write down a more detailed answer.
Answered by juliohm on January 7, 2022
0 Asked on January 28, 2021 by chua-s-yang
0 Asked on January 28, 2021 by visionenthusiast
0 Asked on January 28, 2021 by zqq
lasso machine learning mathematical statistics ridge regression self study
1 Asked on January 28, 2021 by nlapidot
0 Asked on January 28, 2021 by tricostume
1 Asked on January 27, 2021 by mathews24
0 Asked on January 26, 2021 by thecity2
2 Asked on January 26, 2021
0 Asked on January 26, 2021 by bonesones
0 Asked on January 26, 2021 by kplauritzen
2 Asked on January 26, 2021 by user369210
0 Asked on January 26, 2021
causality confounding difference in difference econometrics panel data
1 Asked on January 25, 2021 by kiril-e-proykov
0 Asked on January 24, 2021 by oja-niva
classification machine learning random forest recommender system regression
0 Asked on January 24, 2021 by mephisto73
0 Asked on January 23, 2021 by valjean
1 Asked on January 23, 2021 by tjt
5 Asked on January 23, 2021 by user239903
0 Asked on January 21, 2021 by pooza
Get help from others!
Recent Questions
Recent Answers
© 2022 AnswerBun.com. All rights reserved. Sites we Love: PCI Database, MenuIva, UKBizDB, Menu Kuliner, Sharing RPP