2
$\begingroup$

I'd like to compute confidence intervals in R for quantiles from generalized lambda distribution.

Steve Su (2009) introduces below 2 ways to calculate confidence intervals. I think I could understand method of below (1). But I cannot interpret method (2) clearly. Can anybody help to code in R?

(1) Normal-GLD Approximation Method

# p  : probability point to calculate confidence intervals
# ci : confidence interval such as 0.95 or 0.99

alpha <- 1 - ci

n <- length(data)
# fmkl GLD parameters
lambda1
lambda2
lambda3
lambda4

q <- gld::qgl(
  p = p, lambda1 = lambda1, lambda2 = lambda2, lambda3 = lambda3,
  lambda4 = lambda4, param = "fkml", lambda5 = NULL
)
s <- sqrt(p * (1 - p)) / gld::dgl(
  x = q, lambda1 = lambda1, lambda2 = lambda2, lambda3 = lambda3,
  lambda4 = lambda4, param = "fkml", lambda5 = NULL,
  inverse.eps = .Machine$double.eps, max.iterations = 500
)
z <- qnorm(p = alpha / 2, mean = 0, sd = 1)

# Confidence Intervals
q + c(-1, 1) * z * s / sqrt(n)

I think I can calculate confidence interval by above code.

(2) Analytical-Maximum Likelihood GLD Approach

Consequently, to find the confidence interval analytically, all that required is to solve the following equations:

In above formula, where is Euler's incomplete beta function normalized by the complete Beta function.

I don't interpret above method clearly. What I tried so far is below code but it seems wrong. I hope someone teach me how to calculate in R.

# fmkl GLD parameters
lambda1
lambda2
lambda3
lambda4

n <- length(data)  # does n mean length of input data??
p <-               # probability point to calculate confidence intervals
m <- ceiling(n * p)
intervals <- qbeta(p = c(alpha / 2, 1 - alpha /2), shape1 = m + 1, shape2 = n - m, ncp = 0, lower.tail = TRUE, log.p = FALSE)

q <- gld::qgl(
  p = p, lambda1 = lambda1, lambda2 = lambda2, lambda3 = lambda3,
  lambda4 = lambda4, param = "fkml", lambda5 = NULL
)

# Confidence Intervals
q + intervals     # Correct???
$\endgroup$
2
$\begingroup$

Here, I'm going to reproduce example 3.1.1 in Su (2009) where he calculates 95% confidence intervals for the 99th quantile for the speed of light data from Michelson 1879.

It basically boils down to implementing the formulas (4), (5) and (6) from Su (2009). In the following R code, I used the gld package to fit the generalized lambda distribution (FMKL). The dataset called morley can be found in the datasets package in R. The code is certainly not very optimized but it seems to work.

Here are the formulas (4), (5) and (6):

# Formula (4) in Su (2009)
gx <- function(x, n, p, lambda) {
  m <- n*p
  gamma(n + 1)/(gamma(m + 1)*gamma(n - m))*(pgl(x, lambda1 = lambda[1], lambda2 = lambda[2], lambda3 = lambda[3], lambda4 = lambda[4]))^m*(1 - pgl(x, lambda1 = lambda[1], lambda2 = lambda[2], lambda3 = lambda[3], lambda4 = lambda[4]))^(n - m - 1)*dgl(x, lambda1 = lambda[1], lambda2 = lambda[2], lambda3 = lambda[3], lambda4 = lambda[4])
}
# Formula (5) in Su (2009)
int_fun_up <- function(a, g, n, p, alpha, lambda, lower_lim) {
  integrate(gx, lower = lower_lim, upper = a, lambda = lambda, n = n, p = p, subdivisions = 1e4L, rel.tol = 15e-10)$value - (1 - (alpha/2))
}
# Formula (6) in Su (2009)
int_fun_low <- function(a, g, n, p, alpha, lambda, lower_lim) {
  integrate(gx, lower = lower_lim, upper = a, lambda = lambda, n = n, p = p, subdivisions = 1e4L, rel.tol = 15e-10)$value - (alpha/2)
}

Alternatively, you could use the beta function mentioned in the paper:

# Formula (5) in Su (2009)
int_fun_up <- function(a, g, n, p, alpha, lambda, lower_lim) {
  Fx <- pgl(a, lambda1 = lambda[1], lambda2 = lambda[2], lambda3 = lambda[3], lambda4 = lambda[4])
  (pbeta(Fx, n*p + 1, n - n*p) - pbeta(0, n*p + 1, n - n*p)) - (1 - (alpha/2))
}
# Formula (6) in Su (2009)
int_fun_low <- function(a, g, n, p, alpha, lambda, lower_lim) {
  Fx <- pgl(a, lambda1 = lambda[1], lambda2 = lambda[2], lambda3 = lambda[3], lambda4 = lambda[4])
  (pbeta(Fx, n*p + 1, n - n*p) - pbeta(0, n*p + 1, n - n*p)) - (alpha/2)
}

Now, we're going to fit the FMKL GLD distribution to the data:

library(gld)
data(morley, package = "datasets")

morley$Speed2 <- (morley$Speed + 299000)/1000

fit <- fit.fkml(morley$Speed2, return.data = TRUE)

plot(fit, one.page = TRUE) # Check fit

GLD_fit

The fit seems adequate. Finally, use a root-solver to solve equations (5) and (6) to get the confidence intervals for the 99th quantile:

uniroot(int_fun_low, interval = c(299.5, 301), g = gx, lower_lim = 299.5, n = 100, p = 0.99, alpha = 0.05, lambda = fit$lambda)$root.
[1] 299.9937
uniroot(int_fun_up, interval = c(299.5, 301), g = gx, lower_lim = 299.5, n = 100, p = 0.99, alpha = 0.05, lambda = fit$lambda)$root
[1] 300.141

The confidence interval is $(299.9937;\,300.141)$. This is very close to the values reported in Su (2009) which are $(299.9936;\,300.1412)$.

$\endgroup$
2
  • $\begingroup$ Thank you very much for your answer. Yours are very direct encoding so your code should be fine. I was sticking to how to use beta function. Let me have some time to check. $\endgroup$
    – kabocha
    Apr 9 '21 at 11:07
  • 1
    $\begingroup$ Thank you for your complete explanation. Now I could fully understand how to calculate CIs. $\endgroup$
    – kabocha
    Apr 10 '21 at 16:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.