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: EmpiricalEstimators.r
### Description:
### This file contains a set of procedures in order
### to estimate the empircal covariance or the extreme
### dependence structures for a given dataset.
### Last change: 28/03/2013.
####################################################
### Procedures are in alphabetical order.
EVariogram <- function(data, coordx, coordy=NULL, coordt=NULL, cloud=FALSE, distance="Eucl",
grid=FALSE, gev=c(0,1,0), maxdist=NULL, maxtime=NULL, numbins=NULL,
replicates=1, type='variogram')
{
call <- match.call()
corrmodel <- 'gauss'
### Check the parameters given in input:
if(is.null(type))
type <- 'variogram'
# Checks if its a variogram or a madogram
if(!is.null(type) & all(type!="lorelogram", type!='variogram', type!='madogram', type!='Fmadogram'))
stop('the admitted types are: variogram, madogram (Fmadogram) or lorelogram\n')
if(type=="lorelogram" & cloud)
stop("there is not cloud version of lorelogram\n")
# Set the type of model:
if(type=='variogram'){
model <- 'Gaussian'
fname <- 'Binned_Variogram'}
if(type=="lorelogram"){
model <- "BinaryGauss"
fname <- "Binned_Lorelogram"}
if(type=="madogram" || type=="Fmadogram") model <- 'ExtGauss'
# Checks if its a spatial or spatial-temporal random field:
if(!is.null(coordt))
if(is.numeric(coordt))
if(length(coordt)>1) corrmodel <- 'gneiting'
# Checks the input:
checkinput <- CheckInput(coordx, coordy, coordt, corrmodel, data, distance, "Fitting", NULL, grid,
'None', "Frechet", maxdist, maxtime, model, NULL, 'Nelder-Mead', NULL,
replicates, NULL, NULL,NULL, 0, 'WLeastSquare', FALSE, 'SubSamp', FALSE)
# Checks if there are errors in the input:
if(!is.null(checkinput$error))
stop(checkinput$error)
### START -- Specific checks of the Empirical Variogram:
if(!is.null(cloud) & !is.logical(cloud))
stop('insert a logical value (TRUE/FALSE) for the cloud parameter\n')
if(type=='madogram' & (!is.numeric(gev) || is.null(gev) || !length(gev)==3))
stop('insert a numeric vector with the three GEV parameters\n')
if(!is.null(numbins) & !is.integer(numbins))
if(numbins < 0)
stop('insert a positive integer value for the number of bins\n')
if(is.null(numbins))
numbins <- 13
### END -- Specific checks of the Empirical Variogram
### Initialization parameters:
initparam <- InitParam(coordx, coordy, coordt, corrmodel, data,distance, "Fitting",
NULL, grid, 'None', "Frechet", maxdist,
maxtime, model, NULL, NULL, FALSE, FALSE, replicates,
NULL, NULL, NULL, 0, 'WLeastSquare', 'WLeastSquare', FALSE,
'SubSamp', FALSE, 1, 1)
# Checks if there are inconsistences:
if(!is.null(initparam$error))
stop(initparam$error)
numvario <- numbins-1
if(cloud){
numbins <- numvario <- initparam$numpairs
fname <- 'Cloud_Variogram'}
### Estimation of the empirical spatial or spatial-temporal variogram:
bins <- double(numbins) # spatial bins
moments <- double(numvario) # vector of spatial moments
lenbins <- integer(numvario) # vector of spatial bin sizes
bint <- NULL
lenbinst <- NULL
lenbint <- NULL
variogramst <- NULL
variogramt <- NULL
if(type=='madogram'){ # Trasform to GEV margins:
if(cloud) fname <- 'Cloud_Madogram' else fname <- 'Binned_Madogram'
if(gev[3]>=1)
stop('the shape parameter can not be greater or equal to 1')
initparam$data <- Dist2Dist(initparam$data, to='Gev',
loc=rep(gev[1], initparam$numcoord),
scale=rep(gev[2], initparam$numcoord),
shape=rep(gev[3], initparam$numcoord))}
if(type=='Fmadogram'){ # Transform to Uniform margins:
if(cloud) fname <- 'Cloud_Madogram' else fname <- 'Binned_Madogram'
initparam$data <- Dist2Dist(initparam$data, to='Uniform')}
if(initparam$spacetime){
numbint <- initparam$numtime-1 # number of temporal bins
bint <- double(numbint) # temporal bins
momentt <- double(numbint) # vector of temporal moments
lenbint <- integer(numbint) # vector of temporal bin sizes
numbinst <- numvario*numbint # number of spatial-temporal bins
binst <- double(numbinst) # spatial-temporal bins
momentst <- double(numbinst) # vector of spatial-temporal moments
lenbinst <- integer(numbinst) # vector of spatial-temporal bin sizes
if(cloud) fname <- 'Cloud_Variogram_st' else fname <- 'Binned_Variogram_st'
if(type=="lorelogram") fname <- "Binned_Lorelogram_st"
# Compute the spatial-temporal moments:
EV=.C(fname, bins=bins, bint=bint, as.double(initparam$data),
lenbins=lenbins, lenbinst=lenbinst, lenbint=lenbint,moments=moments, momentst=momentst,momentt=momentt,
as.integer(numbins), as.integer(numbint), DUP=TRUE, NAOK=TRUE)
bins<-EV$bins;bint<-EV$bint;
moments<-EV$moments; momentst<-EV$momentst; momentt<-EV$momentt;
lenbins<-EV$lenbins;lenbint<-EV$lenbint;lenbinst<-EV$lenbinst;
centers <- bins[1:numvario]+diff(bins)/2
if(type=="lorelogram"){
elorel <- matrix(momentst,nrow=numvario,ncol=numbint,byrow=TRUE)
elorel <- rbind(c(0,momentt),cbind(moments,elorel))
a <- rowSums(elorel)==0
b <- colSums(elorel)==0
d <- rowSums(elorel)!=0
f <- colSums(elorel)!=0
if(sum(a)){
elorel <- elorel[-which(a),]
centers <- c(0,centers)[which(d)]
centers <- centers[-1]}
if(sum(b)){
elorel <- elorel[,-which(b)]
bint <- c(0,bint)[which(f)]
bint <- bint[-1]}
elorel[elorel==0] <- NA
moments <- as.vector(elorel[,1][-1])
momentt <- as.vector(elorel[1,][-1])
momentst <- c(t(elorel[-1,-1]))
lenbins <- rep(1,length(moments))
lenbint <- rep(1,length(momentt))
lenbinst <- rep(1,length(momentst))}
indbin <- lenbins>0
indbint <- lenbint>0
indbinst <- lenbinst>0
bins <- bins[indbin]
bint <- bint[indbint]
centers <- centers[indbin]
moments <- moments[indbin]
lenbins <- lenbins[indbin]
momentt <- momentt[indbint]
lenbint <- lenbint[indbint]
momentst <- momentst[indbinst]
lenbinst <- lenbinst[indbinst]
# Computes the spatial marginal variogram:
variograms <- moments/lenbins
# Computes the temporal marginal variogram:
variogramt <- momentt/lenbint
# Computes the spatial-temporal variogram:
variogramst <- momentst/lenbinst}
else{# Computes the spatial moments
EV=.C(fname, bins=bins, as.double(initparam$data), lenbins=lenbins, moments=moments, as.integer(numbins),
DUP=TRUE, NAOK=TRUE)
bins<-EV$bins
moments<-EV$moments
lenbins<-EV$lenbins
# Computes the spatial variogram:
indbin <- lenbins>0
bins <- bins[indbin]
numbins <- length(bins)
# check if cloud or binned variogram:
if(cloud) centers <- bins else centers <- bins[1:(numbins-1)]+diff(bins)/2
moments <- moments[indbin]
lenbins <- lenbins[indbin]
variograms <- moments/lenbins}
# Start --- compute the extremal coefficient
extcoeff <- NULL
if(type == 'madogram'){ # Check if its a madogram:
if(gev[3] == 0)
extcoeff <- exp(variograms / gev[2])
else
extcoeff <- Dist2Dist(gev[1] + variograms / gamma(1 - gev[3]))
extcoeff[extcoeff>2] <- NA}
if(type == 'Fmadogram'){ # Check if its a Fmadogram:
extcoeff <- (1 + 2 * variograms) / (1 - 2 * variograms)
extcoeff[extcoeff>2] <- NA}
.C('DeleteGlobalVar', DUP=TRUE, NAOK=TRUE)
EVariogram <- list(bins=bins,
bint=bint,
cloud=cloud,
centers=centers,
extcoeff=extcoeff,
lenbins=lenbins,
lenbinst=lenbinst,
lenbint=lenbint,
srange=initparam$srange,
variograms=variograms,
variogramst=variogramst,
variogramt=variogramt,
trange=initparam$trange,
type=type)
structure(c(EVariogram, call = call), class = c("Variogram"))
}
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.