R/BOLL-G.R

In ollg: Computes some Measures of OLL-G Family of Distributions

#' The beta Odd log-logistic family of distributions (BOLL-G)
#'
#' Computes the pdf, cdf, hdf, quantile and random numbers
#' of the beta extended distribution due to Cordeiro et al. (2016) specified by the pdf
#' \deqn{f=\frac{\alpha\,g\,G^{a\,\alpha-1}\bar{G}^{b\,\alpha-1}}{B(a,b)[G^\alpha+\bar{G}^\alpha]^{a+b}}}
#' for \eqn{G} any valid continuous cdf , \eqn{\bar{G}=1-G}, \eqn{g} the corresponding pdf, \eqn{B(a, b)}, the beta function, \eqn{a, b > 0}, the shape parameter, \eqn{\alpha > 0}, the first shape parameter.
#'
#' @name BOLLG
#' @param x scaler or vector of values at which the pdf or cdf needs to be computed.
#' @param q scaler or vector of probabilities at which the quantile needs to be computed.
#' @param n number of random numbers to be generated.
#' @param alpha the value of the first shape parameter, must be positive, the default is 1.
#' @param a the value of the shape parameter, must be positive, the default is 1.
#' @param b the value of the shape parameter, must be positive, the default is 1.
#' @param G A baseline continuous cdf.
#' @param ... The baseline cdf parameters.
#' @return  \code{pbollg} gives the distribution function,
#'  \code{dbollg} gives the density,
#'  \code{qbollg} gives the quantile function,
#'  \code{hbollg} gives the hazard function and
#'  \code{rbollg} generates random variables from the The beta Odd log-logistic family of
#'  distributions (BOLL-G) for baseline cdf G.
#' @references Cordeiro, G. M., Alizadeh, M., Tahir, M. H., Mansoor, M., Bourguignon, M., Hamedani, G. G. (2016). The beta odd log-logistic generalized family of distributions. Hacettepe Journal of Mathematics and Statistics, 45(4), 1175-1202.
#' @importFrom stats numericDeriv  pnorm  rbeta pbeta  uniroot  integrate
#' @examples
#' x <- seq(0, 1, length.out = 21)
#' pbollg(x)
#' pbollg(x, alpha = 2, a = 2, b = 2, G = pbeta, shape1 = 1, shape2 = 2)
#' @export
pbollg <- function(x, alpha = 1, a = 1, b = 1, G = pnorm, ...) {
  G <- sapply(x, G, ...)
  u <- G^alpha / (G^alpha + (1 - G)^alpha)
  F0 <- pbeta(u, shape1 = a, shape2 = b) - pbeta(0, shape1 = a, shape2 = b)
  return(F0)
}


#'
#' @name BOLLG
#' @examples
#' dbollg(x, alpha = 2, a = 2, b = 2, G = pbeta, shape1 = 1, shape2 = 2)
#' curve(dbollg, -3, 3)
#' @export
dbollg <- function(x, alpha = 1, a = 1, b = 1, G = pnorm, ...) {
  G0 <- function(y) G(y, ...)
  myenv <- new.env()
  myenv$par <- list(...)
  myenv$x <- as.numeric(x)
  g0 <- numericDeriv(quote(G0(x)), "x", myenv)
  g <- diag(attr(g0, "gradient"))
  G <- sapply(x, G0)
  df <- alpha * g * G^(a * alpha - 1) * (1 - G)^(b * alpha - 1) / (beta(a, b) * ((G^alpha) + (1 - G)^alpha)^(a + b))
  return(df)
}


#'
#' @name BOLLG
#' @examples
#' qbollg(x, alpha = 2, a = 2, b = 2, G = pbeta, shape1 = 1, shape2 = 2)
#' @export
qbollg <- function(q, alpha = 1, a = 1, b = 1, G = pnorm, ...) {
  q0 <- function(x0) {
    if (x0 < 0 || x0 > 1) stop(message = "[Warning] 0 < x < 1.")
    F0 <- function(t) x0 - pbollg(t, alpha, a, b, G, ...)
    F0 <- Vectorize(F0)
    x0 <- uniroot(F0, interval = c(-1e+15, 1e+15))$root
    return(x0)
  }
  return(sapply(q, q0))
}


#'
#' @name BOLLG
#' @examples
#' n <- 10
#' rbollg(n, alpha = 2, a = 2, b = 2, G = pbeta, shape1 = 1, shape2 = 2)
#' @export
rbollg <- function(n, alpha = 1, a = 1, b = 1, G = pnorm, ...) {
  v <- rbeta(n, shape1 = a, shape2 = b)
  Q_G <- function(y) qbollg(y, alpha, a, b, G, ...)
  X <- Q_G(v^(1 / alpha) / (v^(1 / alpha) + (1 - v)^(1 / alpha)))
  return(X)
}


#'
#' @name BOLLG
#' @examples
#' hbollg(x, alpha = 2, a = 2, b = 2, G = pbeta, shape1 = 1, shape2 = 2)
#' curve(hbollg, -3, 3)
#' @export
hbollg <- function(x, alpha = 1, a = 1, b = 1, G = pnorm, ...) {
  G0 <- function(y) G(y, ...)
  myenv <- new.env()
  myenv$par <- list(...)
  myenv$x <- as.numeric(x)
  g0 <- numericDeriv(quote(G0(x)), "x", myenv)
  g <- diag(attr(g0, "gradient"))
  G <- sapply(x, G0)
  h <- alpha * g * G^(a * alpha - 1) * (1 - G)^(b * alpha - 1) / (beta(a, b) * ((G^alpha) + (1 - G)^alpha)^(a + b) * (1 - alpha * g * G^(a * alpha - 1) * (1 - G)^(b * alpha - 1) / (beta(a, b) * ((G^alpha) + (1 - G)^alpha)^(a + b))))
  return(h)
}