R/helper_functions.R

In mos: Simulation and Moment Computation for Order Statistics

#' @importFrom stats pbeta runif rbeta
#' @importFrom hypergeo2 genhypergeo
NULL

#' Check for Natural Numbers
#'
#' @keywords internal
.is.natural <- function(x) {
  is.numeric(x) && all(x > 0 & x %% 1 == 0)
}


#' Incomplete Beta Function
#'
#' @keywords internal
.ibeta <- function(x, a, b) beta(a, b) * pbeta(x, a, b)


#' Quantile Function for Pareto Distribution
#'
#' @keywords internal
.qpareto <- function(p, scale, shape) {
  if (any(p < 0) || any(p > 1)) stop("invalid arguments")
  if (any(c(scale, shape) <= 0)) stop("invalid arguments")
  return(scale * (1 - p)^(-1 / shape))
}


#' Recursive Computation of Moments from Complementary Beta Distribution
#'
#' @description
#' This internal function computes raw moments based on a recursive relation given
#' in Makouei et al. (2021). It is used by the main moment functions.
#'
#' @param a a positive numeric value.
#' @param b a positive integer value.
#' @param k a non-negative integer value.
#' @param s a positive integer value.
#'
#' @references
#' Makouei, R., Khamnei, H. J., & Salehi, M. (2021). \emph{Moments of order statistics and k-record values arising
#' from the complementary beta distribution with application.} Journal of Computational and Applied Mathematics,
#' 390, 113386.
#'
#' @keywords internal
.Ms_recursive <- Vectorize(function(a, b, k, s) {
  if (!(a > 0 && b > 0)) stop('both parameters "a" and "b" must be positive')
  if (!.is.natural(s)) stop("invalid arguments")
  if (!b %% 1 == 0) stop('"b" must be an integer value')
  if (s == 1) {
    hg_value <- genhypergeo(U = c(a + b, 1, a + k + 1), L = c(a + 1, a + b + k + 2), z = 1)
    return((beta(a + k + 1, b + 1) / (a * beta(a, b))) * hg_value)
  } else {
    j <- 0:(b - 1)
    term <- choose(b - 1, j) * (-1)^j * .Ms_recursive(a, b, a + k + j, s - 1)
    return((1 - s / beta(a, b) * sum(term)) / (k + 1))
  }
})


#' Quantile Function for Complementary Beta Distribution
#'
#' @keywords internal
.qcompbeta <- function(p, a, b) {
  if (any(p < 0 | p > 1)) stop('"p" must be in [0, 1]')
  if (a <= 0 || b <= 0) stop('"a" and "b" must be positive')
  return(pbeta(p, shape1 = a, shape2 = b))
}


#' Quantile Function for Kumaraswamy Distribution
#'
#' @keywords internal
.qkumar <- function(p, a, b) {
  if (any(p < 0 | p > 1)) stop('"p" must be in [0, 1]')
  if (a <= 0 || b <= 0) stop('"a" and "b" must be positive')
  return((1 - (1 - p)^(1 / b))^(1 / a))
}