library(GeoModels)
###############################################################
############ Examples of spatial Gaussian random fieldss ################
###############################################################
# Define the spatial-coordinates of the points:
set.seed(3)
N=300 # number of location sites
x <- runif(N, 0, 1)
y <- runif(N, 0, 1)
coords <- cbind(x,y)
# Define spatial matrix covariates and regression parameters
X=cbind(rep(1,N),runif(N))
mean <- 0.2
mean1 <- -0.5
# Set the covariance model's parameters:
corrmodel <- "Matern"
sill <- 1
nugget <- 0
scale <- 0.2/3
smooth=0.5
param<-list(mean=mean,mean1=mean1,sill=sill,nugget=nugget,scale=scale,smooth=smooth)
# Simulation of the spatial Gaussian random fields:
data <- GeoSim(coordx=coords,corrmodel=corrmodel, param=param,X=X)$data
################################################################
###
### Example 0. Maximum independence composite likelihood fitting of
### a Gaussian random fields (no dependence parameters)
###
###############################################################
# setting starting parameters to be estimated
start<-list(mean=mean,mean1=mean1,sill=sill)
fit1 <- GeoFit(data=data,coordx=coords,likelihood="Marginal",
type="Independence", start=start,X=X)
print(fit1)
################################################################
###
### Example 1. Maximum conditional pairwise likelihood fitting of
### a Gaussian random fields using Nelder-Mead
###
###############################################################
# setting fixed and starting parameters to be estimated
fixed<-list(nugget=nugget,smooth=smooth)
start<-list(mean=mean,mean1=mean1,scale=scale,sill=sill)
fit1 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel,
neighb=3,likelihood="Conditional",optimizer="Nelder-Mead",
type="Pairwise", start=start,fixed=fixed,X=X)
print(fit1)
################################################################
###
### Example 2. Standard Maximum likelihood fitting of
### a Gaussian random fields using nlminb
###
###############################################################
# Define the spatial-coordinates of the points:
set.seed(3)
N=250 # number of location sites
x <- runif(N, 0, 1)
y <- runif(N, 0, 1)
coords <- cbind(x,y)
param<-list(mean=mean,sill=sill,nugget=nugget,scale=scale,smooth=smooth)
data <- GeoSim(coordx=coords,corrmodel=corrmodel, param=param)$data
# setting fixed and parameters to be estimated
fixed<-list(nugget=nugget,smooth=smooth)
start<-list(mean=mean,scale=scale,sill=sill)
I=Inf
lower<-list(mean=-I,scale=0,sill=0)
upper<-list(mean=I,scale=I,sill=I)
fit2 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel,
optimizer="nlminb",upper=upper,lower=lower,
likelihood="Full",type="Standard",
start=start,fixed=fixed)
print(fit2)
###############################################################
############ Examples of spatial non-Gaussian random fieldss #############
###############################################################
################################################################
###
### Example 3. Maximum pairwise likelihood fitting of a Weibull random fields
### with Generalized Wendland correlation with Nelder-Mead
###
###############################################################
set.seed(524)
# Define the spatial-coordinates of the points:
N=300
x <- runif(N, 0, 1)
y <- runif(N, 0, 1)
coords <- cbind(x,y)
X=cbind(rep(1,N),runif(N))
mean=1; mean1=2 # regression parameters
nugget=0
shape=2
scale=0.2
smooth=0
model="Weibull"
corrmodel="GenWend"
param=list(mean=mean,mean1=mean1,scale=scale,
shape=shape,nugget=nugget,power2=4,smooth=smooth)
# Simulation of a non stationary weibull random fields:
data <- GeoSim(coordx=coords, corrmodel=corrmodel,model=model,X=X,
param=param)$data
fixed<-list(nugget=nugget,power2=4,smooth=smooth)
start<-list(mean=mean,mean1=mean1,scale=scale,shape=shape)
# Maximum independence likelihood:
fit <- GeoFit(data=data, coordx=coords, X=X,
likelihood="Marginal",type="Independence", corrmodel=corrmodel,
,model=model, start=start, fixed=fixed)
print(unlist(fit$param))
## estimating dependence parameter fixing vector mean parameter
Xb=as.numeric(X %*% c(mean,mean1))
fixed<-list(nugget=nugget,power2=4,smooth=smooth,mean=Xb)
start<-list(scale=scale,shape=shape)
# Maximum conditional composite-likelihood fitting of the random fields:
fit1 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,
neighb=3,likelihood="Conditional",type="Pairwise",
optimizer="Nelder-Mead",
start=start,fixed=fixed)
print(unlist(fit1$param))
### joint estimation of the dependence parameter and mean parameters
fixed<-list(nugget=nugget,power2=4,smooth=smooth)
start<-list(mean=mean,mean1=mean1,scale=scale,shape=shape)
fit2 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,
neighb=3,likelihood="Conditional",type="Pairwise",X=X,
optimizer="Nelder-Mead",
start=start,fixed=fixed)
print(unlist(fit2$param))
################################################################
###
### Example 4. Maximum pairwise likelihood fitting of
### a Skew-Gaussian spatial random fields with Wendland correlation
###
###############################################################
set.seed(261)
model="SkewGaussian"
# Define the spatial-coordinates of the points:
x <- runif(500, 0, 1)
y <- runif(500, 0, 1)
coords <- cbind(x,y)
corrmodel="Wend0"
mean=0;nugget=0
sill=1
skew=-4.5
power2=4
c_supp=0.2
# model parameters
param=list(power2=power2,skew=skew,
mean=mean,sill=sill,scale=c_supp,nugget=nugget)
data <- GeoSim(coordx=coords, corrmodel=corrmodel,model=model, param=param)$data
plot(density(data))
fixed=list(power2=power2,nugget=nugget)
start=list(scale=c_supp,skew=skew,mean=mean,sill=sill)
lower=list(scale=0,skew=-I,mean=-I,sill=0)
upper=list(scale=I,skew=I,mean=I,sill=I)
# Maximum marginal pairwise likelihood:
fit1 <- GeoFit(data=data,coordx=coords,corrmodel=corrmodel, model=model,
neighb=3,likelihood="Marginal",type="Pairwise",
optimizer="bobyqa",lower=lower,upper=upper,
start=start,fixed=fixed)
print(unlist(fit1$param))
################################################################
###
### Example 5. Maximum pairwise likelihood fitting of
### a Bernoulli random fields with exponential correlation
###
###############################################################
set.seed(422)
N=250
x <- runif(N, 0, 1)
y <- runif(N, 0, 1)
coords <- cbind(x,y)
mean=0.1; mean1=0.8; mean2=-0.5 # regression parameters
X=cbind(rep(1,N),runif(N),runif(N)) # marix covariates
corrmodel <- "Wend0"
param=list(mean=mean,mean1=mean1,mean2=mean2,nugget=0,scale=0.2,power2=4)
# Simulation of the spatial Binomial-Gaussian random fields:
data <- GeoSim(coordx=coords, corrmodel=corrmodel, model="Binomial", n=1,X=X,
param=param)$data
## estimating the marginal parameters using independence cl
fixed <- list(power2=4,scale=0.2,nugget=0)
start <- list(mean=mean,mean1=mean1,mean2=mean2)
# Maximum independence likelihood:
fit <- GeoFit(data=data, coordx=coords,n=1, X=X,
likelihood="Marginal",type="Independence",corrmodel=corrmodel,
,model="Binomial", start=start, fixed=fixed)
print(fit)
## estimating dependence parameter fixing vector mean parameter
Xb=as.numeric(X %*% c(mean,mean1,mean2))
fixed <- list(nugget=0,power2=4,mean=Xb)
start <- list(scale=0.2)
# Maximum conditional pairwise likelihood:
fit1 <- GeoFit(data=data, coordx=coords, corrmodel=corrmodel,n=1,
likelihood="Conditional",type="Pairwise", neighb=3
,model="Binomial", start=start, fixed=fixed)
print(fit1)
## estimating jointly marginal and dependnce parameters
fixed <- list(nugget=0,power2=4)
start <- list(mean=mean,mean1=mean1,mean2=mean2,scale=0.2)
# Maximum conditional pairwise likelihood:
fit2 <- GeoFit(data=data, coordx=coords, corrmodel=corrmodel,n=1, X=X,
likelihood="Conditional",type="Pairwise", neighb=3
,model="Binomial", start=start, fixed=fixed)
print(fit2)
###############################################################
######### Examples of Gaussian spatio-temporal random fieldss ###########
###############################################################
set.seed(52)
# Define the temporal sequence:
time <- seq(1, 9, 1)
# Define the spatial-coordinates of the points:
x <- runif(20, 0, 1)
y <- runif(20, 0, 1)
coords=cbind(x,y)
# Set the covariance model's parameters:
scale_s=0.2/3;scale_t=1
smooth_s=0.5;smooth_t=0.5
sill=1
nugget=0
mean=0
param<-list(mean=0,scale_s=scale_s,scale_t=scale_t,
smooth_t=smooth_t, smooth_s=smooth_s ,sill=sill,nugget=nugget)
# Simulation of the spatial-temporal Gaussian random fields:
data <- GeoSim(coordx=coords,coordt=time,corrmodel="Matern_Matern",
param=param)$data
################################################################
###
### Example 6. Maximum pairwise likelihood fitting of a
### space time Gaussian random fields with double-exponential correlation
###
###############################################################
# Fixed parameters
fixed<-list(nugget=nugget,smooth_s=smooth_s,smooth_t=smooth_t)
# Starting value for the estimated parameters
start<-list(mean=mean,scale_s=scale_s,scale_t=scale_t,sill=sill)
# Maximum composite-likelihood fitting of the random fields:
fit <- GeoFit(data=data,coordx=coords,coordt=time,
corrmodel="Matern_Matern",maxtime=1,neighb=3,
likelihood="Marginal",type="Pairwise",
start=start,fixed=fixed)
print(fit)
###############################################################
######### Examples of a bivariate Gaussian random fields ###########
###############################################################
################################################################
### Example 7. Maximum pairwise likelihood fitting of a
### bivariate Gaussian random fields with separable Bivariate matern
### (cross) correlation model
###############################################################
# Define the spatial-coordinates of the points:
set.seed(89)
x <- runif(300, 0, 1)
y <- runif(300, 0, 1)
coords=cbind(x,y)
# parameters
param=list(mean_1=0,mean_2=0,scale=0.1,smooth=0.5,sill_1=1,sill_2=1,
nugget_1=0,nugget_2=0,pcol=0.2)
# Simulation of a spatial bivariate Gaussian random fields:
data <- GeoSim(coordx=coords, corrmodel="Bi_Matern_sep",
param=param)$data
# selecting fixed and estimated parameters
fixed=list(mean_1=0,mean_2=0,nugget_1=0,nugget_2=0,smooth=0.5)
start=list(sill_1=var(data[1,]),sill_2=var(data[2,]),
scale=0.1,pcol=cor(data[1,],data[2,]))
# Maximum marginal pairwise likelihood
fitcl<- GeoFit(data=data, coordx=coords, corrmodel="Bi_Matern_sep",
likelihood="Marginal",type="Pairwise",
start=start,fixed=fixed,
neighb=3)
print(fitcl)
Run the code above in your browser using DataLab