Intervalos de confianza en R para cuantiles de la distribución lambda generalizada

Question

Me gustaría calcular intervalos de confianza en R para cuantiles de la distribución lambda generalizada.

Steve Su (2009) presenta a continuación 2 formas de calcular los intervalos de confianza. Creo que puedo entender el método de abajo (1). Pero no puedo interpretar claramente el método (2). ¿Puede alguien ayudar a codificar en R?

(1) Método de aproximación normal-GLD

%7D)

# 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)

Creo que puedo calcular el intervalo de confianza con el código anterior.

(2) Enfoque GLD de máxima probabilidad analítica

%5E%7Bm%7D(1-F%7BX%7D(x))%5E%7Bn-m-1%7Df%7BX%7D(x))

En consecuencia, para hallar el intervalo de confianza de forma analítica, basta con resolver las siguientes ecuaciones:

En la fórmula anterior, donde es la función beta incompleta de Euler normalizada por la función beta completa.

No interpreto claramente el método anterior. Lo que he intentado hasta ahora es el código de abajo pero parece que está mal. Espero que alguien me enseñe a calcular en 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???

Answer 1

En este caso, voy a reproducir el ejemplo 3.1.1 de Su (2009) en el que calcula los intervalos de confianza del 95% para el cuantil 99 de los datos de la velocidad de la luz de Michelson 1879.

Básicamente se reduce a aplicar las fórmulas (4), (5) y (6) de Su (2009). En la siguiente R código, utilicé el gld para ajustar la distribución lambda generalizada (FMKL). El conjunto de datos denominado morley se puede encontrar en el datasets paquete en R . El código no está ciertamente muy optimizado pero parece funcionar.

Estas son las fórmulas (4), (5) y (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)
}

También se puede utilizar la función beta mencionada en el documento:

# 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)
}

Ahora, vamos a ajustar la distribución FMKL GLD a los datos:

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

El ajuste parece adecuado. Por último, utilice un solucionador de raíces para resolver las ecuaciones (5) y (6) para obtener los intervalos de confianza para el cuantil 99:

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

El intervalo de confianza es $(299.9937;\,300.141)$ . Esto se acerca mucho a los valores reportados en Su (2009), que son $(299.9936;\,300.1412)$ .