Simulación de campos aleatorios

###############################################################################
# Simulación de campos aleatorios
###############################################################################

require(RandomFields)
## Loading required package: RandomFields
## Loading required package: sp
## Loading required package: RandomFieldsUtils
## 
## Attaching package: 'RandomFields'
## The following objects are masked from 'package:base':
## 
##     abs, acosh, asin, asinh, atan, atan2, atanh, cos, cosh, exp,
##     expm1, floor, gamma, lgamma, log, log1p, log2, logb, max, min,
##     round, sin, sinh, sqrt, tan, tanh, trunc
rango_practico<-300

Construcción de los modelos espaciales

# En el caso exponencial el rango práctico es scale=1/3*(rango_práctico)
modeloExp<- RMtrend(mean=10)+RMexp(var=1,scale=rango_practico/3)
# En el caso gaussiano el rango práctico es scale=1/sqrt(3)*(rango_práctico)
modeloGauss<- RMtrend(mean=10)+RMgauss(var=1,scale=rango_practico/sqrt(3))

Definición de la rejilla de los datos.

inicio <- 0
final <- 1000
long <- 201
sec.x <- seq(inicio, final, length=long) 
sec.y <- seq(inicio, final, length=long)

Inicio del campo aleatorio simulado exponencial

RFoptions(seed=12345)
simuExp <- RFsimulate(modeloExp, x=sec.x, y=sec.y, grid=TRUE)
## NOTE: simulation is performed with fixed random seed 12345.
## Set 'RFoptions(seed=NA)' to make the seed arbitrary.
## New output format of RFsimulate: S4 object of class 'RFsp';
## for a bare, but faster array format use 'RFoptions(spConform=FALSE)'.
png("expo1.png")
plot(simuExp,asp=1,col=grey(0:1000/1000))
dev.off()
## png 
##   2

Campo Exponencial

Inicio del campo aleatorio simulado gausiano

RFoptions(seed=12345)
simuGauss <- RFsimulate(modeloGauss, x=sec.x, y=sec.y, grid=TRUE)
## NOTE: simulation is performed with fixed random seed 12345.
## Set 'RFoptions(seed=NA)' to make the seed arbitrary.
png("gaus1.png")
plot(simuGauss,asp=1,col=grey(0:1000/1000))
dev.off()
## png 
##   2

Campo Gaus

Semivariogramas para el modelo exponencial.

# Cálculo del variograma empírico de los datos simulados Exponencial
varioemp<-RFempiricalvariogram(data=simuExp)
plot(varioemp,model=modeloExp)

Semivariogramas para el modelo gausiano.

# Cálculo del variograma empírico de los datos simulados Exponencial
varioemp<-RFempiricalvariogram(data=simuGauss)
plot(varioemp,model=modeloGauss)

Convertir los datos simulados a bases de datos

# Convetir los datos simulado a una base datos 
datosExp<-data.frame(simuExp,coordinates(simuExp))
datosGauss<-data.frame(simuExp,coordinates(simuGauss))

Muestra aleatoria para el caso exponencial

nmuestras<-15
set.seed(NULL)
# Posición en la base de datos de las muestras seleccinadas
idmuestras<-sample(1:nrow(datosExp),nmuestras)
# Extraccion y reordenamiento de las columnas de interés.
datosExpM<-datosExp[idmuestras,c(3,2,1)]

Análisis exploratorio de la muestra tomada

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.1 (built on 2015-04-15) is now loaded
## --------------------------------------------------------------
# Tipo de datos "geodata"
datos1g<-as.geodata(datosExpM)
# Exploratorio inicial
plot(datos1g)

# Resumen de la base de datos
summary(datos1g)
## Number of data points: 15 
## 
## Coordinates summary
##     coords.x2 coords.x1
## min        70        30
## max       955       945
## 
## Distance summary
##        min        max 
##   75.66373 1114.64120 
## 
## Data summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   8.163   9.377   9.831   9.919  10.580  11.760
# Estimación del semivariograma empírico
semivariog1 <- variog(datos1g, uvec = seq(25, 750, by =25))
## variog: computing omnidirectional variogram
# 
plot(semivariog1,type = "b",xlab="Distancia", ylab="Semivarianza",lwd=2)
abline(h=1,lwd=2,col="red")
abline(h=var(datos1g$data),lwd=2,col="red",lty=2)

Detección de estructura espacial

set.seed(NULL)
mc1 <- variog.mc.env(datos1g, obj = semivariog1,nsim=10000)
## variog.env: generating 10000 simulations by permutating data values
## variog.env: computing the empirical variogram for the 10000 simulations
## variog.env: computing the envelops
plot(semivariog1, env = mc1, xlab="Distancia",ylab="Semivarianza",pch=19,lwd=2)

Muestra aleatoria para el caso gausiano.

nmuestras<-300
set.seed(NULL)
# Posición en la base de datos de las muestras seleccinadas
idmuestras<-sample(1:nrow(datosGauss),nmuestras)
# Extraccion y reordenamiento de las columnas de interés.
datosGaussM<-datosExp[idmuestras,c(3,2,1)]

