# R Confidence Intervals for quantiles from Generalized Lambda Distribution

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???  ## 1 Answer 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


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)$$.

• 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. Apr 9 '21 at 11:07
• Thank you for your complete explanation. Now I could fully understand how to calculate CIs. Apr 10 '21 at 16:48