Nothing
R/bias.R
In mev: Modelling of Extreme Values
Defines functions
gpd.abias
gev.abias
gev.Nyr
.gev.postpred
gev.bcor
gpd.bcor
gev.Fscore
gpd.Fscore
gpd.bias
gev.bias
Documented in
gev.abias
gev.bcor
gev.bias
gev.Fscore
gev.Nyr
.gev.postpred
gpd.abias
gpd.bcor
gpd.bias
gpd.Fscore
#' Cox-Snell first order bias for the GEV distribution
#'
#' Bias vector for the GEV distribution based on an \code{n} sample.
#' Due to numerical instability, values of the information matrix and the bias
#' are linearly interpolated when the value of the shape parameter is close to zero.
#' @inheritParams gev
#' @export
#' @keywords internal
#' @seealso \code{\link{gev}}
gev.bias <- function(par, n) {
if (length(n) > 1) {
stop("Invalid argument for sample size")
}
if (length(par) != 3) {
stop("Invalid argument for parameter vector, must be a vector of length 3.")
}
sigma <- par[2]
xi <- par[3]
if (xi < -1/3) {
stop("Cox-Snell correction only valid if the shape is greater than -1/3")
}
zeta3 <- 1.20205690315959
zeta5 <- 1.03692775514337 #gsl::zeta(5)
k111 <- ((1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi))/sigma^3
k11.2 <- -2 * (xi + 1)^2 * gamma(2 * xi + 1)/sigma^3
k11.3 <- 2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1)/sigma^2 + 2 * (xi + 1) * gamma(2 * xi + 1)/sigma^2
k33.2 <- 0
euler_gamma <- -psigamma(1)
# Numerical tolerance 1e-10 (all.equal has this)
if (abs(xi) > 1e-10) {
k112 <- (xi + 1) * (gamma(2 * xi + 2) - (xi + 1) * (4 * xi + 1) * gamma(3 * xi + 1))/(sigma^3 * xi)
k12.2 <- 2 * ((xi + 1)^2 * gamma(2 * xi + 1) - gamma(xi + 2))/(sigma^3 * xi)
} else {
k112 <- (euler_gamma - 3)/sigma^3
k12.2 <- -2 * (euler_gamma - 1)/sigma^3
}
# Numerical tolerance 1e-6 If function is not too numerically unstable
if (abs(xi) > 1e-06) {
k113 <- (1 + xi) * ((1 + xi) * (1 + 4 * xi) * gamma(1 + 3 * xi) - gamma(1 + 2 * xi) * (1 + 2 * xi * (2 + xi) + xi * (1 + 2 *
xi) * psigamma(2 + 2 * xi, 0)))/(sigma^2 * xi^2)
k122 <- ((1 - xi) * gamma(2 + xi) - gamma(3 + 2 * xi) + (1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi))/(sigma^3 * xi^2)
k12.3 <- -(2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1) + 2 * (xi + 1) * gamma(2 * xi + 1) - psigamma(xi + 2) *
gamma(xi + 2))/(sigma^2 * xi) + ((xi + 1)^2 * gamma(2 * xi + 1) - gamma(xi + 2))/(sigma^2 * xi^2)
k13.2 <- -(((xi + 1)^2 * gamma(2 * xi + 1) - xi * ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2))/(sigma^2 * xi^2))
k22.2 <- -2 * ((xi + 1)^2 * gamma(2 * xi + 1) - 2 * gamma(xi + 2) + 1)/(sigma^3 * xi^2)
} else {
# Compute the approximation for xi=0, (limit)
k113_0 <- -1/12 * (36 * euler_gamma - 6 * euler_gamma^2 - pi^2 - 24)/sigma^2
k122_0 <- -1/6 * (36 * euler_gamma - 6 * euler_gamma^2 - pi^2 - 24)/sigma^3
k12.3_0 <- (6 * euler_gamma - 3 * euler_gamma^2 - 1/2 * pi^2 - 2)/(2 * sigma^2)
k13.2_0 <- 1/12 * (12 * euler_gamma - 6 * euler_gamma^2 - pi^2)/sigma^2
k22.2_0 <- 1/3 * (12 * euler_gamma - 6 * euler_gamma^2 - pi^2 - 6)/sigma^3
if (isTRUE(all.equal(xi, 0))) {
# if xi ==0 with roughly precision 1e-8, then set value to this limit
k113 <- k113_0
k122 <- k122_0
k12.3 <- k12.3_0
k13.2 <- k13.2_0
k22.2 <- k22.2_0
} else {
# if xi !=0, but numerical breakdown,
xit <- sign(xi) * 1e-06
k113_l <- sapply(xit, function(xi) {
(1 + xi) * ((1 + xi) * (1 + 4 * xi) * gamma(1 + 3 * xi) - gamma(1 + 2 * xi) * (1 + 2 * xi * (2 + xi) + xi * (1 + 2 *
xi) * psigamma(2 + 2 * xi, 0)))/(sigma^2 * xi^2)
})
k122_l <- sapply(xit, function(xi) {
((1 - xi) * gamma(2 + xi) - gamma(3 + 2 * xi) + (1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi))/(sigma^3 * xi^2)
})
k12.3_l <- sapply(xit, function(xi) {
-(2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1) + 2 * (xi + 1) * gamma(2 * xi + 1) - psigamma(xi + 2) *
gamma(xi + 2))/(sigma^2 * xi) + ((xi + 1)^2 * gamma(2 * xi + 1) - gamma(xi + 2))/(sigma^2 * xi^2)
})
k13.2_l <- sapply(xit, function(xi) {
-(((xi + 1)^2 * gamma(2 * xi + 1) - xi * ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2))/(sigma^2 * xi^2))
})
k22.2_l <- sapply(xit, function(xi) {
-2 * ((xi + 1)^2 * gamma(2 * xi + 1) - 2 * gamma(xi + 2) + 1)/(sigma^3 * xi^2)
})
# use linear interpolation between two values 0 and 1e-6
k113 <- approx(x = c(0, xit), y = c(k113_0, k113_l), xout = xi)$y
k122 <- approx(x = c(0, xit), y = c(k122_0, k122_l), xout = xi)$y
k12.3 <- approx(x = c(0, xit), y = c(k12.3_0, k12.3_l), xout = xi)$y
k13.2 <- approx(x = c(0, xit), y = c(k13.2_0, k13.2_l), xout = xi)$y
k22.2 <- approx(x = c(0, xit), y = c(k22.2_0, k22.2_l), xout = xi)$y
}
}
# Numerical tolerance 1e-4 If function is not too numerically unstable
if (abs(xi) > 1e-04) {
k123 <- ((-gamma(2 + xi)) * (1 + 2 * xi + xi * psigamma(2 + xi, 0)) + (1 + xi) * ((-(1 + 5 * xi + 4 * xi^2)) * gamma(1 + 3 *
xi) + gamma(1 + 2 * xi) * (2 + 7 * xi + 3 * xi^2 + xi * (1 + 2 * xi) * psigamma(2 + 2 * xi, 0))))/(sigma^2 * xi^3)
k222 <- (1 - 3 * xi + 3 * (xi - 1) * gamma(2 + xi) + 1.5 * gamma(3 + 2 * xi) - (1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi))/(sigma^3 *
xi^3)
k13.3 <- -(-(2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1) - xi * ((xi + 1)/xi + psigamma(xi + 1)) * psigamma(xi +
2) * gamma(xi + 2) + xi * ((xi + 1)/xi^2 - 1/xi - psigamma(xi + 1, 1)) * gamma(xi + 2) + 2 * (xi + 1) * gamma(2 * xi +
1) - ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2))/(sigma * xi^2) + 2 * ((xi + 1)^2 * gamma(2 * xi + 1) - xi * ((xi +
1)/xi + psigamma(xi + 1)) * gamma(xi + 2))/(sigma * xi^3))
k22.3 <- 2 * ((xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1) + (xi + 1) * gamma(2 * xi + 1) - psigamma(xi + 2) * gamma(xi +
2))/(sigma^2 * xi^2) - 2 * ((xi + 1)^2 * gamma(2 * xi + 1) - 2 * gamma(xi + 2) + 1)/(sigma^2 * xi^3)
k23.2 <- -(-((digamma(1)) + (xi + 1)^2 * gamma(2 * xi + 1)/xi - ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2) - (gamma(xi +
2) - 1)/xi + 1)/(sigma^2 * xi^2))
} else {
# Compute the approximation for xi=0, (limit)
k123_0 <- 1/12 * (60 * euler_gamma + 6 * euler_gamma^3 - euler_gamma * pi^2 + 4 * pi^2 * (euler_gamma - 1) - 48 * euler_gamma^2 -
4 * pi^2 + 12 * zeta3 - 12)/sigma^2
k222_0 <- 1/4 * (48 * euler_gamma + 4 * euler_gamma^3 + 9 * euler_gamma * pi^2 - 4 * pi^2 * (2 * euler_gamma - 3) + pi^2 *
(euler_gamma - 1) - 36 * euler_gamma^2 - 17 * pi^2 + 8 * zeta3 - 16)/sigma^3
k13.3_0 <- (-6 * euler_gamma - 4 * euler_gamma^3 - 7/2 * euler_gamma * pi^2 + 3/2 * pi^2 * (euler_gamma - 1) + 12 * euler_gamma^2 +
7/2 * pi^2 - 8 * zeta3)/(6 * sigma)
k22.3_0 <- -1/6 * (12 * euler_gamma + 6 * euler_gamma^3 + 4 * euler_gamma * pi^2 - pi^2 * (euler_gamma - 1) - 18 * euler_gamma^2 -
4 * pi^2 + 12 * zeta3)/sigma^2
k23.2_0 <- -1/12 * (12 * euler_gamma + 6 * euler_gamma^3 - euler_gamma * pi^2 + 4 * pi^2 * (euler_gamma - 1) - 18 * euler_gamma^2 +
pi^2 + 12 * zeta3)/sigma^2
if (isTRUE(all.equal(xi, 0))) {
# if xi numerically zero, use the latter
k123 <- k123_0
k222 <- k222_0
k13.3 <- k13.3_0
k22.3 <- k22.3_0
k23.2 <- k23.2_0
} else {
# else if less than tol, but not zero, interpolate linearly
xit <- sign(xi) * 1e-04
k123_l <- sapply(xit, function(xi) {
((-gamma(2 + xi)) * (1 + 2 * xi + xi * psigamma(2 + xi, 0)) + (1 + xi) * ((-(1 + 5 * xi + 4 * xi^2)) * gamma(1 + 3 *
xi) + gamma(1 + 2 * xi) * (2 + 7 * xi + 3 * xi^2 + xi * (1 + 2 * xi) * psigamma(2 + 2 * xi, 0))))/(sigma^2 * xi^3)
})
k222_l <- sapply(xit, function(xi) {
(1 - 3 * xi + 3 * (xi - 1) * gamma(2 + xi) + 1.5 * gamma(3 + 2 * xi) - (1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi))/(sigma^3 *
xi^3)
})
k13.3_l <- sapply(xit, function(xi) {
-(-(2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1) - xi * ((xi + 1)/xi + psigamma(xi + 1)) * psigamma(xi +
2) * gamma(xi + 2) + xi * ((xi + 1)/xi^2 - 1/xi - psigamma(xi + 1, 1)) * gamma(xi + 2) + 2 * (xi + 1) * gamma(2 *
xi + 1) - ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2))/(sigma * xi^2) + 2 * ((xi + 1)^2 * gamma(2 * xi + 1) -
xi * ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2))/(sigma * xi^3))
})
k22.3_l <- sapply(xit, function(xi) {
2 * ((xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1) + (xi + 1) * gamma(2 * xi + 1) - psigamma(xi + 2) * gamma(xi +
2))/(sigma^2 * xi^2) - 2 * ((xi + 1)^2 * gamma(2 * xi + 1) - 2 * gamma(xi + 2) + 1)/(sigma^2 * xi^3)
})
k23.2_l <- sapply(xit, function(xi) {
-(-((digamma(1)) + (xi + 1)^2 * gamma(2 * xi + 1)/xi - ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2) - (gamma(xi +
2) - 1)/xi + 1)/(sigma^2 * xi^2))
})
# Linear interpolation
k123 <- approx(x = c(0, xit), y = c(k123_0, k123_l), xout = xi)$y
k222 <- approx(x = c(0, xit), y = c(k222_0, k222_l), xout = xi)$y
k13.3 <- approx(x = c(0, xit), y = c(k13.3_0, k13.3_l), xout = xi)$y
k22.3 <- approx(x = c(0, xit), y = c(k22.3_0, k22.3_l), xout = xi)$y
k23.2 <- approx(x = c(0, xit), y = c(k23.2_0, k23.2_l), xout = xi)$y
}
}
# Numerical tolerance 1e-3 If function is not too numerically unstable
if (abs(xi) > 0.001) {
k23.3 <- -((2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1)/xi - ((xi + 1)/xi + psigamma(xi + 1)) * psigamma(xi +
2) * gamma(xi + 2) + ((xi + 1)/xi^2 - 1/xi - psigamma(xi + 1, 1)) * gamma(xi + 2) - (xi + 1)^2 * gamma(2 * xi + 1)/xi^2 +
2 * (xi + 1) * gamma(2 * xi + 1)/xi - psigamma(xi + 2) * gamma(xi + 2)/xi + (gamma(xi + 2) - 1)/xi^2)/(sigma * xi^2) -
2 * ((digamma(1)) + (xi + 1)^2 * gamma(2 * xi + 1)/xi - ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2) - (gamma(xi +
2) - 1)/xi + 1)/(sigma * xi^3))
k133 <- ((1 + xi)^2 * (1 + 6 * xi + 8 * xi^2) * gamma(1 + 3 * xi) - gamma(3 + 2 * xi) * (1 + 5 * xi + 3 * xi^2 + xi * (1 +
2 * xi) * psigamma(2 + 2 * xi, 0)) + (1 + 2 * xi) * gamma(1 + xi) * (1 + 6 * xi + 5 * xi^2 + 2 * xi^3 + 2 * xi * (1 +
3 * xi + 2 * xi^2) * psigamma(2 + xi, 0) + xi^2 * (1 + xi) * psigamma(2 + xi, 0)^2 + (xi^2) * (1 + xi) * psigamma(2 +
xi, 1)))/(sigma * xi^4 * (1 + 2 * xi))
k223 <- -(1 + 2 * xi + digamma(1) * xi - xi^2 - digamma(1) * xi^2 - 27 * xi^3 * gamma(3 * xi) - 12 * xi^4 * gamma(3 * xi) -
3 * gamma(2 + xi) - 5 * xi * gamma(2 + xi) + 3 * (1 + xi) * gamma(1 + 2 * xi) + 10 * xi * (1 + xi) * gamma(1 + 2 * xi) +
4 * xi^2 * (1 + xi) * gamma(1 + 2 * xi) - gamma(1 + 3 * xi) - 6 * xi * gamma(1 + 3 * xi) - 2 * xi * gamma(2 + xi) * psigamma(2 +
xi, 0) + xi * (1 + xi) * (1 + 2 * xi) * gamma(1 + 2 * xi) * psigamma(2 + 2 * xi, 0))/(sigma^2 * xi^4)
k233 <- (1 + 7 * xi + 2 * digamma(1) * xi + 4 * xi^2 + 6 * digamma(1) * xi^2 + digamma(1)^2 * xi^2 + (pi^2 * xi^2)/6 - 3 *
xi^3 * (9 + 4 * xi) * gamma(3 * xi) + 3 * gamma(1 + 2 * xi) + 17 * xi * gamma(1 + 2 * xi) + 22 * xi^2 * gamma(1 + 2 *
xi) + 8 * xi^3 * gamma(1 + 2 * xi) - (1 + 6 * xi) * gamma(1 + 3 * xi) + (1 + 3 * xi + 2 * xi^2) * 2 * xi * gamma(1 + 2 *
xi) * psigamma(2 + 2 * xi, 0) - gamma(1 + xi) * (3 + 16 * xi + 13 * xi^2 + 4 * xi^3 + 2 * xi * (2 + 5 * xi + 3 * xi^2) *
psigamma(2 + xi, 0) + xi^2 * (1 + xi) * psigamma(2 + xi, 0)^2 + xi^2 * (1 + xi) * psigamma(2 + xi, 1)))/(sigma * xi^5)
} else {
k23.3_0 <- -3.70965809351906/sigma
k133_0 <- 0.10683192718888/sigma
k223_0 <- 1/40 * (20 * euler_gamma^4 + 3 * pi^4 - 200 * euler_gamma^3 + 20 * euler_gamma^2 * (pi^2 + 18) + 60 * pi^2 - 20 *
euler_gamma * (5 * pi^2 - 8 * zeta3 + 8) - 400 * zeta3)/sigma^2
k233_0 <- 1/48 * (12 * euler_gamma^5 - 140 * euler_gamma^4 - 21 * pi^4 + 20 * euler_gamma^3 * (pi^2 + 16) - 4 * euler_gamma^2 *
(35 * pi^2 - 60 * zeta3 + 48) + 8 * pi^2 * (5 * zeta3 - 4) + euler_gamma * (9 * pi^4 + 160 * pi^2 - 1120 * zeta3) + 288 *
zeta5 + 640 * zeta3)/sigma
if (isTRUE(all.equal(xi, 0))) {
k23.3 <- k23.3_0
k133 <- k133_0
k223 <- k223_0
k233 <- k233_0
} else {
xit <- sign(xi) * 0.01
k23.3_l <- sapply(xit, function(xi) {
-((2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1)/xi - ((xi + 1)/xi + psigamma(xi + 1)) * psigamma(xi +
2) * gamma(xi + 2) + ((xi + 1)/xi^2 - 1/xi - psigamma(xi + 1, 1)) * gamma(xi + 2) - (xi + 1)^2 * gamma(2 * xi +
1)/xi^2 + 2 * (xi + 1) * gamma(2 * xi + 1)/xi - psigamma(xi + 2) * gamma(xi + 2)/xi + (gamma(xi + 2) - 1)/xi^2)/(sigma *
xi^2) - 2 * ((digamma(1)) + (xi + 1)^2 * gamma(2 * xi + 1)/xi - ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2) -
(gamma(xi + 2) - 1)/xi + 1)/(sigma * xi^3))
})
k133_l <- sapply(xit, function(xi) {
((1 + xi)^2 * (1 + 6 * xi + 8 * xi^2) * gamma(1 + 3 * xi) - gamma(3 + 2 * xi) * (1 + 5 * xi + 3 * xi^2 + xi * (1 +
2 * xi) * psigamma(2 + 2 * xi, 0)) + (1 + 2 * xi) * gamma(1 + xi) * (1 + 6 * xi + 5 * xi^2 + 2 * xi^3 + 2 * xi *
(1 + 3 * xi + 2 * xi^2) * psigamma(2 + xi, 0) + xi^2 * (1 + xi) * psigamma(2 + xi, 0)^2 + (xi^2) * (1 + xi) * psigamma(2 +
xi, 1)))/(sigma * xi^4 * (1 + 2 * xi))
})
k223_l <- sapply(xit, function(xi) {
-(1 + 2 * xi + digamma(1) * xi - xi^2 - digamma(1) * xi^2 - 27 * xi^3 * gamma(3 * xi) - 12 * xi^4 * gamma(3 * xi) -
3 * gamma(2 + xi) - 5 * xi * gamma(2 + xi) + 3 * (1 + xi) * gamma(1 + 2 * xi) + 10 * xi * (1 + xi) * gamma(1 + 2 *
xi) + 4 * xi^2 * (1 + xi) * gamma(1 + 2 * xi) - gamma(1 + 3 * xi) - 6 * xi * gamma(1 + 3 * xi) - 2 * xi * gamma(2 +
xi) * psigamma(2 + xi, 0) + xi * (1 + xi) * (1 + 2 * xi) * gamma(1 + 2 * xi) * psigamma(2 + 2 * xi, 0))/(sigma^2 *
xi^4)
})
k233_l <- sapply(xit, function(xi) {
(1 + 7 * xi + 2 * digamma(1) * xi + 4 * xi^2 + 6 * digamma(1) * xi^2 + digamma(1)^2 * xi^2 + (pi^2 * xi^2)/6 - 3 *
xi^3 * (9 + 4 * xi) * gamma(3 * xi) + 3 * gamma(1 + 2 * xi) + 17 * xi * gamma(1 + 2 * xi) + 22 * xi^2 * gamma(1 +
2 * xi) + 8 * xi^3 * gamma(1 + 2 * xi) - (1 + 6 * xi) * gamma(1 + 3 * xi) + (1 + 3 * xi + 2 * xi^2) * 2 * xi * gamma(1 +
2 * xi) * psigamma(2 + 2 * xi, 0) - gamma(1 + xi) * (3 + 16 * xi + 13 * xi^2 + 4 * xi^3 + 2 * xi * (2 + 5 * xi +
3 * xi^2) * psigamma(2 + xi, 0) + xi^2 * (1 + xi) * psigamma(2 + xi, 0)^2 + xi^2 * (1 + xi) * psigamma(2 + xi, 1)))/(sigma *
xi^5)
})
k133 <- approx(x = c(0, xit), y = c(k133_0, k133_l), xout = xi)$y
k223 <- approx(x = c(0, xit), y = c(k223_0, k223_l), xout = xi)$y
k233 <- approx(x = c(0, xit), y = c(k233_0, k233_l), xout = xi)$y
k23.3 <- approx(x = c(0, xit), y = c(k23.3_0, k23.3_l), xout = xi)$y
}
}
# Numerical tolerance 1e-2 If function is not too numerically unstable
if (abs(xi) > 0.01) {
k333 <- (-3 * gamma(3 + 2 * xi) * (1 + 6 * xi + 4 * xi^2 + xi * (1 + 2 * xi) * psigamma(2 + 2 * xi, 0)) + 6 * (1 + 2 * xi) *
gamma(1 + xi) * (1 + 8 * xi + 7 * xi^2 + 4 * xi^3 + 2 * xi * (1 + 4 * xi + 3 * xi^2) * psigamma(2 + xi, 0) + xi^2 * (1 +
xi) * psigamma(2 + xi, 0)^2 + xi^2 * (1 + xi) * psigamma(2 + xi, 1)) + (1 + 2 * xi) * (-2 - 6 * (4 + digamma(1)) * xi -
(30 + 42 * digamma(1) + 6 * digamma(1)^2 + pi^2) * xi^2 + 2 * (1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi) + xi^3 * (-8 -
18 * digamma(1)^2 - 2 * digamma(1)^3 - 3 * pi^2 - digamma(1) * (24 + pi^2) + 4 * zeta3)))/(xi^6 * (2 + 4 * xi))
k33.3 <- 2 * ((xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1)/xi^2 - ((xi + 1)/xi + psigamma(xi + 1)) * psigamma(xi +
2) * gamma(xi + 2)/xi + ((xi + 1)/xi^2 - 1/xi - psigamma(xi + 1, 1)) * gamma(xi + 2)/xi - (xi + 1)^2 * gamma(2 * xi +
1)/xi^3 + (xi + 1) * gamma(2 * xi + 1)/xi^2 + ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2)/xi^2 - (digamma(1) + 1/xi +
1)/xi^2)/xi^2 - 1/3 * (pi^2 + 6 * (digamma(1) + 1/xi + 1)^2 + 6 * (xi + 1)^2 * gamma(2 * xi + 1)/xi^2 - 12 * ((xi + 1)/xi +
psigamma(xi + 1)) * gamma(xi + 2)/xi)/xi^3
} else {
k333_0 <- -20.8076715595589
k33.3_0 <- -5.45021409786022
if (isTRUE(all.equal(xi, 0))) {
k333 <- k333_0
k33.3 <- k33.3_0
} else {
xit <- sign(xi) * 0.01
k333_l <- sapply(xit, function(xi) {
(-3 * gamma(3 + 2 * xi) * (1 + 6 * xi + 4 * xi^2 + xi * (1 + 2 * xi) * psigamma(2 + 2 * xi, 0)) + 6 * (1 + 2 * xi) *
gamma(1 + xi) * (1 + 8 * xi + 7 * xi^2 + 4 * xi^3 + 2 * xi * (1 + 4 * xi + 3 * xi^2) * psigamma(2 + xi, 0) + xi^2 *
(1 + xi) * psigamma(2 + xi, 0)^2 + xi^2 * (1 + xi) * psigamma(2 + xi, 1)) + (1 + 2 * xi) * (-2 - 6 * (4 + digamma(1)) *
xi - (30 + 42 * digamma(1) + 6 * digamma(1)^2 + pi^2) * xi^2 + 2 * (1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi) +
xi^3 * (-8 - 18 * digamma(1)^2 - 2 * digamma(1)^3 - 3 * pi^2 - digamma(1) * (24 + pi^2) + 4 * zeta3)))/(xi^6 * (2 +
4 * xi))
})
k33.3_l <- sapply(xit, function(xi) {
2 * ((xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1)/xi^2 - ((xi + 1)/xi + psigamma(xi + 1)) * psigamma(xi +
2) * gamma(xi + 2)/xi + ((xi + 1)/xi^2 - 1/xi - psigamma(xi + 1, 1)) * gamma(xi + 2)/xi - (xi + 1)^2 * gamma(2 *
xi + 1)/xi^3 + (xi + 1) * gamma(2 * xi + 1)/xi^2 + ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2)/xi^2 - (digamma(1) +
1/xi + 1)/xi^2)/xi^2 - 1/3 * (pi^2 + 6 * (digamma(1) + 1/xi + 1)^2 + 6 * (xi + 1)^2 * gamma(2 * xi + 1)/xi^2 - 12 *
((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2)/xi)/xi^3
})
k333 <- approx(x = c(0, xit), y = c(k333_0, k333_l), xout = xi)$y
k33.3 <- approx(x = c(0, xit), y = c(k33.3_0, k33.3_l), xout = xi)$y
}
}
# Derivatives of information matrix
A1 <- 0.5 * cbind(c(k111, k112, k113), c(k112, k122, k123), c(k113, k123, k133))
A2 <- -cbind(c(k11.2, k12.2, k13.2), c(k12.2, k22.2, k23.2), c(k13.2, k23.2, k33.2)) + 0.5 * cbind(c(k112, k122, k123), c(k122,
k222, k223), c(k123, k223, k233))
A3 <- -cbind(c(k11.3, k12.3, k13.3), c(k12.3, k22.3, k23.3), c(k13.3, k23.3, k33.3)) + 0.5 * cbind(c(k113, k123, k133), c(k123,
k223, k233), c(k133, k233, k333))
# Information matrix
infomat <- gev.infomat(par = c(0, sigma, xi), dat = 1, method = "exp", nobs = 1)
infoinv <- solve(infomat)
return(infoinv %*% cbind(A1, A2, A3) %*% c(infoinv)/n)
}
#' Cox-Snell first order bias expression for the generalized Pareto distribution
#'
#' Bias vector for the GP distribution based on an \code{n} sample.
#' @inheritParams gpd
#' @references Coles, S. (2001). \emph{An Introduction to Statistical Modeling of Extreme Values}, Springer, 209 p.
#'@references Cox, D. R. and E. J. Snell (1968). A general definition of residuals, \emph{Journal of the Royal Statistical Society: Series B (Methodological)}, \strong{30}, 248--275.
#' @references Cordeiro, G. M. and R. Klein (1994). Bias correction in ARMA models, \emph{Statistics and Probability Letters}, \strong{19}(3), 169--176.
#' @references Giles, D. E., Feng, H. and R. T. Godwin (2016). Bias-corrected maximum likelihood estimation of the parameters of the generalized Pareto distribution, \emph{Communications in Statistics - Theory and Methods}, \strong{45}(8), 2465--2483.
#' @export
#' @keywords internal
#' @seealso \code{\link{gpd}}, \code{\link{gpd.bcor}}
gpd.bias <- function(par, n) {
# scale, shape
if (length(par) != 2) {
stop("Invalid input for correction")
}
if (length(n) > 1) {
stop("Invalid argument for sample size")
}
if (3 * par[2] < -1) {
stop("Invalid bias correction for GPD; need shape > -1/3")
}
c(par[1] * (3 + 5 * par[2] + 4 * par[2]^2), -(1 + par[2]) * (3 + par[2]))/(n + 3 * n * par[2])
}
#' Firth's modified score equation for the generalized Pareto distribution
#'
#' @inheritParams gpd
#' @references Firth, D. (1993). Bias reduction of maximum likelihood estimates, \emph{Biometrika}, \strong{80}(1), 27--38.
#' @seealso \code{\link{gpd}}, \code{\link{gpd.bcor}}
#' @export
#' @keywords internal
gpd.Fscore <- function(par, dat, method = c("obs", "exp")) {
if (missing(method) || method != "exp") {
method <- "obs"
}
gpd.score(par, dat) - gpd.infomat(par, dat, method) %*% gpd.bias(par, length(dat))
}
#' Firth's modified score equation for the generalized extreme value distribution
#'
#' @seealso \code{\link{gev}}
#' @inheritParams gev
#' @param method string indicating whether to use the expected ('exp') or the observed ('obs' - the default) information matrix.
#' @references Firth, D. (1993). Bias reduction of maximum likelihood estimates, \emph{Biometrika}, \strong{80}(1), 27--38.
#' @export
#' @keywords internal
gev.Fscore <- function(par, dat, method = "obs") {
if (missing(method) || method != "exp") {
method <- "obs"
}
if (par[3] < 0) {
mdat <- max(dat)
} else if (par[3] >= 0) {
mdat <- min(dat)
}
if (((mdat - par[1]) * par[3] + par[2]) < 0) {
# warning('Error in \"gev.Fscore\": data outside of range specified by parameter, yielding a zero likelihood')
return(rep(1e+08, 3))
} else {
gev.score(par, dat) - gev.infomat(par, dat, method) %*% gev.bias(par, length(dat))
}
}
#' Bias correction for GP distribution
#'
#' Bias corrected estimates for the generalized Pareto distribution using
#' Firth's modified score function or implicit bias subtraction.
#'
#'
#' Method \code{subtract} solves
#' \deqn{\tilde{\boldsymbol{\theta}} = \hat{\boldsymbol{\theta}} + b(\tilde{\boldsymbol{\theta}}}
#' for \eqn{\tilde{\boldsymbol{\theta}}}, using the first order term in the bias expansion as given by \code{\link{gpd.bias}}.
#'
#' The alternative is to use Firth's modified score and find the root of
#' \deqn{U(\tilde{\boldsymbol{\theta}})-i(\tilde{\boldsymbol{\theta}})b(\tilde{\boldsymbol{\theta}}),}
#' where \eqn{U} is the score vector, \eqn{b} is the first order bias and \eqn{i} is either the observed or Fisher information.
#'
#' The routine uses the MLE as starting value and proceeds
#' to find the solution using a root finding algorithm.
#' Since the bias-correction is not valid for \eqn{\xi < -1/3}, any solution that is unbounded
#' will return a vector of \code{NA} as the bias correction does not exist then.
#'
#' @importFrom alabama constrOptim.nl
#' @importFrom nleqslv nleqslv
#' @param par parameter vector (\code{scale}, \code{shape})
#' @param dat sample of observations
#' @param corr string indicating which correction to employ either \code{subtract} or \code{firth}
#' @param method string indicating whether to use the expected (\code{'exp'}) or the observed (\code{'obs'} --- the default) information matrix. Used only if \code{corr='firth'}
#' @return vector of bias-corrected parameters
#' @export
#' @examples
#' set.seed(1)
#' dat <- rgp(n=40, scale=1, shape=-0.2)
#' par <- gp.fit(dat, threshold=0, show=FALSE)$estimate
#' gpd.bcor(par,dat, 'subtract')
#' gpd.bcor(par,dat, 'firth') #observed information
#' gpd.bcor(par,dat, 'firth','exp')
gpd.bcor <- function(par, dat, corr = c("subtract", "firth"), method = c("obs", "exp")) {
corr <- match.arg(corr)
method <- match.arg(method)
# Basic bias correction - substract bias at MLE parbc=par-bias(par) Other bias correction - find bias corrected that solves
# implicit eqn parbc=par-bias(parbc)
if (length(par) != 2) {
stop("Invalid \"par\" argument.")
}
# Basic bias correction - substract bias at MLE parbc=par-bias(par) bcor1 <- function(par, dat){ par-gpd.bias(par,length(dat))}
# Other bias correction - find bias corrected that solves implicit eqn parbc=par-bias(parbc)
bcor <- function(par, dat) {
mdat <- max(dat)
st.opt <- par
st.opt[2] <- max(-0.25, st.opt[2])
# Constrained optimization on L2 norm squared to find decent starting value
bcor.st <- try(alabama::constrOptim.nl(par = st.opt, fn = function(parbc, para, dat) {
sum((parbc - para + gpd.bias(parbc, length(dat)))^2)
}, hin = function(x, ...) {
c(x[1], x[2] + 1/3, -x[2] + 1, x[2] * mdat + x[1])
}, dat = dat, para = par, control.optim = list(maxit = 1000), control.outer = list(trace = FALSE)), silent = TRUE)
# Check if convergence, use the starting value (otherwise use par)
if (isTRUE(try(bcor.st$value < 0.01))) {
st <- bcor.st$par
} else {
st <- par
}
if (st[2] < -1/3 || st[2] > 1) {
# If value is the MLE, make sure that it is a sensible starting value i.e. does not violate the constraints for the moments.
warning("Error in \"gpd.bcor\": invalid starting value for the shape (less than -1/3). Aborting")
return(rep(NA, 2))
}
# Root finding
bcor.rf <- try(nleqslv::nleqslv(x = st, fn = function(parbc, par, dat) {
parbc - par + gpd.bias(parbc, length(dat))
}, par = par, dat = dat), silent = TRUE)
if (!inherits(bcor.rf, what = "try-error")) {
if (bcor.rf$termcd == 1 || (bcor.rf$termcd == 2 && isTRUE(all.equal(bcor.rf$fvec, rep(0, 2), tolerance = 1e-06)))) {
return(bcor.rf$x)
}
}
return(rep(NA, 2))
}
bcorF <- function(par, dat, method = c("obs", "exp")) {
method <- match.arg(method)
if (par[2] < 0) {
mdat <- max(dat)
} else {
mdat <- min(dat)
}
st.opt <- par
st.opt[2] <- max(-0.25, st.opt[2])
bcor.st <- try(alabama::constrOptim.nl(par = st.opt, fn = function(par, dat, met) {
sum(gpd.Fscore(par = par, dat = dat, method = met)^2)
}, hin = function(x, ...) {
c(x[1], x[2] + 1/3, -x[2] + 1, x[2] * mdat + x[1])
}, dat = dat, met = method, control.optim = list(maxit = 1000), control.outer = list(trace = FALSE)), silent = TRUE)
if (isTRUE(try(bcor.st$value < 0.001))) {
st <- bcor.st$par
} else {
st <- st.opt
}
# Starting values for score-based methods, depending on feasibility par.firth <-
# try(suppressWarnings(rootSolve::multiroot(gpd.Fscore, start = st, dat = dat, method = method, positive = FALSE, maxiter =
# 1000)), silent = TRUE)
gpd.Fscoren <- function(par, met, dat) {
gpd.Fscore(par = par, method = met, dat = dat)
}
par.firth <- try(suppressWarnings(nleqslv::nleqslv(fn = gpd.Fscoren, x = st, dat = dat, met = method, control = list(maxit = 1000))),
silent = TRUE)
if (!inherits(par.firth, what = "try-error")) {
if (par.firth$termcd == 1 || (par.firth$termcd == 2 && isTRUE(all.equal(par.firth$fvec, rep(0, 2), tolerance = 1e-06)))) {
return(par.firth$x)
}
}
return(rep(NA, 2))
}
# Return values
if (corr == "subtract") {
return(bcor(par = par, dat = dat))
}
if (corr == "firth") {
return(bcorF(par = par, dat = dat, method = method))
}
}
#' Bias correction for GEV distribution
#'
#' Bias corrected estimates for the generalized extreme value distribution using
#' Firth's modified score function or implicit bias subtraction.
#'
#' Method \code{subtract}solves
#' \deqn{\tilde{\boldsymbol{\theta}} = \hat{\boldsymbol{\theta}} + b(\tilde{\boldsymbol{\theta}}}
#' for \eqn{\tilde{\boldsymbol{\theta}}}, using the first order term in the bias expansion as given by \code{\link{gev.bias}}.
#'
#' The alternative is to use Firth's modified score and find the root of
#' \deqn{U(\tilde{\boldsymbol{\theta}})-i(\tilde{\boldsymbol{\theta}})b(\tilde{\boldsymbol{\theta}}),}
#' where \eqn{U} is the score vector, \eqn{b} is the first order bias and \eqn{i} is either the observed or Fisher information.
#'
#' The routine uses the MLE (bias-corrected) as starting values and proceeds
#' to find the solution using a root finding algorithm.
#' Since the bias-correction is not valid for \eqn{\xi < -1/3}, any solution that is unbounded
#' will return a vector of \code{NA} as the solution does not exist then.
#'
#' @param par parameter vector (\code{scale}, \code{shape})
#' @param dat sample of observations
#' @param corr string indicating which correction to employ either \code{subtract} or \code{firth}
#' @param method string indicating whether to use the expected (\code{'exp'}) or the observed (\code{'obs'} --- the default) information matrix. Used only if \code{corr='firth'}
#' @return vector of bias-corrected parameters
#' @export
#' @examples
#' set.seed(1)
#' dat <- mev::rgev(n=40, loc = 1, scale=1, shape=-0.2)
#' par <- mev::fit.gev(dat)$estimate
#' gev.bcor(par, dat, 'subtract')
#' gev.bcor(par, dat, 'firth') #observed information
#' gev.bcor(par, dat, 'firth','exp')
gev.bcor <- function(par, dat, corr = c("subtract", "firth"), method = c("obs", "exp")) {
corr <- match.arg(corr)
# Basic bias correction - substract bias at MLE parbc=par-bias(par) bcor1 <- function(par, dat){ par-gpd.bias(par,length(dat))}
# Other bias correction - find bias corrected that solves implicit eqn parbc=par-bias(parbc)
bcor <- function(par, dat) {
maxdat <- max(dat)
mindat <- min(dat)
st.opt <- par
st.opt[3] <- max(-0.25, st.opt[3])
# Constrained optimization on L2 norm squared to find decent starting value
bcor.st <- try(alabama::constrOptim.nl(par = st.opt, fn = function(parbc, para, dat) {
sum((parbc - para + gev.bias(parbc, length(dat)))^2)
}, hin = function(x, ...) {
c(x[2], x[3] + 1/3, -x[3] + 1, x[3] * (maxdat - x[1]) + x[2], x[3] * (mindat - x[1]) + x[2])
}, dat = dat, para = par, control.optim = list(maxit = 1000), control.outer = list(trace = FALSE)), silent = TRUE)
# Check if convergence, use the starting value (otherwise use par)
if (isTRUE(try(bcor.st$value < 0.01))) {
st <- bcor.st$par
} else {
st <- par
}
if (st[3] < -1/3 || st[3] > 1) {
# If value is the MLE, make sure that it is a sensible starting value i.e. does not violate the constraints for the moments.
warning("Error in \"gev.bcor\": invalid starting value for the shape (less than -1/3). Aborting")
return(rep(NA, 3))
}
# Root finding
bcor.rf <- try(nleqslv::nleqslv(x = st, fn = function(parbc, par, dat) {
parbc - par + gev.bias(parbc, length(dat))
}, par = par, dat = dat, control = list(maxit = 1000, xtol = 1e-10)), silent = TRUE)
if (!inherits(bcor.rf, what = "try-error")) {
if (bcor.rf$termcd == 1 || (bcor.rf$termcd %in% c(2, 3) && isTRUE(all.equal(bcor.rf$fvec, rep(0, 3), tolerance = 1e-06)))) {
return(bcor.rf$x)
}
bcor.rf <- try(nleqslv::nleqslv(x = st, fn = function(parbc, par, dat) {
parbc - par + gev.bias(parbc, length(dat))
}, global = "none", par = par, dat = dat, control = list(maxit = 1000, xtol = 1e-10)), silent = TRUE)
if (!inherits(bcor.rf, what = "try-error")) {
if (bcor.rf$termcd == 1 || (bcor.rf$termcd %in% c(2, 3) && isTRUE(all.equal(bcor.rf$fvec, rep(0, 3), tolerance = 1e-06)))) {
return(bcor.rf$x)
} else if (abs(bcor.rf$x[3]) < 0.025 && bcor.rf$termcd == 2 && isTRUE(all.equal(bcor.rf$fvec, rep(0, 3), tolerance = 0.01))) {
warning(paste0("Approximate solution for implicit bias correction - the shape is close to zero"))
return(bcor.rf$x)
}
}
}
return(rep(NA, 3))
}
bcorF <- function(par, dat, method = c("obs", "exp")) {
method <- match.arg(method)
maxdat <- max(dat)
mindat <- min(dat)
st.opt <- par
st.opt[3] <- max(-0.25, st.opt[3])
bcor.st <- try(alabama::constrOptim.nl(par = st.opt, fn = function(par, dat, met) {
sum(gev.Fscore(par = par, dat = dat, method = met)^2)
}, hin = function(x, ...) {
c(x[2], x[3] + 0.3, -x[3] + 1, x[3] * (maxdat - x[1]) + x[2], x[3] * (mindat - x[1]) + x[2])
}, dat = dat, met = method, control.optim = list(maxit = 1000), control.outer = list(trace = FALSE)), silent = TRUE)
if (isTRUE(try(bcor.st$value < 0.001))) {
st <- bcor.st$par
} else {
st <- st.opt
}
# Starting values for score-based methods, depending on feasibility par.firth <-
# try(suppressWarnings(rootSolve::multiroot(gev.Fscore, start = st, dat = dat, method = method, positive = FALSE, maxiter =
# 1000)), silent = TRUE)
gev.Fscoren <- function(par, met, dat) {
gev.Fscore(par = par, method = met, dat = dat)
}
par.firth <- try(suppressWarnings(nleqslv::nleqslv(fn = gev.Fscoren, x = st, dat = dat, met = method, control = list(maxit = 1000,
xtol = 1e-10))), silent = TRUE)
if (!inherits(par.firth, what = "try-error")) {
# Try finding the root with default options
if (par.firth$termcd == 1 || (par.firth$termcd %in% c(2, 3) && isTRUE(all.equal(par.firth$fvec, rep(0, 3), tolerance = 1e-06)))) {
return(par.firth$x)
}
# Sometimes, method fails for values of xi close to zero - try with a full Broyden or Newton search
par.firth <- try(suppressWarnings(nleqslv::nleqslv(fn = gev.Fscoren, x = st, dat = dat, met = method, control = list(maxit = 1000,
xtol = 1e-10), global = "none")), silent = TRUE)
if (!inherits(par.firth, what = "try-error")) {
if (par.firth$termcd == 1 || (par.firth$termcd %in% c(2, 3) && isTRUE(all.equal(par.firth$fvec, rep(0, 3), tolerance = 1e-06)))) {
return(par.firth$x)
} else if (abs(par.firth$x[3]) < 0.025 && par.firth$termcd == 2 && isTRUE(all.equal(par.firth$fvec, rep(0, 3), tolerance = 0.05))) {
warning(paste0("Approximate solution for Firth's score equation with method \"", method, "\" - the shape is close to zero"))
return(par.firth$x)
}
}
}
return(rep(NA, 3))
}
# Return values
if (corr == "subtract") {
return(bcor(par = par, dat = dat))
}
if (corr == "firth") {
return(bcorF(par = par, dat = dat, method = method))
}
}
#' Posterior predictive distribution and density for the GEV distribution
#'
#' This function calculates the posterior predictive density at points x
#' based on a matrix of posterior density parameters
#'
#' @param x \code{n} vector of points
#' @param posterior \code{n} by \code{3} matrix of posterior samples
#' @param Nyr number of years to extrapolate
#' @param type string indicating whether to return the posterior \code{density} or the \code{quantile}.
#' @return a vector of values for the posterior predictive density or quantile at \code{x}, depending on the value of \code{type}
#' @export
#' @keywords internal
.gev.postpred <- function(x, posterior, Nyr = 100, type = c("density", "quantile")) {
rowMeans(cbind(apply(rbind(posterior), 1, function(par) {
switch(type, density = mev::dgev(x = x, loc = par[1] - par[2] * (1 - Nyr^par[3])/par[3], scale = par[2] * Nyr^par[3], shape = par[3]),
quantile = mev::qgev(p = x, loc = par[1] - par[2] * (1 - Nyr^par[3])/par[3], scale = par[2] * Nyr^par[3], shape = par[3]))
})))
}
#' N-year return levels, median and mean estimate
#'
#' @param par vector of location, scale and shape parameters for the GEV distribution
#' @param nobs integer number of observation on which the fit is based
#' @param N integer number of observations for return level. See \strong{Details}
#' @param type string indicating the statistic to be calculated (can be abbreviated).
#' @param p probability indicating the return level, corresponding to the quantile at 1-1/p
#'
#' @details If there are \eqn{n_y} observations per year, the \code{L}-year return level is obtained by taking
#' \code{N} equal to \eqn{n_yL}.
#' @export
#' @return a list with components
#' \itemize{
#' \item{est} point estimate
#' \item{var} variance estimate based on delta-method
#' \item{type} statistic
#' }
gev.Nyr <- function(par, nobs, N, type = c("retlev", "median", "mean"), p = 1/N) {
# Create copy of parameters
par <- as.numeric(par)
mu <- par[1]
sigma <- par[2]
xi <- par[3]
type <- match.arg(type)
# Check whether arguments are well defined
if (type == "retlev") {
stopifnot(p >= 0, p < 1)
yp <- -log(1 - p)
} else {
stopifnot(N >= 1)
}
# Euler-Masc. constant :
emcst <- -psigamma(1)
# Return levels, N-year median and mean for GEV
estimate <- as.numeric(switch(type, retlev = ifelse(xi == 0, mu - sigma * log(yp), mu - sigma/xi * (1 - yp^(-xi))), median = ifelse(xi ==
0, mu + sigma * (log(N) - log(log(2))), mu - sigma/xi * (1 - (N/log(2))^xi)), mean = ifelse(xi == 0, mu + sigma * (log(N) +
emcst), mu - sigma/xi * (1 - N^xi * gamma(1 - xi)))))
if (type == "retlev") {
if (p > 0) {
if (xi == 0) {
grad_retlev <- c(1, -log(yp), 0.5 * sigma * log(yp)^2)
} else {
grad_retlev <- c(1, -(1 - yp^(-xi))/xi, sigma * (1 - yp^(-xi))/xi^2 - sigma/xi * yp^(-xi) * log(yp))
}
}
if (p == 0) {
if (xi < 0) {
grad_retlev <- c(1, -1/xi, sigma/xi^2)
} else {
stop("Invalid argument; maximum likelihood estimate of the uptyper endpoint is Inf")
}
}
# Variance estimate based on delta-method
var_retlev <- t(grad_retlev) %*% solve(gev.infomat(par = par, dat = 1, method = "exp", nobs = nobs)) %*% grad_retlev
} else if (type == "median") {
# Gradient of N-years maxima median
if (xi == 0) {
grad_Nmed <- c(1, log(N/log(2)), 0.5 * sigma * log(N/log(2))^2)
} else {
grad_Nmed <- c(1, ((N/log(2))^xi - 1)/xi, sigma * (N/log(2))^xi * log(N/log(2))/xi - sigma * ((N/log(2))^xi - 1)/xi^2)
}
# Delta-method covariance matrix
var_Nmed <- t(grad_Nmed) %*% solve(gev.infomat(par = par, dat = 1, method = "exp", nobs = nobs)) %*% grad_Nmed
} else if (type == "mean") {
if (xi == 0) {
grad_Nmean <- c(1, log(N) + emcst, 0.5 * sigma * (emcst^2 + pi^2/6 + 2 * emcst * log(N) + log(N)^2))
} else {
grad_Nmean <- c(1, (N^xi * gamma(-xi + 1) - 1)/xi, (N^xi * log(N) * gamma(-xi + 1) - N^xi * digamma(-xi + 1) * gamma(-xi +
1)) * sigma/xi - (N^xi * gamma(-xi + 1) - 1) * sigma/xi^2)
}
var_Nmean <- t(grad_Nmean) %*% solve(gev.infomat(par = par, dat = 1, method = "exp", nobs = nobs)) %*% grad_Nmean
}
var_est <- switch(type, retlev = var_retlev, median = var_Nmed, mean = var_Nmean)
return(list(est = estimate, var = var_est[1, 1], type = type))
}
#' Asymptotic bias of block maxima for fixed sample sizes
#'
#' @param shape shape parameter
#' @param rho second-order parameter, non-positive
#' @references Dombry, C. and A. Ferreira (2017). Maximum likelihood estimators based on the block maxima method. \code{https://arxiv.org/abs/1705.00465}
#' @export
#' @return a vector of length three containing the bias for location, scale and shape (in this order)
gev.abias <- function(shape, rho) {
stopifnot(rho <= 0, shape > -0.5)
if (shape != 0 && rho < 0) {
bmu <- (1 + shape)/(shape * rho * (shape + rho)) * (-(shape + rho) * gamma(1 + shape) + (1 + shape) * rho * gamma(1 + 2 *
shape) + shape * (1 - rho) * gamma(1 + shape - rho))
bsigma <- (-shape - rho + (1 + shape) * (shape + 2 * rho) * gamma(1 + shape) - (1 + shape)^2 * rho * gamma(1 + 2 * shape) +
shape * gamma(2 - rho) - shape * (1 + shape) * (1 - rho) * gamma(1 + shape - rho))/(shape^2 * rho * (shape + rho))
bshape <- ((shape + rho) * (1 + shape + psigamma(1) * shape) - (shape + shape^2 * (1 + rho) + 2 * rho * (1 + shape)) * gamma(1 +
shape) + (1 + shape)^2 * rho * gamma(1 + 2 * shape) + shape^2 * gamma(1 - rho) - shape * (1 + shape) * gamma(2 - rho) +
shape * (1 + shape) * (1 - rho) * gamma(1 + shape - rho) - shape * rho * psigamma(2 + shape) * gamma(2 + shape) - shape^2 *
psigamma(2 - rho) * gamma(2 - rho))/(shape^3 * rho * (shape + rho))
} else if (rho == 0 && shape != 0) {
bmu <- (1 + shape)/shape^2 * ((1 + shape) * gamma(1 + 2 * shape) - gamma(2 + shape) - shape * psigamma(1 + shape) * gamma(1 +
shape))
bsigma <- (shape^2 * gamma(1 + shape) * psigamma(1 + shape) - (shape + 1)^2 * gamma(1 + 2 * shape) + shape^2 * gamma(1 + shape) +
shape * gamma(1 + shape) * psigamma(1 + shape) - (psigamma(1) + 1) * shape + 3 * shape * gamma(1 + shape) + 2 * gamma(1 +
shape) - 1)/shape^3
bshape <- ((1 + shape + shape * psigamma(1))^2 + shape^2 * pi^2/6 + (1 + shape)^2 * gamma(1 + 2 * shape) - 2 * (1 + shape) *
((1 + shape) * gamma(1 + shape) + shape * psigamma(1 + shape) * gamma(1 + shape)))/shape^4
} else if (rho < 0 && shape == 0) {
bmu <- (-1 + rho + psigamma(1) * rho + (1 - rho) * gamma(1 - rho))/rho^2
bsigma <- (6 - 6 * rho - 6 * psigamma(1)^2 * rho - pi^2 * rho + 6 * psigamma(1) * (1 - 2 * rho) - 6 * (1 - rho) * gamma(1 -
rho) * (1 + psigamma(1 - rho)))/(6 * rho^2)
bshape <- -1/(12 * rho^3) * (-12 * gamma(2 - rho) * psigamma(2 - rho) + 6 * gamma(1 - rho) * (2 + 2 * (-1 + rho)^2 * psigamma(1 -
rho) + (-1 + rho) * rho * psigamma(1 - rho)^2 + (-1 + rho) * rho * psigamma(1 - rho, deriv = 1)) + rho * (psigamma(1)^2 *
(6 - 18 * rho) + pi^2 * (1 - 3 * rho) + 6 * (-psigamma(1)^3) * rho + 3 * (-psigamma(1)) * (-4 + (4 + pi^2) * rho) + 6 *
rho * (-2 - 2 * psigamma(1, deriv = 2) + psigamma(2, 2))))
} else if (rho == 0 && shape == 0) {
# Last row of information matrix with (0, 1, shape)
bmu <- 0.41184033042644
bsigma <- 0.332484907160274
bshape <- 2.42360605517703
}
c(solve(gev.infomat(c(0, 1, shape), dat = 1, method = "exp", nobs = 1)) %*% c(bmu, bsigma, bshape))
}
#' Asymptotic bias of threshold exceedances for k order statistics
#'
#' The formula given in de Haan and Ferreira, 2007 (Springer). Note that the latter differs from that found in Drees, Ferreira and de Haan.
#' @references Dombry, C. and A. Ferreira (2017). Maximum likelihood estimators based on the block maxima method. \code{https://arxiv.org/abs/1705.00465}
#' @param shape shape parameter
#' @param rho second-order parameter, non-positive
#' @export
#' @return a vector of length containing the bias for scale and shape (in this order)
gpd.abias <- function(shape, rho) {
stopifnot(rho <= 0, shape > -0.5)
c(-rho/((1 - rho) * (1 + shape - rho)), (1 + shape)/((1 - rho) * (1 + shape - rho)))
}
Try the mev package in your browser
Any scripts or data that you put into this service are public.
mev documentation built on May 29, 2024, 9:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions gpd.abias gev.abias gev.Nyr .gev.postpred gev.bcor gpd.bcor gev.Fscore gpd.Fscore gpd.bias gev.bias
Documented in gev.abias gev.bcor gev.bias gev.Fscore gev.Nyr .gev.postpred gpd.abias gpd.bcor gpd.bias gpd.Fscore
#' Cox-Snell first order bias for the GEV distribution
#'
#' Bias vector for the GEV distribution based on an \code{n} sample.
#' Due to numerical instability, values of the information matrix and the bias
#' are linearly interpolated when the value of the shape parameter is close to zero.
#' @inheritParams gev
#' @export
#' @keywords internal
#' @seealso \code{\link{gev}}
gev.bias <- function(par, n) {
if (length(n) > 1) {
stop("Invalid argument for sample size")
}
if (length(par) != 3) {
stop("Invalid argument for parameter vector, must be a vector of length 3.")
}
sigma <- par[2]
xi <- par[3]
if (xi < -1/3) {
stop("Cox-Snell correction only valid if the shape is greater than -1/3")
}
zeta3 <- 1.20205690315959
zeta5 <- 1.03692775514337 #gsl::zeta(5)
k111 <- ((1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi))/sigma^3
k11.2 <- -2 * (xi + 1)^2 * gamma(2 * xi + 1)/sigma^3
k11.3 <- 2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1)/sigma^2 + 2 * (xi + 1) * gamma(2 * xi + 1)/sigma^2
k33.2 <- 0
euler_gamma <- -psigamma(1)
# Numerical tolerance 1e-10 (all.equal has this)
if (abs(xi) > 1e-10) {
k112 <- (xi + 1) * (gamma(2 * xi + 2) - (xi + 1) * (4 * xi + 1) * gamma(3 * xi + 1))/(sigma^3 * xi)
k12.2 <- 2 * ((xi + 1)^2 * gamma(2 * xi + 1) - gamma(xi + 2))/(sigma^3 * xi)
} else {
k112 <- (euler_gamma - 3)/sigma^3
k12.2 <- -2 * (euler_gamma - 1)/sigma^3
}
# Numerical tolerance 1e-6 If function is not too numerically unstable
if (abs(xi) > 1e-06) {
k113 <- (1 + xi) * ((1 + xi) * (1 + 4 * xi) * gamma(1 + 3 * xi) - gamma(1 + 2 * xi) * (1 + 2 * xi * (2 + xi) + xi * (1 + 2 *
xi) * psigamma(2 + 2 * xi, 0)))/(sigma^2 * xi^2)
k122 <- ((1 - xi) * gamma(2 + xi) - gamma(3 + 2 * xi) + (1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi))/(sigma^3 * xi^2)
k12.3 <- -(2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1) + 2 * (xi + 1) * gamma(2 * xi + 1) - psigamma(xi + 2) *
gamma(xi + 2))/(sigma^2 * xi) + ((xi + 1)^2 * gamma(2 * xi + 1) - gamma(xi + 2))/(sigma^2 * xi^2)
k13.2 <- -(((xi + 1)^2 * gamma(2 * xi + 1) - xi * ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2))/(sigma^2 * xi^2))
k22.2 <- -2 * ((xi + 1)^2 * gamma(2 * xi + 1) - 2 * gamma(xi + 2) + 1)/(sigma^3 * xi^2)
} else {
# Compute the approximation for xi=0, (limit)
k113_0 <- -1/12 * (36 * euler_gamma - 6 * euler_gamma^2 - pi^2 - 24)/sigma^2
k122_0 <- -1/6 * (36 * euler_gamma - 6 * euler_gamma^2 - pi^2 - 24)/sigma^3
k12.3_0 <- (6 * euler_gamma - 3 * euler_gamma^2 - 1/2 * pi^2 - 2)/(2 * sigma^2)
k13.2_0 <- 1/12 * (12 * euler_gamma - 6 * euler_gamma^2 - pi^2)/sigma^2
k22.2_0 <- 1/3 * (12 * euler_gamma - 6 * euler_gamma^2 - pi^2 - 6)/sigma^3
if (isTRUE(all.equal(xi, 0))) {
# if xi ==0 with roughly precision 1e-8, then set value to this limit
k113 <- k113_0
k122 <- k122_0
k12.3 <- k12.3_0
k13.2 <- k13.2_0
k22.2 <- k22.2_0
} else {
# if xi !=0, but numerical breakdown,
xit <- sign(xi) * 1e-06
k113_l <- sapply(xit, function(xi) {
(1 + xi) * ((1 + xi) * (1 + 4 * xi) * gamma(1 + 3 * xi) - gamma(1 + 2 * xi) * (1 + 2 * xi * (2 + xi) + xi * (1 + 2 *
xi) * psigamma(2 + 2 * xi, 0)))/(sigma^2 * xi^2)
})
k122_l <- sapply(xit, function(xi) {
((1 - xi) * gamma(2 + xi) - gamma(3 + 2 * xi) + (1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi))/(sigma^3 * xi^2)
})
k12.3_l <- sapply(xit, function(xi) {
-(2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1) + 2 * (xi + 1) * gamma(2 * xi + 1) - psigamma(xi + 2) *
gamma(xi + 2))/(sigma^2 * xi) + ((xi + 1)^2 * gamma(2 * xi + 1) - gamma(xi + 2))/(sigma^2 * xi^2)
})
k13.2_l <- sapply(xit, function(xi) {
-(((xi + 1)^2 * gamma(2 * xi + 1) - xi * ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2))/(sigma^2 * xi^2))
})
k22.2_l <- sapply(xit, function(xi) {
-2 * ((xi + 1)^2 * gamma(2 * xi + 1) - 2 * gamma(xi + 2) + 1)/(sigma^3 * xi^2)
})
# use linear interpolation between two values 0 and 1e-6
k113 <- approx(x = c(0, xit), y = c(k113_0, k113_l), xout = xi)$y
k122 <- approx(x = c(0, xit), y = c(k122_0, k122_l), xout = xi)$y
k12.3 <- approx(x = c(0, xit), y = c(k12.3_0, k12.3_l), xout = xi)$y
k13.2 <- approx(x = c(0, xit), y = c(k13.2_0, k13.2_l), xout = xi)$y
k22.2 <- approx(x = c(0, xit), y = c(k22.2_0, k22.2_l), xout = xi)$y
}
}
# Numerical tolerance 1e-4 If function is not too numerically unstable
if (abs(xi) > 1e-04) {
k123 <- ((-gamma(2 + xi)) * (1 + 2 * xi + xi * psigamma(2 + xi, 0)) + (1 + xi) * ((-(1 + 5 * xi + 4 * xi^2)) * gamma(1 + 3 *
xi) + gamma(1 + 2 * xi) * (2 + 7 * xi + 3 * xi^2 + xi * (1 + 2 * xi) * psigamma(2 + 2 * xi, 0))))/(sigma^2 * xi^3)
k222 <- (1 - 3 * xi + 3 * (xi - 1) * gamma(2 + xi) + 1.5 * gamma(3 + 2 * xi) - (1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi))/(sigma^3 *
xi^3)
k13.3 <- -(-(2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1) - xi * ((xi + 1)/xi + psigamma(xi + 1)) * psigamma(xi +
2) * gamma(xi + 2) + xi * ((xi + 1)/xi^2 - 1/xi - psigamma(xi + 1, 1)) * gamma(xi + 2) + 2 * (xi + 1) * gamma(2 * xi +
1) - ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2))/(sigma * xi^2) + 2 * ((xi + 1)^2 * gamma(2 * xi + 1) - xi * ((xi +
1)/xi + psigamma(xi + 1)) * gamma(xi + 2))/(sigma * xi^3))
k22.3 <- 2 * ((xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1) + (xi + 1) * gamma(2 * xi + 1) - psigamma(xi + 2) * gamma(xi +
2))/(sigma^2 * xi^2) - 2 * ((xi + 1)^2 * gamma(2 * xi + 1) - 2 * gamma(xi + 2) + 1)/(sigma^2 * xi^3)
k23.2 <- -(-((digamma(1)) + (xi + 1)^2 * gamma(2 * xi + 1)/xi - ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2) - (gamma(xi +
2) - 1)/xi + 1)/(sigma^2 * xi^2))
} else {
# Compute the approximation for xi=0, (limit)
k123_0 <- 1/12 * (60 * euler_gamma + 6 * euler_gamma^3 - euler_gamma * pi^2 + 4 * pi^2 * (euler_gamma - 1) - 48 * euler_gamma^2 -
4 * pi^2 + 12 * zeta3 - 12)/sigma^2
k222_0 <- 1/4 * (48 * euler_gamma + 4 * euler_gamma^3 + 9 * euler_gamma * pi^2 - 4 * pi^2 * (2 * euler_gamma - 3) + pi^2 *
(euler_gamma - 1) - 36 * euler_gamma^2 - 17 * pi^2 + 8 * zeta3 - 16)/sigma^3
k13.3_0 <- (-6 * euler_gamma - 4 * euler_gamma^3 - 7/2 * euler_gamma * pi^2 + 3/2 * pi^2 * (euler_gamma - 1) + 12 * euler_gamma^2 +
7/2 * pi^2 - 8 * zeta3)/(6 * sigma)
k22.3_0 <- -1/6 * (12 * euler_gamma + 6 * euler_gamma^3 + 4 * euler_gamma * pi^2 - pi^2 * (euler_gamma - 1) - 18 * euler_gamma^2 -
4 * pi^2 + 12 * zeta3)/sigma^2
k23.2_0 <- -1/12 * (12 * euler_gamma + 6 * euler_gamma^3 - euler_gamma * pi^2 + 4 * pi^2 * (euler_gamma - 1) - 18 * euler_gamma^2 +
pi^2 + 12 * zeta3)/sigma^2
if (isTRUE(all.equal(xi, 0))) {
# if xi numerically zero, use the latter
k123 <- k123_0
k222 <- k222_0
k13.3 <- k13.3_0
k22.3 <- k22.3_0
k23.2 <- k23.2_0
} else {
# else if less than tol, but not zero, interpolate linearly
xit <- sign(xi) * 1e-04
k123_l <- sapply(xit, function(xi) {
((-gamma(2 + xi)) * (1 + 2 * xi + xi * psigamma(2 + xi, 0)) + (1 + xi) * ((-(1 + 5 * xi + 4 * xi^2)) * gamma(1 + 3 *
xi) + gamma(1 + 2 * xi) * (2 + 7 * xi + 3 * xi^2 + xi * (1 + 2 * xi) * psigamma(2 + 2 * xi, 0))))/(sigma^2 * xi^3)
})
k222_l <- sapply(xit, function(xi) {
(1 - 3 * xi + 3 * (xi - 1) * gamma(2 + xi) + 1.5 * gamma(3 + 2 * xi) - (1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi))/(sigma^3 *
xi^3)
})
k13.3_l <- sapply(xit, function(xi) {
-(-(2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1) - xi * ((xi + 1)/xi + psigamma(xi + 1)) * psigamma(xi +
2) * gamma(xi + 2) + xi * ((xi + 1)/xi^2 - 1/xi - psigamma(xi + 1, 1)) * gamma(xi + 2) + 2 * (xi + 1) * gamma(2 *
xi + 1) - ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2))/(sigma * xi^2) + 2 * ((xi + 1)^2 * gamma(2 * xi + 1) -
xi * ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2))/(sigma * xi^3))
})
k22.3_l <- sapply(xit, function(xi) {
2 * ((xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1) + (xi + 1) * gamma(2 * xi + 1) - psigamma(xi + 2) * gamma(xi +
2))/(sigma^2 * xi^2) - 2 * ((xi + 1)^2 * gamma(2 * xi + 1) - 2 * gamma(xi + 2) + 1)/(sigma^2 * xi^3)
})
k23.2_l <- sapply(xit, function(xi) {
-(-((digamma(1)) + (xi + 1)^2 * gamma(2 * xi + 1)/xi - ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2) - (gamma(xi +
2) - 1)/xi + 1)/(sigma^2 * xi^2))
})
# Linear interpolation
k123 <- approx(x = c(0, xit), y = c(k123_0, k123_l), xout = xi)$y
k222 <- approx(x = c(0, xit), y = c(k222_0, k222_l), xout = xi)$y
k13.3 <- approx(x = c(0, xit), y = c(k13.3_0, k13.3_l), xout = xi)$y
k22.3 <- approx(x = c(0, xit), y = c(k22.3_0, k22.3_l), xout = xi)$y
k23.2 <- approx(x = c(0, xit), y = c(k23.2_0, k23.2_l), xout = xi)$y
}
}
# Numerical tolerance 1e-3 If function is not too numerically unstable
if (abs(xi) > 0.001) {
k23.3 <- -((2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1)/xi - ((xi + 1)/xi + psigamma(xi + 1)) * psigamma(xi +
2) * gamma(xi + 2) + ((xi + 1)/xi^2 - 1/xi - psigamma(xi + 1, 1)) * gamma(xi + 2) - (xi + 1)^2 * gamma(2 * xi + 1)/xi^2 +
2 * (xi + 1) * gamma(2 * xi + 1)/xi - psigamma(xi + 2) * gamma(xi + 2)/xi + (gamma(xi + 2) - 1)/xi^2)/(sigma * xi^2) -
2 * ((digamma(1)) + (xi + 1)^2 * gamma(2 * xi + 1)/xi - ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2) - (gamma(xi +
2) - 1)/xi + 1)/(sigma * xi^3))
k133 <- ((1 + xi)^2 * (1 + 6 * xi + 8 * xi^2) * gamma(1 + 3 * xi) - gamma(3 + 2 * xi) * (1 + 5 * xi + 3 * xi^2 + xi * (1 +
2 * xi) * psigamma(2 + 2 * xi, 0)) + (1 + 2 * xi) * gamma(1 + xi) * (1 + 6 * xi + 5 * xi^2 + 2 * xi^3 + 2 * xi * (1 +
3 * xi + 2 * xi^2) * psigamma(2 + xi, 0) + xi^2 * (1 + xi) * psigamma(2 + xi, 0)^2 + (xi^2) * (1 + xi) * psigamma(2 +
xi, 1)))/(sigma * xi^4 * (1 + 2 * xi))
k223 <- -(1 + 2 * xi + digamma(1) * xi - xi^2 - digamma(1) * xi^2 - 27 * xi^3 * gamma(3 * xi) - 12 * xi^4 * gamma(3 * xi) -
3 * gamma(2 + xi) - 5 * xi * gamma(2 + xi) + 3 * (1 + xi) * gamma(1 + 2 * xi) + 10 * xi * (1 + xi) * gamma(1 + 2 * xi) +
4 * xi^2 * (1 + xi) * gamma(1 + 2 * xi) - gamma(1 + 3 * xi) - 6 * xi * gamma(1 + 3 * xi) - 2 * xi * gamma(2 + xi) * psigamma(2 +
xi, 0) + xi * (1 + xi) * (1 + 2 * xi) * gamma(1 + 2 * xi) * psigamma(2 + 2 * xi, 0))/(sigma^2 * xi^4)
k233 <- (1 + 7 * xi + 2 * digamma(1) * xi + 4 * xi^2 + 6 * digamma(1) * xi^2 + digamma(1)^2 * xi^2 + (pi^2 * xi^2)/6 - 3 *
xi^3 * (9 + 4 * xi) * gamma(3 * xi) + 3 * gamma(1 + 2 * xi) + 17 * xi * gamma(1 + 2 * xi) + 22 * xi^2 * gamma(1 + 2 *
xi) + 8 * xi^3 * gamma(1 + 2 * xi) - (1 + 6 * xi) * gamma(1 + 3 * xi) + (1 + 3 * xi + 2 * xi^2) * 2 * xi * gamma(1 + 2 *
xi) * psigamma(2 + 2 * xi, 0) - gamma(1 + xi) * (3 + 16 * xi + 13 * xi^2 + 4 * xi^3 + 2 * xi * (2 + 5 * xi + 3 * xi^2) *
psigamma(2 + xi, 0) + xi^2 * (1 + xi) * psigamma(2 + xi, 0)^2 + xi^2 * (1 + xi) * psigamma(2 + xi, 1)))/(sigma * xi^5)
} else {
k23.3_0 <- -3.70965809351906/sigma
k133_0 <- 0.10683192718888/sigma
k223_0 <- 1/40 * (20 * euler_gamma^4 + 3 * pi^4 - 200 * euler_gamma^3 + 20 * euler_gamma^2 * (pi^2 + 18) + 60 * pi^2 - 20 *
euler_gamma * (5 * pi^2 - 8 * zeta3 + 8) - 400 * zeta3)/sigma^2
k233_0 <- 1/48 * (12 * euler_gamma^5 - 140 * euler_gamma^4 - 21 * pi^4 + 20 * euler_gamma^3 * (pi^2 + 16) - 4 * euler_gamma^2 *
(35 * pi^2 - 60 * zeta3 + 48) + 8 * pi^2 * (5 * zeta3 - 4) + euler_gamma * (9 * pi^4 + 160 * pi^2 - 1120 * zeta3) + 288 *
zeta5 + 640 * zeta3)/sigma
if (isTRUE(all.equal(xi, 0))) {
k23.3 <- k23.3_0
k133 <- k133_0
k223 <- k223_0
k233 <- k233_0
} else {
xit <- sign(xi) * 0.01
k23.3_l <- sapply(xit, function(xi) {
-((2 * (xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1)/xi - ((xi + 1)/xi + psigamma(xi + 1)) * psigamma(xi +
2) * gamma(xi + 2) + ((xi + 1)/xi^2 - 1/xi - psigamma(xi + 1, 1)) * gamma(xi + 2) - (xi + 1)^2 * gamma(2 * xi +
1)/xi^2 + 2 * (xi + 1) * gamma(2 * xi + 1)/xi - psigamma(xi + 2) * gamma(xi + 2)/xi + (gamma(xi + 2) - 1)/xi^2)/(sigma *
xi^2) - 2 * ((digamma(1)) + (xi + 1)^2 * gamma(2 * xi + 1)/xi - ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2) -
(gamma(xi + 2) - 1)/xi + 1)/(sigma * xi^3))
})
k133_l <- sapply(xit, function(xi) {
((1 + xi)^2 * (1 + 6 * xi + 8 * xi^2) * gamma(1 + 3 * xi) - gamma(3 + 2 * xi) * (1 + 5 * xi + 3 * xi^2 + xi * (1 +
2 * xi) * psigamma(2 + 2 * xi, 0)) + (1 + 2 * xi) * gamma(1 + xi) * (1 + 6 * xi + 5 * xi^2 + 2 * xi^3 + 2 * xi *
(1 + 3 * xi + 2 * xi^2) * psigamma(2 + xi, 0) + xi^2 * (1 + xi) * psigamma(2 + xi, 0)^2 + (xi^2) * (1 + xi) * psigamma(2 +
xi, 1)))/(sigma * xi^4 * (1 + 2 * xi))
})
k223_l <- sapply(xit, function(xi) {
-(1 + 2 * xi + digamma(1) * xi - xi^2 - digamma(1) * xi^2 - 27 * xi^3 * gamma(3 * xi) - 12 * xi^4 * gamma(3 * xi) -
3 * gamma(2 + xi) - 5 * xi * gamma(2 + xi) + 3 * (1 + xi) * gamma(1 + 2 * xi) + 10 * xi * (1 + xi) * gamma(1 + 2 *
xi) + 4 * xi^2 * (1 + xi) * gamma(1 + 2 * xi) - gamma(1 + 3 * xi) - 6 * xi * gamma(1 + 3 * xi) - 2 * xi * gamma(2 +
xi) * psigamma(2 + xi, 0) + xi * (1 + xi) * (1 + 2 * xi) * gamma(1 + 2 * xi) * psigamma(2 + 2 * xi, 0))/(sigma^2 *
xi^4)
})
k233_l <- sapply(xit, function(xi) {
(1 + 7 * xi + 2 * digamma(1) * xi + 4 * xi^2 + 6 * digamma(1) * xi^2 + digamma(1)^2 * xi^2 + (pi^2 * xi^2)/6 - 3 *
xi^3 * (9 + 4 * xi) * gamma(3 * xi) + 3 * gamma(1 + 2 * xi) + 17 * xi * gamma(1 + 2 * xi) + 22 * xi^2 * gamma(1 +
2 * xi) + 8 * xi^3 * gamma(1 + 2 * xi) - (1 + 6 * xi) * gamma(1 + 3 * xi) + (1 + 3 * xi + 2 * xi^2) * 2 * xi * gamma(1 +
2 * xi) * psigamma(2 + 2 * xi, 0) - gamma(1 + xi) * (3 + 16 * xi + 13 * xi^2 + 4 * xi^3 + 2 * xi * (2 + 5 * xi +
3 * xi^2) * psigamma(2 + xi, 0) + xi^2 * (1 + xi) * psigamma(2 + xi, 0)^2 + xi^2 * (1 + xi) * psigamma(2 + xi, 1)))/(sigma *
xi^5)
})
k133 <- approx(x = c(0, xit), y = c(k133_0, k133_l), xout = xi)$y
k223 <- approx(x = c(0, xit), y = c(k223_0, k223_l), xout = xi)$y
k233 <- approx(x = c(0, xit), y = c(k233_0, k233_l), xout = xi)$y
k23.3 <- approx(x = c(0, xit), y = c(k23.3_0, k23.3_l), xout = xi)$y
}
}
# Numerical tolerance 1e-2 If function is not too numerically unstable
if (abs(xi) > 0.01) {
k333 <- (-3 * gamma(3 + 2 * xi) * (1 + 6 * xi + 4 * xi^2 + xi * (1 + 2 * xi) * psigamma(2 + 2 * xi, 0)) + 6 * (1 + 2 * xi) *
gamma(1 + xi) * (1 + 8 * xi + 7 * xi^2 + 4 * xi^3 + 2 * xi * (1 + 4 * xi + 3 * xi^2) * psigamma(2 + xi, 0) + xi^2 * (1 +
xi) * psigamma(2 + xi, 0)^2 + xi^2 * (1 + xi) * psigamma(2 + xi, 1)) + (1 + 2 * xi) * (-2 - 6 * (4 + digamma(1)) * xi -
(30 + 42 * digamma(1) + 6 * digamma(1)^2 + pi^2) * xi^2 + 2 * (1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi) + xi^3 * (-8 -
18 * digamma(1)^2 - 2 * digamma(1)^3 - 3 * pi^2 - digamma(1) * (24 + pi^2) + 4 * zeta3)))/(xi^6 * (2 + 4 * xi))
k33.3 <- 2 * ((xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1)/xi^2 - ((xi + 1)/xi + psigamma(xi + 1)) * psigamma(xi +
2) * gamma(xi + 2)/xi + ((xi + 1)/xi^2 - 1/xi - psigamma(xi + 1, 1)) * gamma(xi + 2)/xi - (xi + 1)^2 * gamma(2 * xi +
1)/xi^3 + (xi + 1) * gamma(2 * xi + 1)/xi^2 + ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2)/xi^2 - (digamma(1) + 1/xi +
1)/xi^2)/xi^2 - 1/3 * (pi^2 + 6 * (digamma(1) + 1/xi + 1)^2 + 6 * (xi + 1)^2 * gamma(2 * xi + 1)/xi^2 - 12 * ((xi + 1)/xi +
psigamma(xi + 1)) * gamma(xi + 2)/xi)/xi^3
} else {
k333_0 <- -20.8076715595589
k33.3_0 <- -5.45021409786022
if (isTRUE(all.equal(xi, 0))) {
k333 <- k333_0
k33.3 <- k33.3_0
} else {
xit <- sign(xi) * 0.01
k333_l <- sapply(xit, function(xi) {
(-3 * gamma(3 + 2 * xi) * (1 + 6 * xi + 4 * xi^2 + xi * (1 + 2 * xi) * psigamma(2 + 2 * xi, 0)) + 6 * (1 + 2 * xi) *
gamma(1 + xi) * (1 + 8 * xi + 7 * xi^2 + 4 * xi^3 + 2 * xi * (1 + 4 * xi + 3 * xi^2) * psigamma(2 + xi, 0) + xi^2 *
(1 + xi) * psigamma(2 + xi, 0)^2 + xi^2 * (1 + xi) * psigamma(2 + xi, 1)) + (1 + 2 * xi) * (-2 - 6 * (4 + digamma(1)) *
xi - (30 + 42 * digamma(1) + 6 * digamma(1)^2 + pi^2) * xi^2 + 2 * (1 + xi)^2 * (1 + 4 * xi) * gamma(1 + 3 * xi) +
xi^3 * (-8 - 18 * digamma(1)^2 - 2 * digamma(1)^3 - 3 * pi^2 - digamma(1) * (24 + pi^2) + 4 * zeta3)))/(xi^6 * (2 +
4 * xi))
})
k33.3_l <- sapply(xit, function(xi) {
2 * ((xi + 1)^2 * psigamma(2 * xi + 1) * gamma(2 * xi + 1)/xi^2 - ((xi + 1)/xi + psigamma(xi + 1)) * psigamma(xi +
2) * gamma(xi + 2)/xi + ((xi + 1)/xi^2 - 1/xi - psigamma(xi + 1, 1)) * gamma(xi + 2)/xi - (xi + 1)^2 * gamma(2 *
xi + 1)/xi^3 + (xi + 1) * gamma(2 * xi + 1)/xi^2 + ((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2)/xi^2 - (digamma(1) +
1/xi + 1)/xi^2)/xi^2 - 1/3 * (pi^2 + 6 * (digamma(1) + 1/xi + 1)^2 + 6 * (xi + 1)^2 * gamma(2 * xi + 1)/xi^2 - 12 *
((xi + 1)/xi + psigamma(xi + 1)) * gamma(xi + 2)/xi)/xi^3
})
k333 <- approx(x = c(0, xit), y = c(k333_0, k333_l), xout = xi)$y
k33.3 <- approx(x = c(0, xit), y = c(k33.3_0, k33.3_l), xout = xi)$y
}
}
# Derivatives of information matrix
A1 <- 0.5 * cbind(c(k111, k112, k113), c(k112, k122, k123), c(k113, k123, k133))
A2 <- -cbind(c(k11.2, k12.2, k13.2), c(k12.2, k22.2, k23.2), c(k13.2, k23.2, k33.2)) + 0.5 * cbind(c(k112, k122, k123), c(k122,
k222, k223), c(k123, k223, k233))
A3 <- -cbind(c(k11.3, k12.3, k13.3), c(k12.3, k22.3, k23.3), c(k13.3, k23.3, k33.3)) + 0.5 * cbind(c(k113, k123, k133), c(k123,
k223, k233), c(k133, k233, k333))
# Information matrix
infomat <- gev.infomat(par = c(0, sigma, xi), dat = 1, method = "exp", nobs = 1)
infoinv <- solve(infomat)
return(infoinv %*% cbind(A1, A2, A3) %*% c(infoinv)/n)
}
#' Cox-Snell first order bias expression for the generalized Pareto distribution
#'
#' Bias vector for the GP distribution based on an \code{n} sample.
#' @inheritParams gpd
#' @references Coles, S. (2001). \emph{An Introduction to Statistical Modeling of Extreme Values}, Springer, 209 p.
#'@references Cox, D. R. and E. J. Snell (1968). A general definition of residuals, \emph{Journal of the Royal Statistical Society: Series B (Methodological)}, \strong{30}, 248--275.
#' @references Cordeiro, G. M. and R. Klein (1994). Bias correction in ARMA models, \emph{Statistics and Probability Letters}, \strong{19}(3), 169--176.
#' @references Giles, D. E., Feng, H. and R. T. Godwin (2016). Bias-corrected maximum likelihood estimation of the parameters of the generalized Pareto distribution, \emph{Communications in Statistics - Theory and Methods}, \strong{45}(8), 2465--2483.
#' @export
#' @keywords internal
#' @seealso \code{\link{gpd}}, \code{\link{gpd.bcor}}
gpd.bias <- function(par, n) {
# scale, shape
if (length(par) != 2) {
stop("Invalid input for correction")
}
if (length(n) > 1) {
stop("Invalid argument for sample size")
}
if (3 * par[2] < -1) {
stop("Invalid bias correction for GPD; need shape > -1/3")
}
c(par[1] * (3 + 5 * par[2] + 4 * par[2]^2), -(1 + par[2]) * (3 + par[2]))/(n + 3 * n * par[2])
}
#' Firth's modified score equation for the generalized Pareto distribution
#'
#' @inheritParams gpd
#' @references Firth, D. (1993). Bias reduction of maximum likelihood estimates, \emph{Biometrika}, \strong{80}(1), 27--38.
#' @seealso \code{\link{gpd}}, \code{\link{gpd.bcor}}
#' @export
#' @keywords internal
gpd.Fscore <- function(par, dat, method = c("obs", "exp")) {
if (missing(method) || method != "exp") {
method <- "obs"
}
gpd.score(par, dat) - gpd.infomat(par, dat, method) %*% gpd.bias(par, length(dat))
}
#' Firth's modified score equation for the generalized extreme value distribution
#'
#' @seealso \code{\link{gev}}
#' @inheritParams gev
#' @param method string indicating whether to use the expected ('exp') or the observed ('obs' - the default) information matrix.
#' @references Firth, D. (1993). Bias reduction of maximum likelihood estimates, \emph{Biometrika}, \strong{80}(1), 27--38.
#' @export
#' @keywords internal
gev.Fscore <- function(par, dat, method = "obs") {
if (missing(method) || method != "exp") {
method <- "obs"
}
if (par[3] < 0) {
mdat <- max(dat)
} else if (par[3] >= 0) {
mdat <- min(dat)
}
if (((mdat - par[1]) * par[3] + par[2]) < 0) {
# warning('Error in \"gev.Fscore\": data outside of range specified by parameter, yielding a zero likelihood')
return(rep(1e+08, 3))
} else {
gev.score(par, dat) - gev.infomat(par, dat, method) %*% gev.bias(par, length(dat))
}
}
#' Bias correction for GP distribution
#'
#' Bias corrected estimates for the generalized Pareto distribution using
#' Firth's modified score function or implicit bias subtraction.
#'
#'
#' Method \code{subtract} solves
#' \deqn{\tilde{\boldsymbol{\theta}} = \hat{\boldsymbol{\theta}} + b(\tilde{\boldsymbol{\theta}}}
#' for \eqn{\tilde{\boldsymbol{\theta}}}, using the first order term in the bias expansion as given by \code{\link{gpd.bias}}.
#'
#' The alternative is to use Firth's modified score and find the root of
#' \deqn{U(\tilde{\boldsymbol{\theta}})-i(\tilde{\boldsymbol{\theta}})b(\tilde{\boldsymbol{\theta}}),}
#' where \eqn{U} is the score vector, \eqn{b} is the first order bias and \eqn{i} is either the observed or Fisher information.
#'
#' The routine uses the MLE as starting value and proceeds
#' to find the solution using a root finding algorithm.
#' Since the bias-correction is not valid for \eqn{\xi < -1/3}, any solution that is unbounded
#' will return a vector of \code{NA} as the bias correction does not exist then.
#'
#' @importFrom alabama constrOptim.nl
#' @importFrom nleqslv nleqslv
#' @param par parameter vector (\code{scale}, \code{shape})
#' @param dat sample of observations
#' @param corr string indicating which correction to employ either \code{subtract} or \code{firth}
#' @param method string indicating whether to use the expected (\code{'exp'}) or the observed (\code{'obs'} --- the default) information matrix. Used only if \code{corr='firth'}
#' @return vector of bias-corrected parameters
#' @export
#' @examples
#' set.seed(1)
#' dat <- rgp(n=40, scale=1, shape=-0.2)
#' par <- gp.fit(dat, threshold=0, show=FALSE)$estimate
#' gpd.bcor(par,dat, 'subtract')
#' gpd.bcor(par,dat, 'firth') #observed information
#' gpd.bcor(par,dat, 'firth','exp')
gpd.bcor <- function(par, dat, corr = c("subtract", "firth"), method = c("obs", "exp")) {
corr <- match.arg(corr)
method <- match.arg(method)
# Basic bias correction - substract bias at MLE parbc=par-bias(par) Other bias correction - find bias corrected that solves
# implicit eqn parbc=par-bias(parbc)
if (length(par) != 2) {
stop("Invalid \"par\" argument.")
}
# Basic bias correction - substract bias at MLE parbc=par-bias(par) bcor1 <- function(par, dat){ par-gpd.bias(par,length(dat))}
# Other bias correction - find bias corrected that solves implicit eqn parbc=par-bias(parbc)
bcor <- function(par, dat) {
mdat <- max(dat)
st.opt <- par
st.opt[2] <- max(-0.25, st.opt[2])
# Constrained optimization on L2 norm squared to find decent starting value
bcor.st <- try(alabama::constrOptim.nl(par = st.opt, fn = function(parbc, para, dat) {
sum((parbc - para + gpd.bias(parbc, length(dat)))^2)
}, hin = function(x, ...) {
c(x[1], x[2] + 1/3, -x[2] + 1, x[2] * mdat + x[1])
}, dat = dat, para = par, control.optim = list(maxit = 1000), control.outer = list(trace = FALSE)), silent = TRUE)
# Check if convergence, use the starting value (otherwise use par)
if (isTRUE(try(bcor.st$value < 0.01))) {
st <- bcor.st$par
} else {
st <- par
}
if (st[2] < -1/3 || st[2] > 1) {
# If value is the MLE, make sure that it is a sensible starting value i.e. does not violate the constraints for the moments.
warning("Error in \"gpd.bcor\": invalid starting value for the shape (less than -1/3). Aborting")
return(rep(NA, 2))
}
# Root finding
bcor.rf <- try(nleqslv::nleqslv(x = st, fn = function(parbc, par, dat) {
parbc - par + gpd.bias(parbc, length(dat))
}, par = par, dat = dat), silent = TRUE)
if (!inherits(bcor.rf, what = "try-error")) {
if (bcor.rf$termcd == 1 || (bcor.rf$termcd == 2 && isTRUE(all.equal(bcor.rf$fvec, rep(0, 2), tolerance = 1e-06)))) {
return(bcor.rf$x)
}
}
return(rep(NA, 2))
}
bcorF <- function(par, dat, method = c("obs", "exp")) {
method <- match.arg(method)
if (par[2] < 0) {
mdat <- max(dat)
} else {
mdat <- min(dat)
}
st.opt <- par
st.opt[2] <- max(-0.25, st.opt[2])
bcor.st <- try(alabama::constrOptim.nl(par = st.opt, fn = function(par, dat, met) {
sum(gpd.Fscore(par = par, dat = dat, method = met)^2)
}, hin = function(x, ...) {
c(x[1], x[2] + 1/3, -x[2] + 1, x[2] * mdat + x[1])
}, dat = dat, met = method, control.optim = list(maxit = 1000), control.outer = list(trace = FALSE)), silent = TRUE)
if (isTRUE(try(bcor.st$value < 0.001))) {
st <- bcor.st$par
} else {
st <- st.opt
}
# Starting values for score-based methods, depending on feasibility par.firth <-
# try(suppressWarnings(rootSolve::multiroot(gpd.Fscore, start = st, dat = dat, method = method, positive = FALSE, maxiter =
# 1000)), silent = TRUE)
gpd.Fscoren <- function(par, met, dat) {
gpd.Fscore(par = par, method = met, dat = dat)
}
par.firth <- try(suppressWarnings(nleqslv::nleqslv(fn = gpd.Fscoren, x = st, dat = dat, met = method, control = list(maxit = 1000))),
silent = TRUE)
if (!inherits(par.firth, what = "try-error")) {
if (par.firth$termcd == 1 || (par.firth$termcd == 2 && isTRUE(all.equal(par.firth$fvec, rep(0, 2), tolerance = 1e-06)))) {
return(par.firth$x)
}
}
return(rep(NA, 2))
}
# Return values
if (corr == "subtract") {
return(bcor(par = par, dat = dat))
}
if (corr == "firth") {
return(bcorF(par = par, dat = dat, method = method))
}
}
#' Bias correction for GEV distribution
#'
#' Bias corrected estimates for the generalized extreme value distribution using
#' Firth's modified score function or implicit bias subtraction.
#'
#' Method \code{subtract}solves
#' \deqn{\tilde{\boldsymbol{\theta}} = \hat{\boldsymbol{\theta}} + b(\tilde{\boldsymbol{\theta}}}
#' for \eqn{\tilde{\boldsymbol{\theta}}}, using the first order term in the bias expansion as given by \code{\link{gev.bias}}.
#'
#' The alternative is to use Firth's modified score and find the root of
#' \deqn{U(\tilde{\boldsymbol{\theta}})-i(\tilde{\boldsymbol{\theta}})b(\tilde{\boldsymbol{\theta}}),}
#' where \eqn{U} is the score vector, \eqn{b} is the first order bias and \eqn{i} is either the observed or Fisher information.
#'
#' The routine uses the MLE (bias-corrected) as starting values and proceeds
#' to find the solution using a root finding algorithm.
#' Since the bias-correction is not valid for \eqn{\xi < -1/3}, any solution that is unbounded
#' will return a vector of \code{NA} as the solution does not exist then.
#'
#' @param par parameter vector (\code{scale}, \code{shape})
#' @param dat sample of observations
#' @param corr string indicating which correction to employ either \code{subtract} or \code{firth}
#' @param method string indicating whether to use the expected (\code{'exp'}) or the observed (\code{'obs'} --- the default) information matrix. Used only if \code{corr='firth'}
#' @return vector of bias-corrected parameters
#' @export
#' @examples
#' set.seed(1)
#' dat <- mev::rgev(n=40, loc = 1, scale=1, shape=-0.2)
#' par <- mev::fit.gev(dat)$estimate
#' gev.bcor(par, dat, 'subtract')
#' gev.bcor(par, dat, 'firth') #observed information
#' gev.bcor(par, dat, 'firth','exp')
gev.bcor <- function(par, dat, corr = c("subtract", "firth"), method = c("obs", "exp")) {
corr <- match.arg(corr)
# Basic bias correction - substract bias at MLE parbc=par-bias(par) bcor1 <- function(par, dat){ par-gpd.bias(par,length(dat))}
# Other bias correction - find bias corrected that solves implicit eqn parbc=par-bias(parbc)
bcor <- function(par, dat) {
maxdat <- max(dat)
mindat <- min(dat)
st.opt <- par
st.opt[3] <- max(-0.25, st.opt[3])
# Constrained optimization on L2 norm squared to find decent starting value
bcor.st <- try(alabama::constrOptim.nl(par = st.opt, fn = function(parbc, para, dat) {
sum((parbc - para + gev.bias(parbc, length(dat)))^2)
}, hin = function(x, ...) {
c(x[2], x[3] + 1/3, -x[3] + 1, x[3] * (maxdat - x[1]) + x[2], x[3] * (mindat - x[1]) + x[2])
}, dat = dat, para = par, control.optim = list(maxit = 1000), control.outer = list(trace = FALSE)), silent = TRUE)
# Check if convergence, use the starting value (otherwise use par)
if (isTRUE(try(bcor.st$value < 0.01))) {
st <- bcor.st$par
} else {
st <- par
}
if (st[3] < -1/3 || st[3] > 1) {
# If value is the MLE, make sure that it is a sensible starting value i.e. does not violate the constraints for the moments.
warning("Error in \"gev.bcor\": invalid starting value for the shape (less than -1/3). Aborting")
return(rep(NA, 3))
}
# Root finding
bcor.rf <- try(nleqslv::nleqslv(x = st, fn = function(parbc, par, dat) {
parbc - par + gev.bias(parbc, length(dat))
}, par = par, dat = dat, control = list(maxit = 1000, xtol = 1e-10)), silent = TRUE)
if (!inherits(bcor.rf, what = "try-error")) {
if (bcor.rf$termcd == 1 || (bcor.rf$termcd %in% c(2, 3) && isTRUE(all.equal(bcor.rf$fvec, rep(0, 3), tolerance = 1e-06)))) {
return(bcor.rf$x)
}
bcor.rf <- try(nleqslv::nleqslv(x = st, fn = function(parbc, par, dat) {
parbc - par + gev.bias(parbc, length(dat))
}, global = "none", par = par, dat = dat, control = list(maxit = 1000, xtol = 1e-10)), silent = TRUE)
if (!inherits(bcor.rf, what = "try-error")) {
if (bcor.rf$termcd == 1 || (bcor.rf$termcd %in% c(2, 3) && isTRUE(all.equal(bcor.rf$fvec, rep(0, 3), tolerance = 1e-06)))) {
return(bcor.rf$x)
} else if (abs(bcor.rf$x[3]) < 0.025 && bcor.rf$termcd == 2 && isTRUE(all.equal(bcor.rf$fvec, rep(0, 3), tolerance = 0.01))) {
warning(paste0("Approximate solution for implicit bias correction - the shape is close to zero"))
return(bcor.rf$x)
}
}
}
return(rep(NA, 3))
}
bcorF <- function(par, dat, method = c("obs", "exp")) {
method <- match.arg(method)
maxdat <- max(dat)
mindat <- min(dat)
st.opt <- par
st.opt[3] <- max(-0.25, st.opt[3])
bcor.st <- try(alabama::constrOptim.nl(par = st.opt, fn = function(par, dat, met) {
sum(gev.Fscore(par = par, dat = dat, method = met)^2)
}, hin = function(x, ...) {
c(x[2], x[3] + 0.3, -x[3] + 1, x[3] * (maxdat - x[1]) + x[2], x[3] * (mindat - x[1]) + x[2])
}, dat = dat, met = method, control.optim = list(maxit = 1000), control.outer = list(trace = FALSE)), silent = TRUE)
if (isTRUE(try(bcor.st$value < 0.001))) {
st <- bcor.st$par
} else {
st <- st.opt
}
# Starting values for score-based methods, depending on feasibility par.firth <-
# try(suppressWarnings(rootSolve::multiroot(gev.Fscore, start = st, dat = dat, method = method, positive = FALSE, maxiter =
# 1000)), silent = TRUE)
gev.Fscoren <- function(par, met, dat) {
gev.Fscore(par = par, method = met, dat = dat)
}
par.firth <- try(suppressWarnings(nleqslv::nleqslv(fn = gev.Fscoren, x = st, dat = dat, met = method, control = list(maxit = 1000,
xtol = 1e-10))), silent = TRUE)
if (!inherits(par.firth, what = "try-error")) {
# Try finding the root with default options
if (par.firth$termcd == 1 || (par.firth$termcd %in% c(2, 3) && isTRUE(all.equal(par.firth$fvec, rep(0, 3), tolerance = 1e-06)))) {
return(par.firth$x)
}
# Sometimes, method fails for values of xi close to zero - try with a full Broyden or Newton search
par.firth <- try(suppressWarnings(nleqslv::nleqslv(fn = gev.Fscoren, x = st, dat = dat, met = method, control = list(maxit = 1000,
xtol = 1e-10), global = "none")), silent = TRUE)
if (!inherits(par.firth, what = "try-error")) {
if (par.firth$termcd == 1 || (par.firth$termcd %in% c(2, 3) && isTRUE(all.equal(par.firth$fvec, rep(0, 3), tolerance = 1e-06)))) {
return(par.firth$x)
} else if (abs(par.firth$x[3]) < 0.025 && par.firth$termcd == 2 && isTRUE(all.equal(par.firth$fvec, rep(0, 3), tolerance = 0.05))) {
warning(paste0("Approximate solution for Firth's score equation with method \"", method, "\" - the shape is close to zero"))
return(par.firth$x)
}
}
}
return(rep(NA, 3))
}
# Return values
if (corr == "subtract") {
return(bcor(par = par, dat = dat))
}
if (corr == "firth") {
return(bcorF(par = par, dat = dat, method = method))
}
}
#' Posterior predictive distribution and density for the GEV distribution
#'
#' This function calculates the posterior predictive density at points x
#' based on a matrix of posterior density parameters
#'
#' @param x \code{n} vector of points
#' @param posterior \code{n} by \code{3} matrix of posterior samples
#' @param Nyr number of years to extrapolate
#' @param type string indicating whether to return the posterior \code{density} or the \code{quantile}.
#' @return a vector of values for the posterior predictive density or quantile at \code{x}, depending on the value of \code{type}
#' @export
#' @keywords internal
.gev.postpred <- function(x, posterior, Nyr = 100, type = c("density", "quantile")) {
rowMeans(cbind(apply(rbind(posterior), 1, function(par) {
switch(type, density = mev::dgev(x = x, loc = par[1] - par[2] * (1 - Nyr^par[3])/par[3], scale = par[2] * Nyr^par[3], shape = par[3]),
quantile = mev::qgev(p = x, loc = par[1] - par[2] * (1 - Nyr^par[3])/par[3], scale = par[2] * Nyr^par[3], shape = par[3]))
})))
}
#' N-year return levels, median and mean estimate
#'
#' @param par vector of location, scale and shape parameters for the GEV distribution
#' @param nobs integer number of observation on which the fit is based
#' @param N integer number of observations for return level. See \strong{Details}
#' @param type string indicating the statistic to be calculated (can be abbreviated).
#' @param p probability indicating the return level, corresponding to the quantile at 1-1/p
#'
#' @details If there are \eqn{n_y} observations per year, the \code{L}-year return level is obtained by taking
#' \code{N} equal to \eqn{n_yL}.
#' @export
#' @return a list with components
#' \itemize{
#' \item{est} point estimate
#' \item{var} variance estimate based on delta-method
#' \item{type} statistic
#' }
gev.Nyr <- function(par, nobs, N, type = c("retlev", "median", "mean"), p = 1/N) {
# Create copy of parameters
par <- as.numeric(par)
mu <- par[1]
sigma <- par[2]
xi <- par[3]
type <- match.arg(type)
# Check whether arguments are well defined
if (type == "retlev") {
stopifnot(p >= 0, p < 1)
yp <- -log(1 - p)
} else {
stopifnot(N >= 1)
}
# Euler-Masc. constant :
emcst <- -psigamma(1)
# Return levels, N-year median and mean for GEV
estimate <- as.numeric(switch(type, retlev = ifelse(xi == 0, mu - sigma * log(yp), mu - sigma/xi * (1 - yp^(-xi))), median = ifelse(xi ==
0, mu + sigma * (log(N) - log(log(2))), mu - sigma/xi * (1 - (N/log(2))^xi)), mean = ifelse(xi == 0, mu + sigma * (log(N) +
emcst), mu - sigma/xi * (1 - N^xi * gamma(1 - xi)))))
if (type == "retlev") {
if (p > 0) {
if (xi == 0) {
grad_retlev <- c(1, -log(yp), 0.5 * sigma * log(yp)^2)
} else {
grad_retlev <- c(1, -(1 - yp^(-xi))/xi, sigma * (1 - yp^(-xi))/xi^2 - sigma/xi * yp^(-xi) * log(yp))
}
}
if (p == 0) {
if (xi < 0) {
grad_retlev <- c(1, -1/xi, sigma/xi^2)
} else {
stop("Invalid argument; maximum likelihood estimate of the uptyper endpoint is Inf")
}
}
# Variance estimate based on delta-method
var_retlev <- t(grad_retlev) %*% solve(gev.infomat(par = par, dat = 1, method = "exp", nobs = nobs)) %*% grad_retlev
} else if (type == "median") {
# Gradient of N-years maxima median
if (xi == 0) {
grad_Nmed <- c(1, log(N/log(2)), 0.5 * sigma * log(N/log(2))^2)
} else {
grad_Nmed <- c(1, ((N/log(2))^xi - 1)/xi, sigma * (N/log(2))^xi * log(N/log(2))/xi - sigma * ((N/log(2))^xi - 1)/xi^2)
}
# Delta-method covariance matrix
var_Nmed <- t(grad_Nmed) %*% solve(gev.infomat(par = par, dat = 1, method = "exp", nobs = nobs)) %*% grad_Nmed
} else if (type == "mean") {
if (xi == 0) {
grad_Nmean <- c(1, log(N) + emcst, 0.5 * sigma * (emcst^2 + pi^2/6 + 2 * emcst * log(N) + log(N)^2))
} else {
grad_Nmean <- c(1, (N^xi * gamma(-xi + 1) - 1)/xi, (N^xi * log(N) * gamma(-xi + 1) - N^xi * digamma(-xi + 1) * gamma(-xi +
1)) * sigma/xi - (N^xi * gamma(-xi + 1) - 1) * sigma/xi^2)
}
var_Nmean <- t(grad_Nmean) %*% solve(gev.infomat(par = par, dat = 1, method = "exp", nobs = nobs)) %*% grad_Nmean
}
var_est <- switch(type, retlev = var_retlev, median = var_Nmed, mean = var_Nmean)
return(list(est = estimate, var = var_est[1, 1], type = type))
}
#' Asymptotic bias of block maxima for fixed sample sizes
#'
#' @param shape shape parameter
#' @param rho second-order parameter, non-positive
#' @references Dombry, C. and A. Ferreira (2017). Maximum likelihood estimators based on the block maxima method. \code{https://arxiv.org/abs/1705.00465}
#' @export
#' @return a vector of length three containing the bias for location, scale and shape (in this order)
gev.abias <- function(shape, rho) {
stopifnot(rho <= 0, shape > -0.5)
if (shape != 0 && rho < 0) {
bmu <- (1 + shape)/(shape * rho * (shape + rho)) * (-(shape + rho) * gamma(1 + shape) + (1 + shape) * rho * gamma(1 + 2 *
shape) + shape * (1 - rho) * gamma(1 + shape - rho))
bsigma <- (-shape - rho + (1 + shape) * (shape + 2 * rho) * gamma(1 + shape) - (1 + shape)^2 * rho * gamma(1 + 2 * shape) +
shape * gamma(2 - rho) - shape * (1 + shape) * (1 - rho) * gamma(1 + shape - rho))/(shape^2 * rho * (shape + rho))
bshape <- ((shape + rho) * (1 + shape + psigamma(1) * shape) - (shape + shape^2 * (1 + rho) + 2 * rho * (1 + shape)) * gamma(1 +
shape) + (1 + shape)^2 * rho * gamma(1 + 2 * shape) + shape^2 * gamma(1 - rho) - shape * (1 + shape) * gamma(2 - rho) +
shape * (1 + shape) * (1 - rho) * gamma(1 + shape - rho) - shape * rho * psigamma(2 + shape) * gamma(2 + shape) - shape^2 *
psigamma(2 - rho) * gamma(2 - rho))/(shape^3 * rho * (shape + rho))
} else if (rho == 0 && shape != 0) {
bmu <- (1 + shape)/shape^2 * ((1 + shape) * gamma(1 + 2 * shape) - gamma(2 + shape) - shape * psigamma(1 + shape) * gamma(1 +
shape))
bsigma <- (shape^2 * gamma(1 + shape) * psigamma(1 + shape) - (shape + 1)^2 * gamma(1 + 2 * shape) + shape^2 * gamma(1 + shape) +
shape * gamma(1 + shape) * psigamma(1 + shape) - (psigamma(1) + 1) * shape + 3 * shape * gamma(1 + shape) + 2 * gamma(1 +
shape) - 1)/shape^3
bshape <- ((1 + shape + shape * psigamma(1))^2 + shape^2 * pi^2/6 + (1 + shape)^2 * gamma(1 + 2 * shape) - 2 * (1 + shape) *
((1 + shape) * gamma(1 + shape) + shape * psigamma(1 + shape) * gamma(1 + shape)))/shape^4
} else if (rho < 0 && shape == 0) {
bmu <- (-1 + rho + psigamma(1) * rho + (1 - rho) * gamma(1 - rho))/rho^2
bsigma <- (6 - 6 * rho - 6 * psigamma(1)^2 * rho - pi^2 * rho + 6 * psigamma(1) * (1 - 2 * rho) - 6 * (1 - rho) * gamma(1 -
rho) * (1 + psigamma(1 - rho)))/(6 * rho^2)
bshape <- -1/(12 * rho^3) * (-12 * gamma(2 - rho) * psigamma(2 - rho) + 6 * gamma(1 - rho) * (2 + 2 * (-1 + rho)^2 * psigamma(1 -
rho) + (-1 + rho) * rho * psigamma(1 - rho)^2 + (-1 + rho) * rho * psigamma(1 - rho, deriv = 1)) + rho * (psigamma(1)^2 *
(6 - 18 * rho) + pi^2 * (1 - 3 * rho) + 6 * (-psigamma(1)^3) * rho + 3 * (-psigamma(1)) * (-4 + (4 + pi^2) * rho) + 6 *
rho * (-2 - 2 * psigamma(1, deriv = 2) + psigamma(2, 2))))
} else if (rho == 0 && shape == 0) {
# Last row of information matrix with (0, 1, shape)
bmu <- 0.41184033042644
bsigma <- 0.332484907160274
bshape <- 2.42360605517703
}
c(solve(gev.infomat(c(0, 1, shape), dat = 1, method = "exp", nobs = 1)) %*% c(bmu, bsigma, bshape))
}
#' Asymptotic bias of threshold exceedances for k order statistics
#'
#' The formula given in de Haan and Ferreira, 2007 (Springer). Note that the latter differs from that found in Drees, Ferreira and de Haan.
#' @references Dombry, C. and A. Ferreira (2017). Maximum likelihood estimators based on the block maxima method. \code{https://arxiv.org/abs/1705.00465}
#' @param shape shape parameter
#' @param rho second-order parameter, non-positive
#' @export
#' @return a vector of length containing the bias for scale and shape (in this order)
gpd.abias <- function(shape, rho) {
stopifnot(rho <= 0, shape > -0.5)
c(-rho/((1 - rho) * (1 + shape - rho)), (1 + shape)/((1 - rho) * (1 + shape - rho)))
}
Try the mev package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.