Análisis exploratorio de la muestra tomada

require(geoR)
# Tipo de datos "geodata"
datos1g<-as.geodata(datosGaussM)
# Exploratorio inicial
plot(datos1g)

# Resumen de la base de datos
summary(datos1g)
## Number of data points: 300 
## 
## Coordinates summary
##     coords.x2 coords.x1
## min         0         5
## max      1000      1000
## 
## Distance summary
##      min      max 
##    5.000 1347.266 
## 
## Data summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   6.959   9.351   9.939   9.933  10.490  12.280
# Estimación del semivariograma empírico
semivariog1 <- variog(datos1g, uvec = seq(25, 750, by =25))
## variog: computing omnidirectional variogram
# 
plot(semivariog1,type = "b",xlab="Distancia", ylab="Semivarianza",lwd=2)
abline(h=1,lwd=2,col="red")
abline(h=var(datos1g$data),lwd=2,col="red",lty=2)

Detección de estructura espacial

set.seed(NULL)
mc1 <- variog.mc.env(datos1g, obj = semivariog1,nsim=10000)
## variog.env: generating 10000 simulations by permutating data values
## variog.env: computing the empirical variogram for the 10000 simulations
## variog.env: computing the envelops
plot(semivariog1, env = mc1, xlab="Distancia",ylab="Semivarianza",pch=19,lwd=2)

Esquema de simulación para detectar tamaño de muestra

# Construye una realizaciòn de un campo aleatorio
rango_practico <- 300
modeloGauss2<- RMtrend(mean=10) +
              RMgauss(var=50,scale=rango_practico/sqrt(3)) +
              RMnugget(var=0)

RFoptions(seed=21684323)
simuGauss2  <- RFsimulate(modeloGauss2, x=sec.x, y=sec.y, grid=TRUE)
## NOTE: simulation is performed with fixed random seed 21684323.
## Set 'RFoptions(seed=NA)' to make the seed arbitrary.
datosGauss2 <- data.frame(simuExp,coordinates(simuGauss2))

# El campo aletorio se toma como la poblaciòn

nmuestras <- 80
nsimul <- 10
RFoptions(seed=NULL)
set.seed(NULL)
l1 <- NULL
for (i in 1:nsimul){
  idmuestras<-sample(1:nrow(datosGauss2),nmuestras)
  datosGaussM2<-datosExp[idmuestras,c(3,2,1)]
  datosGaussGD <- as.geodata(datosGaussM2)
  semivariog1 <- variog(datosGaussGD, uvec = seq(25, 750, by =25))
  mc1 <- variog.mc.env(datosGaussGD, obj = semivariog1, nsim=10000)
  r1 <- any(mc1$v.lower > semivariog1$v)
  l1 <- c(l1,r1)
}
## variog: computing omnidirectional variogram
## variog.env: generating 10000 simulations by permutating data values
## variog.env: computing the empirical variogram for the 10000 simulations
## variog.env: computing the envelops
## variog: computing omnidirectional variogram
## variog.env: generating 10000 simulations by permutating data values
## variog.env: computing the empirical variogram for the 10000 simulations
## variog.env: computing the envelops
## variog: computing omnidirectional variogram
## variog.env: generating 10000 simulations by permutating data values
## variog.env: computing the empirical variogram for the 10000 simulations
## variog.env: computing the envelops
## variog: computing omnidirectional variogram
## variog.env: generating 10000 simulations by permutating data values
## variog.env: computing the empirical variogram for the 10000 simulations
## variog.env: computing the envelops
## variog: computing omnidirectional variogram
## variog.env: generating 10000 simulations by permutating data values
## variog.env: computing the empirical variogram for the 10000 simulations
## variog.env: computing the envelops
## variog: computing omnidirectional variogram
## variog.env: generating 10000 simulations by permutating data values
## variog.env: computing the empirical variogram for the 10000 simulations
## variog.env: computing the envelops
## variog: computing omnidirectional variogram
## variog.env: generating 10000 simulations by permutating data values
## variog.env: computing the empirical variogram for the 10000 simulations
## variog.env: computing the envelops
## variog: computing omnidirectional variogram
## variog.env: generating 10000 simulations by permutating data values
## variog.env: computing the empirical variogram for the 10000 simulations
## variog.env: computing the envelops
## variog: computing omnidirectional variogram
## variog.env: generating 10000 simulations by permutating data values
## variog.env: computing the empirical variogram for the 10000 simulations
## variog.env: computing the envelops
## variog: computing omnidirectional variogram
## variog.env: generating 10000 simulations by permutating data values
## variog.env: computing the empirical variogram for the 10000 simulations
## variog.env: computing the envelops
l1
##  [1] FALSE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE
table(l1)
## l1
## FALSE  TRUE 
##     8     2
table(l1)/sum(table(l1))
## l1
## FALSE  TRUE 
##   0.8   0.2