Nothing
####################################################
### Authors: Simone Padoan and Moreno Bevilacqua.
### Emails: simone.padoan@unibocconi.it,
### moreno.bevilacqua@uv.cl
### Institutions: Department of Decision Sciences,
### University Bocconi of Milan and
### Departamento de Estadistica
### Universidad de Valparaiso
### File name: Covariogram.r
### Description:
### This file contains a set of procedures
### to compute and plot the estimated covariance
### function and the variogram after fitting a
### random field by composite-likelihood.
### Last change: 28/03/2013.
####################################################
### Procedures are in alphabetical order.
### Compute and plot the (estimated) covariance function and the variogram
### from a fitted model obtain from the FitComposite or the WLeastSquare procedure
Covariogram <- function(fitted, lags=NULL, lagt=NULL, answer.cov=FALSE, answer.vario=FALSE,
answer.extc=FALSE, answer.range=FALSE, fix.lags=NULL, fix.lagt=NULL,
show.cov=FALSE, show.vario=FALSE, show.extc=FALSE, show.range=FALSE,
add.cov=FALSE, add.vario=FALSE, add.extc=FALSE, pract.range=95,
vario=NULL, ...)
{
result <- NULL
# define the bivariate Gaussian distribution:
vpbnorm <- function(corrmodel,lags,lagt,nuisance,numlags,numlagt,param,threshold)
{
rho=double(numlags*numlagt)
p=.C("vpbnorm", as.integer(corrmodel), as.double(lags), as.double(lagt),
as.integer(numlags), as.integer(numlagt), as.double(nuisance),
as.double(param), rho=double(numlags*numlagt), as.double(threshold),
DUP=TRUE, NAOK=TRUE)$rho
return(p)
}
# define the correlation function:
CorrelationFct <- function(corrmodel, lags, lagt, numlags, numlagt, param)
{
corr=double(numlags*numlagt)
p=.C('VectCorrelation', corr=double(numlags*numlagt), as.integer(corrmodel), as.double(lags),
as.integer(numlags), as.integer(numlagt), as.double(param),
as.double(lagt), DUP=TRUE, NAOK=TRUE)$corr
return(p)
}
# define the extremal coefficient function:
ExtremalCoeff <- function(corrmodel, lags, model, nuisance, numlags, param)
{
p=.C('ExtCoeff', as.integer(corrmodel), extc=double(numlags), as.double(lags),
as.integer(model), as.integer(numlags), as.double(nuisance),
as.double(param), DUP=TRUE, NAOK=TRUE)$extc
return(p)
}
# Pratical range in the Gaussian case:
PracRangeNorm <- function(corrmodel, lags, lagt, nuisance, numlags, numlagt, param, pract.range)
return(nuisance["nugget"]+nuisance["sill"]*CorrelationFct(corrmodel, lags, lagt, numlags, numlagt, param)-
(nuisance["nugget"]+nuisance["sill"]*(1-pract.range)))
# Pratical range in the max-stable case:
PracRangeExt <- function(corrmodel, lags, model, nuisance, numlags, param, pract.range)
return(ExtremalCoeff(corrmodel, lags, model, nuisance, numlags, param)-(2*pract.range))
isextc <- FALSE
isvario <- !is.null(vario) # is empirical variogram is passed?
if(isvario) isextc <- !is.null(vario$extcoeff)
ispatim <- fitted$numtime>1# is a space-time random field?
# START ---- check input --- #
if(!class(fitted)=='FitComposite' & !class(fitted)=='WLS')
stop("Enter an object obtained FitComposite or WLeastSquare")
if(isvario & !class(vario)=='Variogram')
stop("Enter an object obtained from the function EVariogram")
if(!is.numeric(pract.range) & answer.range)
stop("Enter a number for the parameter % of sill")
else{
if(pract.range < 0 || pract.range > 100)
stop("Entered an incorrect value for the % of sill")
else
pract.range <- pract.range / 100}
# Set the type of process:
model <- CheckModel(fitted$model)
gaussian <- model==1
binary <- model==2
maxstab <- model>=3
# determine the spatial(temporal) distances:
if(binary) {slow <- 1e-3; zero <- NA} else {slow <- 0; zero <- 0}
if(is.null(lags)) lags <- seq(slow,fitted$srange[2],length=100)
# determine the spatio-temporal distances:
if(is.null(lagt)) lagt <- seq(0,fitted$trange[2],by=1)
if(isvario & ispatim){
lags <- c(0,vario$centers)
lagt <- c(0,vario$bint)}
numlags <- length(lags)
numlagt <- length(lagt)
# check the fixed lags:
if(!is.null(fix.lags)){
if(!is.numeric(fix.lags) || fix.lags<0 || fix.lags>numlags){
cat("invalid spatial fixed lag\n")
return(result)}}
if(!is.null(fix.lagt))
if(!is.numeric(fix.lagt) || fix.lagt<0 || fix.lagt>numlagt){
cat("invalid temporal fixed lag\n")
return(result)}
# END ---- check input --- #
# set range (for practical range)
if(fitted$corrmodel=="matern"){
lower <- 1e-10
upper <- 1e20}
else{
lower <- 0
upper <- 1e100}
param <- c(fitted$fixed,fitted$param)[CorrelationParam(fitted$corrmodel)]
nuisance <- c(fitted$fixed,fitted$param)[NuisanceParam(fitted$model)]
corrmodel <- CheckCorrModel(fitted$corrmodel)
correlation <- CorrelationFct(corrmodel, lags, lagt, numlags, numlagt, param)
# Gaussian random field:
if(gaussian){
# covariance function
covariance <- nuisance["nugget"]+nuisance["sill"]*correlation
# variogram function
variogram <- nuisance["nugget"]+nuisance["sill"]*(1-correlation)
vario.main <- "Spatial variogram"
vario.ylab <- "Variogram"
if(ispatim){
dim(covariance) <- c(numlags, numlagt)
dim(variogram) <- c(numlags, numlagt)
vario.main <- "Space-time variogram"
vario.zlab <- "Variogram"
}
# compute the practical range:
if(show.range || answer.range){
Range <- uniroot(PracRangeNorm, c(lower, upper), corrmodel=corrmodel,
nuisance=nuisance, numlags=1, numlagt=1, lagt=lagt,
param=param, pract.range=pract.range)$root}}
# max-stable random field:
if(maxstab){
# covariance function
covariance <- 1-nuisance["sill"]+nuisance["sill"]*correlation
# extremal coefficient function
extrcoeff <- ExtremalCoeff(corrmodel, lags, model, nuisance, numlags, param)
# compute the practical range:
Range <- uniroot(PracRangeExt,c(lower, upper),corrmodel=corrmodel,model=model,
nuisance=nuisance,numlags=1,param=param,pract.range=pract.range)$root}
# binary random field:
if(binary){
covariance <- nuisance["nugget"]+nuisance["sill"]*correlation
p <- pnorm(nuisance["mean"])
q <- vpbnorm(corrmodel, lags, lagt, nuisance,
numlags, numlagt, param, fitted$threshold)
variogram <- log(q*(1-2*p+q)/(p-q)^2)
vario.main <- "Spatial lorelogram"
vario.ylab <- "Lorelogram"
if(ispatim){
dim(covariance) <- c(numlags, numlagt)
dim(variogram) <- c(numlags, numlagt)
vario.main <- "Space-time lorelogram"
vario.zlab <- "Lorelogram"}}
# display the covariance function
if(show.cov){
if(ispatim){# spatio-temporal case:
# build the covariance matrix:
plagt <- !is.null(fix.lags)
plags <- !is.null(fix.lagt)
numplot <- 1+plags+plagt
par(mfrow=c(1,numplot))
# temporal section
par(mai=c(.2,.2,.2,.2))
persp(lags, lagt, covariance, xlab="Distance", ylab="Time",
zlab="Covariance", ltheta=90,
shade=0.75, ticktype="detailed", phi=30,
theta=30,main="Space-time covariance",
, cex.axis=.8, cex.lab=.8,zlim=c(0,max(covariance))) #
if(plagt){
par(mai=c(.5,.5,.3,.3),mgp=c(1.6,.6,0))
plot(lagt, covariance[fix.lags,], xlab="Time",
ylab="Covariance", type="l",cex.axis=.8,cex.lab=.8,
main="Space-time cov: temporal profile",...)}
if(plags){
par(mai=c(.5,.5,.3,.3),mgp=c(1.6,.6,0))
plot(lags, covariance[,fix.lagt], xlab="Distance",
ylab="Covariance", type="l",cex.axis=.8,cex.lab=.8,
main="Space-time cov: spatial profile",...)}
}
else{# spatial case:
if(add.cov & dev.cur()!=1){
lines(lags, covariance,...)
if(show.range) abline(v=Range)}
else{
plot(lags, covariance, type='l', ylim=c(min(covariance),
max(covariance)), main="Spatial covariance",
xlab="Distance", ylab="Covariance",...)
if(show.range) abline(v=Range)}}}
# display the variogram function
if(show.vario){
if(ispatim){# spatio-temporal case:
plagt <- !is.null(fix.lags)
plags <- !is.null(fix.lagt)
nrowp <- 1
ncolp <- 1
if(isvario){
ncolp <- ncolp+1
if(plags || plagt) nrowp <- nrowp+1}
else ncolp <- ncolp+plags+plagt
par(mfrow=c(nrowp,ncolp))
sup <- 0
tup <- 0
if(isvario){
nbins <- length(vario$centers)
nbint <- length(vario$bint)
evario <- matrix(vario$variogramst,nrow=nbins,
ncol=nbint,byrow=TRUE)
evario <- rbind(c(zero,vario$variogramt),cbind(vario$variograms,evario))
evario.grid <- expand.grid(c(0,vario$centers),c(0,vario$bint))
scatterplot3d::scatterplot3d(evario.grid[,1],evario.grid[,2], c(evario),
type="h",highlight.3d=TRUE,cex.axis=.7,cex.lab=.7,
main=paste("Empirical",vario.main),xlab="Distance",
ylab="Time",zlab=vario.zlab,mar=c(2,2,2,2),mgp=c(0,0,0))
if(plagt) tup <- max(evario[fix.lags,],na.rm=TRUE)
if(plags) sup <- max(evario[,fix.lagt],na.rm=TRUE)
}
par(mai=c(.2,.2,.2,.2))
persp(lags, lagt, variogram, xlab="Distance",
ylab="Time", zlab=vario.zlab, ltheta=90,
shade=0.75, ticktype="detailed", phi=30,
theta=30,main=vario.main, cex.axis=.8,
cex.lab=.8) #zlim=c(0,max(variogram))
if(plagt){
par(mai=c(.5,.5,.3,.3),mgp=c(1.6,.6,0))
plot(lagt, variogram[fix.lags,], xlab="Time",cex.axis=.8,cex.lab=.8,
ylab=vario.ylab, type="l", ylim=c(min(0,min(variogram)),max(nuisance["nugget"]+
nuisance["sill"],tup)), main=paste(vario.ylab,": temporal profile",
sep=""),...)
if(isvario) points(lagt, evario[fix.lags,],...)}
if(plags){
par(mai=c(.5,.5,.3,.3),mgp=c(1.6,.6,0))
plot(lags, variogram[,fix.lagt], xlab="Distance",cex.axis=.8,cex.lab=.8,
ylab=vario.ylab, type="l", ylim=c(min(0,min(variogram)),max(nuisance["nugget"]+
nuisance["sill"],sup)), main=paste(vario.ylab,": spatial profile",
sep=""),...)
if(isvario) points(lags, evario[,fix.lagt],...)}}
else{# spatial case:
if(add.vario & dev.cur()!=1){
points(vario$centers, vario$variograms,...)
lines(lags, variogram,...)
if(show.range) abline(v=Range)}
else{
bnds <- range(variogram)
bnds[1] <- min(bnds[1], vario$variograms)
bnds[2] <- max(bnds[2], vario$variograms)
plot(lags, variogram, type='l', ylim=c(bnds[1], bnds[2]),
main=vario.main,xlab="Distance",
ylab=vario.ylab,...)
points(vario$centers, vario$variograms,...)
if(show.range) abline(v=Range)}}}
# display the extremal function
if(show.extc){
if(add.extc & dev.cur()!=1){
if(isvario & isextc) points(vario$centers, vario$extcoeff,...)
lines(lags, extrcoeff,...)
if(show.range) abline(v=Range)}
else{
bnds <- range(extrcoeff)
if(isvario & isextc){
bnds[1] <- min(bnds[1], vario$extcoeff)
bnds[2] <- max(bnds[2], vario$extcoeff)}
plot(lags, extrcoeff, type='l',ylim=c(bnds[1], bnds[2]),
main="Spatial extremal coefficient",ylab="Extremal coefficient",
xlab="Distance",...)
if(isvario & isextc)
points(vario$centers, vario$extcoeff,...)
if(show.range) abline(v=Range)}}
# return the estimated covariance function
if(answer.cov) result <- list(lags=lags, covariance=covariance)
# return the estimated variogram function
if(answer.vario)
if(!is.list(result)) result <- list(lags=lags, variogram=variogram)
else result$variogram <- variogram
# return the range
if(answer.range)
if(!is.list(result)) result <- list(range=Range)
else result$range <- Range
# return the estimated extremal coefficient function
if(answer.extc)
if(!is.list(result)) result <- list(lags=lags, extrcoeff=extrcoeff)
else result$extrcoeff <- extrcoeff
if(!is.null(result))
return(result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.