Nothing
## ======================================================================== ##
## Miguel de Carvalho ##
## Copyright (C) 2018 ##
## ------------------------------------------------------------------------ ##
## This program is free software; you can redistribute it and/or modify ##
## it under the terms of the GNU General Public License as published by ##
## the Free Software Foundation; either version 2 of the License, or ##
## (at your option) any later version. ##
## ##
## This program is distributed in the hope that it will be useful, ##
## but WITHOUT ANY WARRANTY; without even the implied warranty of ##
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ##
## GNU General Public License for more details. ##
## ##
## You should have received a copy of the GNU General Public License ##
## along with this program; if not, a copy is available at ##
## http://www.r-project.org/Licenses/ ##
## ======================================================================== ##
angdensity <- function(Y, tau = 0.95, nu, grid = seq(0.01, 0.99, length = 2^8),
method = "euclidean", raw = TRUE)
UseMethod("angdensity")
angdensity.default <- function(Y, tau = 0.95, nu,
grid = seq(0.01, 0.99, length = 2^8),
method = "euclidean", raw = TRUE) {
## Basic input validation
if (tau >= 1)
stop('tau must be smaller than 1')
if(nu < 0)
stop('nu must be positive')
if (dim(Y)[2] > 2)
stop('this function only applies to bivariate extremes')
## Convert data to common margins, if raw == TRUE
if(raw == TRUE) {
n <- dim(Y)[1]
FX <- ecdf(Y[, 1])
FY <- ecdf(Y[, 2])
x <- -1 / log(n / (n + 1) * FX(Y[, 1]))
y <- -1 / log(n / (n + 1) * FY(Y[, 2]))
}
if(raw == FALSE) {
x <- Y[, 1]
y <- Y[, 2]
}
## Compute angles and threshold radius
w <- x / (x + y)
u <- quantile(x + y, tau)
w <- w[which(x + y > u)]
## Compute angular distribution function
k <- length(w)
if(method == "euclidean") {
v <- var(w) * (k / (k - 1))
p <- 1 / k * (1 - (mean(w) - 1/2) * v^{-1} * (w - 1/2))
}
if(method == "empirical") {
lambda <- el.test(w, 1 / 2)$lambda
p <- 1 / k * (1 + (w - 1/2) * lambda)^-1
}
smooth <- function(W)
sum(p * dbeta(W, w * nu, (1 - w) * nu))
h <- sapply(grid, smooth)
## Organize and return outputs
outputs <- list(h = h, w = w, nu = nu, grid = grid, Y = Y)
class(outputs) <- "angdensity"
return(outputs)
}
plot.angdensity <- function(x, xlab = "w", ...) {
plot(x$grid, x$h, xlab = "w", ylab = "Angular Density", type = "l", ...)
}
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.