La variable propuesta tiene una función \(\rho(u)\) Gaussiana, por lo tanto \(\kappa \to \infty\), se propone trabajar con \(\kappa = \{1.5,3,5,9\}\), con un rango práctico fijo de \(120\).
require(geoR)
# Función de Matérn para igualar a 0
# Nivel de rango práctico rho = 0.05
fM <- function(x,kappa, rho = 0.05, phi = 1)
cov.spatial(x,cov.pars=c(1,phi),kappa=kappa,cov.model="matern")-rho
# Encontrar el factor multiplicador de phi
r1<-uniroot(fM,c(2,10),kappa=1.5)$root
r2<-uniroot(fM,c(2,10),kappa=3)$root
r3<-uniroot(fM,c(2,10),kappa=5)$root
r4<-uniroot(fM,c(2,20),kappa=9)$root
# Rango práctico variable propuesta
RPract<-120
# Phi teórico de la distribución gaussiana
Phi<-RPract/sqrt(3)
fGaus<-function(u,p) exp(-(u/p)^2)
# Valores phi adecuados modelo matérn
phi1<-RPract/r1
phi2<-RPract/r2
phi3<-RPract/r3
phi4<-RPract/r4
# Gráfico Función Gaussiana propuesta vs
# Estimación Función Matérn.
x <- seq(0,200,l=1000)
#Gráfico Función Gaussiana
plot(x,fGaus(x,Phi),type='l',ylim=c(0,1),col=2,lwd=2,
xlab='Distancias',ylab='Correlación',
main = "Función Matérn vs Función Gaussiana")
lines(x,cov.spatial(x,cov.model="matern",kappa=1.5,cov.pars=c(1,phi1)),lty=2,lwd=2)
lines(x,cov.spatial(x,cov.model="matern",kappa=3,cov.pars=c(1,phi2)),lty=3,lwd=2)
lines(x,cov.spatial(x,cov.model="matern",kappa=5,cov.pars=c(1,phi3)),lty=4,lwd=2)
lines(x,cov.spatial(x,cov.model="matern",kappa=9,cov.pars=c(1,phi4)),lty=5,lwd=2)
# Rangos prácticos
abline(h=0.05,lty=3)
abline(v=120,lty=3,col=2)
legend(122,1,
legend=c("Función Gaussiana",
"kappa= 1.5, phi= 25.3",
"kappa= 3.0, phi= 18.7",
"kappa= 5.0, phi= 14.8",
"kappa= 9.0, phi= 11.3"),
lty=c(1,2,3,4,5),lwd=2,col=c(2,1,1,1,1))
En resúmen el modelo Matérn con los parámetros \(\kappa\) y \(\phi\) de la variable propuesta se presentan en la siguiente tabla:
| Kappa = 1.5 | Kappa = 3 | Kappa = 5 | Kappa = 9 | |
|---|---|---|---|---|
| Phi | 25.3 | 18.7 | 14.83 | 11.26 |
La función Matérn que estima mejor la función Gaussiana dado el rango práctico, es aquella con \(\kappa = 9\) y \(\phi = 11.25667\).
require(RandomFields)
xx<-seq(0,200,l=500)
ModelMat<-RMmatern(nu=9,scale=phi4,var=1.2)+RMtrend(mean=5)
RFoptions(seed = 8)
SimuMat <- RFsimulate(ModelMat,x=xx,y=xx)
plot(SimuMat, asp=1)
VarGramEmp<-RFempiricalvariogram(data=SimuMat)
plot(VarGramEmp,model=ModelMat)