R/functions.R

In mos: Simulation and Moment Computation for Order Statistics

#' Generate Random Data from Order Statistics
#'
#' This function generates random data from the order statistics of a specified distribution.
#' The user can specify a known distribution in R or provide a custom quantile function.
#'
#' @param size number of observations.
#' @param r rank(s) of the desired order statistic(s) (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param dist a character string specifying the name of a known distribution in R (e.g. \code{'norm'}, \code{'exp'}). Default is \code{NULL}.
#' @param qf a custom quantile function, either as a name (string) or directly as a function. Default is \code{NULL}.
#' @param ... further arguments to be passed to \code{dist} or \code{qf}.
#'
#' @return
#' A numeric vector or matrix containing the generated random data from the specified order statistics.
#' If a single rank is provided (i.e., scalar \code{r}), a numeric vector of size \code{size} is returned.
#' If multiple ranks are provided (i.e., vector \code{r}), a matrix is returned with \code{size} rows and
#' \code{length(r)} columns, where each row corresponds to a simulation and each column to an order statistic.
#'
#' @details The \code{ros} function generates random data from order statistics using two approaches:
#' \enumerate{
#'  \item \bold{Using a Known Distribution:} When \code{dist} is provided, random data is generated from a known distribution in R.
#' \item \bold{Using a Custom Quantile Function:} When \code{qf} is provided, \code{ros} applies
#' the user-provided quantile function to generate random data.
#'   }
#'
#' @examples
#' # Example 1: Generate from the normal distribution
#' ros(5, 3, 15, "norm", mean = 4, sd = 2)
#'
#' # Example 2: Using a custom quantile function for the Pareto distribution
#' ros(3, 2, 10, qf = function(p, scale, shape) scale * (1 - p)^(-1 / shape), scale = 3, shape = 2)
#'
#' # Example 3: Generate multiple order statistics from the uniform distribution
#' # In this example, first through 5th order statistics are generated from a sample size of 5.
#' ros(3, 1:5, 5, dist = "unif")
#'
#' @export
ros <- function(size, r, n, dist = NULL, qf = NULL, ...) {
  if (is.null(size)) stop('argument "size" is missing, with no default')
  if (is.null(n)) stop('argument "n" is missing, with no default')
  if (length(size) > 1 || length(n) > 1) stop('"size" and "n" must be scalar')
  if (!all(.is.natural(c(size, r, n)))) stop("invalid arguments")
  if (any(r > n)) stop('"r" must be less than or equal to "n"')
  if (is.null(dist) && is.null(qf)) stop('either "dist" or "qf" must be provided')
  if (!is.null(dist) && !is.null(qf)) stop('specify only one of "dist" or "qf"')
  if (!is.null(qf)) {
    if (is.character(qf)) {
      if (!exists(qf, mode = "function")) {
        stop(sprintf('the quantile function "%s" does not exist in the current environment', qf))
      }
      qfun <- get(qf, mode = "function")
    } else if (is.function(qf)) {
      qfun <- qf
    } else {
      stop('"qf" must be either a function or a charachter string na,ing a function')
    }
  }
  if (length(r) == 1) {
    if (!is.null(dist)) {
      return(do.call(paste0("q", dist), list(rbeta(size, r, n + 1 - r), ...)))
    } else {
      return(do.call(qfun, list(rbeta(size, r, n + 1 - r), ...)))
    }
  } else {
    args <- list(...)
    mat <- replicate(size, expr = {
      if (!is.null(dist)) {
        sort(do.call(paste0("r", dist), c(list(n), args)))
      } else {
        sort(do.call(qfun, c(list(runif(n)), args)))
      }
    })
    out <- t(mat)[, r, drop = FALSE]
    colnames(out) <- paste0("r", r)
    return(out)
  }
}




#' Variance of Order Statistics
#'
#' This function computes the variance of order statistics for a given distribution.
#'
#' @param r rank(s) of the desired order statistic(s) (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param dist a character string specifying the name of a distribution. Supported values are:
#'   \itemize{
#'     \item \code{"unif"}: Uniform distribution
#'     \item \code{"exp"}: Exponential distribution
#'     \item \code{"weibull"}: Weibull distribution
#'     \item \code{"tri"}: Triangular distribution
#'     \item \code{"topple"}: Topp-Leon distribution
#'   }
#' @param ... further arguments to be passed to \code{dist}.
#'
#' @details This function computes the variance of the \eqn{r}th order statistic (\eqn{X_{r:n}})
#' for a given sample size (\code{n}) and distribution (\code{dist}). The variance is calculated using:
#' \deqn{\text{Var}(X_{r:n}) = \text{E}(X^2_{r:n}) - (\text{E}(X_{r:n}))^2}
#'
#' @return The variance of the \eqn{r}th order statistic.
#'
#' @examples
#' # Variance of the 3rd order statistic in a sample of size 10 from a uniform distribution
#' varOS(r = 3, n = 10, dist = "unif")
#'
#' @export
varOS <- Vectorize(function(r, n, dist = c("unif", "exp", "weibull", "tri", "topple"), ...) {
  if (!all(.is.natural(c(r, n)))) stop("invalid arguments")
  if (any(r > n)) stop('"r" must be less than or equal to "n"')
  dist <- match.arg(dist)
  second_moment <- do.call(paste0("mo_", dist), list(r, n, k = 2, ...))
  first_moment <- do.call(paste0("mo_", dist), list(r, n, k = 1, ...))
  return(second_moment - first_moment^2)
}, vectorize.args = "r")




#' Skewness of Order Statistics
#'
#' This function computes the skewness of the order statistics for a given distribution.
#'
#' @param r rank(s) of the desired order statistic(s) (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param dist a character string specifying the name of a distribution. Supported values are:
#'   \itemize{
#'     \item \code{"unif"}: Uniform distribution
#'     \item \code{"exp"}: Exponential distribution
#'     \item \code{"weibull"}: Weibull distribution
#'     \item \code{"tri"}: Triangular distribution
#'   }
#' @param ... further arguments to be passed to \code{dist}.
#'
#' @details
#' The skewness of the \eqn{r}th order statistic is calculated using the formula:
#' \deqn{\text{Skewness}(X_{r:n}) = \text{E}(\frac{X_{r:n}-\mu_{r:n}} {\sigma_{r:n}})^3}
#' where \eqn{\mu_{r:n}} and \eqn{\sigma_{r:n}} are the mean and standard deviation of the \eqn{r}th order statistic, respectively.
#'
#' @return The skewness of the \eqn{r}th order statistic.
#'
#' @seealso \link{varOS}, \link{kurtOS}
#'
#' @examples
#' # Skewness of the 3rd order statistic for a uniform distribution
#' skewOS(3, 10, "unif")
#'
#' @export
skewOS <- Vectorize(function(r, n, dist = c("unif", "exp", "weibull", "tri"), ...) {
  if (!all(.is.natural(c(r, n)))) stop("invalid arguments")
  if (any(r > n)) stop('"r" must be less than or equal to "n"')
  dist <- match.arg(dist)
  mu <- do.call(paste0("mo_", dist), list(r, n, k = 1, ...))
  var_r <- varOS(r, n, dist, ...)
  third_moment <- do.call(paste0("mo_", dist), list(r, n, k = 3, ...))
  return((third_moment - 3 * mu * var_r - mu^3) / var_r^(3 / 2))
}, vectorize.args = "r")




#' Kurtosis of Order Statistics
#'
#' This function computes the kurtosis of order statistics for a given distribution.
#'
#' @param r rank(s) of the desired order statistic(s) (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param dist a character string specifying the name of a distribution. Supported values are:
#'   \itemize{
#'     \item \code{"unif"}: Uniform distribution
#'     \item \code{"exp"}: Exponential distribution
#'     \item \code{"weibull"}: Weibull distribution
#'     \item \code{"tri"}: Triangular distribution
#'   }
#' @param ... further arguments to be passed to \code{dist}.
#'
#' @details
#' The kurtosis of the \eqn{r}th order statistic is calculated using the formula:
#' \deqn{\text{kurtosis}(X_{r:n}) = \text{E}(\frac{X_{r:n}-\mu_{r:n}} {\sigma_{r:n}})^4}
#' where \eqn{\mu_{r:n}} and \eqn{\sigma_{r:n}} are the mean and standard deviation of the \eqn{r}th order statistic, respectively.
#'
#' @return The kurtosis of the \eqn{r}th order statistic.
#'
#' @seealso \link{varOS}, \link{skewOS}
#'
#' @examples
#' # Compute the kurtosis of the 3rd order statistic from a sample of size 10
#' kurtOS(r = 3, n = 10, dist = "unif")
#'
#' @export
kurtOS <- Vectorize(function(r, n, dist = c("unif", "exp", "weibull", "tri"), ...) {
  if (!all(.is.natural(c(r, n)))) stop("invalid arguments")
  if (any(r > n)) stop('"r" must be less than or equal to "n"')
  dist <- match.arg(dist)
  mu <- do.call(paste0("mo_", dist), list(r, n, k = 1, ...))
  j <- 1:4
  s <- choose(4, j) * do.call(paste0("mo_", dist), list(r, n, k = j, ...)) * (-mu)^(4 - j)
  return((sum(mu^4, s)) / varOS(r, n, dist, ...)^2)
}, vectorize.args = "r")




#' Moments of Order Statistics from the Uniform Distribution
#'
#' This function computes the moments of order statistics for the uniform distribution
#' based on the relationship described by Arnold and Balakrishnan (2012).
#'
#' @param r rank(s) of the desired order statistic(s) (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param k order of the moment to compute (default is \code{1}).
#' @param a,b lower and upper limits of the distribution. Must be finite.
#'
#' @details
#' The function calculates the \eqn{k}th moment based on the formula:
#' \deqn{\text{E}[U_{r,n}^k] = \frac{B(k + r, n - r + 1)}{B(r, n - r + 1)},}
#' where \eqn{B(a, b)} is the complete beta function. When \eqn{a \neq 0} or \eqn{b \neq 1},
#' the transformation \eqn{U^* = a + (b - a)U} is used.
#'
#' @return
#' The \eqn{k}th moment of the \eqn{r}th order statistic from a uniform distribution.
#'
#' @references
#' Arnold, B. C., & Balakrishnan, N. (2012). \emph{Relations, bounds and approximations for
#' order statistics} (Vol. 53). Springer Science & Business Media.
#'
#' @examples
#' # Example 1: First moment (mean) of the 2nd order statistic from a sample of size 5
#' mo_unif(2, 5, k = 1, a = 0, b = 1)
#'
#' # Example 2: Second moment of the 3rd order statistic from a uniform distribution on [2, 5]
#' mo_unif(3, 7, k = 2, a = 2, b = 5)
#'
#' @export
mo_unif <- Vectorize(function(r, n, k = 1, a = 0, b = 1) {
  if (!all(.is.natural(c(r, n, k)))) stop("invalid arguments")
  if (any(r > n)) stop('"r" must be less than or equal to "n"')
  if (a == 0 && b == 1) {
    return(beta(k + r, n - r + 1) / beta(r, n - r + 1))
  } else {
    j <- 0:(k - 1)
    s <- choose(k, j) * a^j * (b - a)^(k - j) * mo_unif(r, n, k = k - j, a = 0, b = 1)
    return(sum(s, 1))
  }
}, vectorize.args = "r")




#' Moments Of Order Statistics from the Exponential Distribution
#'
#' This function computes the moments of order statistics from the exponential distribution.
#'
#' @param r rank(s) of the desired order statistic(s) (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param k order of the moment to compute (default is \code{1}).
#' @param mu location parameter of the exponential distribution (default is 0).
#' @param sigma scale parameter of the exponential distribution (default is 1).
#'
#' @details
#' The function calculates the \eqn{k}th moment using the following relationship:
#'
#' \deqn{
#' \text{E}[X_{r:n}^k] = \frac{n!}{(r-1)!(n-r)!} \cdot \sum_{j=0}^{r-1} (-1)^j \binom{r-1}{j}
#' \frac{\Gamma(k+1)}{(n-r+j+1)^{k+1}}.
#' }
#'
#' For non-standard exponential distributions with \eqn{\mu} and \eqn{\sigma} parameters,
#' the transformation \eqn{X^*=\mu + \sigma X} is used.
#'
#' @return
#' The \eqn{k}th moment of the \eqn{r}th order statistic from an exponential distribution.
#'
#' @references
#' Ahsanullah, M., Nevzorov, V. B., & Shakil, M. (2013). \emph{An introduction to order statistics} (Vol. 8). Paris: Atlantis Press.
#'
#' @examples
#' # First moment (mean) of the 2nd order statistic from a sample of size 5
#' mo_exp(2, 5, k = 1, mu = 0, sigma = 1)
#'
#' # Second moment of the 3rd order statistic from an exponential distribution
#' # with mu = 2 and sigma = 3
#' mo_exp(3, 7, k = 2, mu = 2, sigma = 3)
#'
#' @export
mo_exp <- Vectorize(function(r, n, k = 1, mu = 0, sigma = 1) {
  if (!all(.is.natural(c(r, n, k)))) stop("invalid arguments")
  if (any(r > n)) stop('"r" must be less than or equal to "n"')
  if (mu == 0 && sigma == 1) {
    j <- 0:(r - 1)
    s <- (-1)^j * choose(r - 1, j) * (gamma(k + 1) / (n - r + j + 1)^(k + 1))
    return(r * choose(n, r) * sum(s))
  } else {
    j <- 0:(k - 1)
    s <- choose(k, j) * mu^j * sigma^(k - j) * mo_exp(r, n, k = k - j, mu = 0, sigma = 1)
    return(sum(s, 1))
  }
}, vectorize.args = "r")




#' Moments of Order Statistics from the Weibull Distribution
#'
#' This function computes the moments of order statistics from the weibull distribution.
#'
#' @param r rank(s) of the desired order statistic(s) (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param k order of the moment to compute (default is \code{1}).
#' @param shape shape parameter of the weibull distribution.
#' @param scale scale parameter of the weibull distribution (default is \code{1}).
#'
#' @details
#' The function calculates the \eqn{k}th moment using the formula:
#' \deqn{
#' \text{E}[X_{r:n}^k] = \frac{n!}{(r-1)!(n-r)!} \Gamma\left(1 + \frac{k}{\text{shape}}\right)
#' \sum_{j=0}^{r-1} (-1)^j \binom{r-1}{j} \frac{1}{(n-r+1+j)^{1 + \frac{k}{\text{shape}}}}
#' }
#' For non-standard weibull distributions (\code{scale} not equal to 1), the relationship
#' \eqn{\text{E}[Z_{r:n}^k] = \text{scale}^k \text{E}[X_{r:n}^k]} is used.
#'
#' @return The \eqn{k}th moment of the \eqn{r}th order statistic from a weibull distribution.
#'
#' @references
#' Harter, H. L., & Balakrishnan, N. (1996). \emph{CRC handbook of tables for the use of order statistics in estimation.} CRC press.
#'
#' @examples
#' # Example 1: Standard weibull distribution (shape = 2, scale = 1)
#' mo_weibull(r = 2, n = 5, k = 1, shape = 2)
#'
#' # Example 2: Non-standard weibull distribution (shape = 2, scale = 3)
#' mo_weibull(r = 3, n = 6, k = 2, shape = 2, scale = 3)
#'
#' @export
mo_weibull <- Vectorize(function(r, n, k = 1, shape, scale = 1) {
  if (!all(.is.natural(c(r, n, k)))) stop("invalid arguments")
  if (any(r > n)) stop('"r" must be less than or equal to "n"')
  if (scale == 1) {
    j <- 0:(r - 1)
    s <- (-1)^j * choose(r - 1, j) / (n - r + 1 + j)^(1 + (k / shape))
    return(r * choose(n, r) * gamma(1 + (k / shape)) * sum(s))
  } else {
    return(scale^k * mo_weibull(r, n, k, shape, scale = 1))
  }
}, vectorize.args = "r")




#' Moments of Order Statistics from the Symmetric Triangular Distribution
#'
#' This function computes the moments of order statistics from the symmetric triangular distribution.
#'
#' @param r rank(s) of the desired order statistic(s) (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param k order of the moment to compute (default is \code{1}).
#'
#' @return The \eqn{k}th moment of the \eqn{r}th order statistic from a symmetric triangular distribution.
#'
#' @details
#' The function implements the following relationship from Nagaraja (2013)
#' for the symmetric triangular distribution:
#'
#' \deqn{
#' \text{E}[X_{r:n}^k] = \frac{n!}{(r-1)!(n-r)!} \left\{
#' (\frac{1}{2})^{k/2} B\left(\frac{1}{2}; \frac{k}{2} + r, n - r + 1\right) +
#' \sum_{j=0}^k (-1)^j \binom{k}{j} (\frac{1}{2})^{j/2}
#' B\left(\frac{1}{2}; \frac{j}{2} + n - r + 1, r\right)
#' \right\}.
#' }
#'
#' Here, \eqn{B(x; a, b)} is the incomplete Beta function.
#'
#' @references
#' Nagaraja, H. N. (2013). \emph{Moments of order statistics and L-moments for the symmetric
#' triangular distribution.} Statistics & Probability Letters, 83(10), 2357-2363.
#'
#' @examples
#' # Compute the 2nd moment of the 3rd order statistic for n=5
#' mo_tri(3, 5, 2)
#'
#' @export
mo_tri <- Vectorize(function(r, n, k = 1) {
  if (!all(.is.natural(c(r, n, k)))) stop("invalid arguments")
  if (any(r > n)) stop('"r" must be less than or equal to "n"')
  j <- 0:k
  s <- (-1)^j * choose(k, j) * 1 / 2^(j / 2) * .ibeta(1 / 2, j / 2 + n - r + 1, r)
  return(r * choose(n, r) * ((1 / 2^(k / 2) * .ibeta(1 / 2, k / 2 + r, n - r + 1)) + sum(s)))
}, vectorize.args = "r")




#' Moments of Order Statistics from the Topp-Leone Distribution
#'
#' This function computes the moments of order statistic from the topp-leone distribution, based on the formula
#' presented in Genç, A. İ. (2012).
#'
#' @param r rank(s) of the desired order statistic(s) (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param k order of the moment to compute (default is \code{1}).
#' @param a shape parameter of the topp-leone distribution (\code{a} > \code{0}).
#' @param b scale parameter of the topp-leone distribution (default is \code{1}, \code{b} > \code{0}).
#'
#' @return The \eqn{k}th moment of the \eqn{r}th order statistic from a topp-leone distribution.
#'
#' @details
#' This function implements the exact formula for moments of order statistics from the
#' topp-leone distribution as provided in Genç, A. İ. (2012):
#' \deqn{
#' \text{E}[X_{r:n}^k] = \frac{n!(ab^k)}{(r-1)!(n-r)!} \sum_{j=0}^{n-r} \binom{n-r}{j} (-1)^j
#' 2^{k + 2a(r+j)} \left[ B_{1/2}(k + a(r+j), a(r+j)) - 2 B_{1/2}(k + a(r+j) + 1, a(r+j)) \right]
#' }
#' Here, \eqn{B_x(., .)} is the incomplete Beta function.
#'
#' @references
#' Genç, A. İ. (2012). \emph{Moments of order statistics of Topp–Leone distribution.}
#' Statistical Papers, 53, 117-131.
#'
#' @examples
#' # Compute the first moment of the first order statistic for n=5, a=2, b=1
#' mo_topple(1, 5, 1, 2)
#'
#' # Compute the second moment of the second order statistic for n=10, a=1.5, b=2
#' mo_topple(2, 10, 2, 1.5, 2)
#'
#' @export
mo_topple <- Vectorize(function(r, n, k = 1, a, b = 1) {
  if (!all(.is.natural(c(r, n, k)))) stop("invalid arguments")
  if (any(r > n)) stop('"r" must be less than or equal to "n"')
  if (!(a > 0 && b > 0)) stop("invalid arguments")
  j <- 0:(n - r)
  term <- choose(n - r, j) * (-1)^j * 2^(k + 2 * a * (r + j)) * (.ibeta(1 / 2, k + a * (r + j), a * (r + j)) - 2 * .ibeta(1 / 2, k + a * (r + j) + 1, a * (r + j)))
  return((factorial(n) * (a * b^k)) / (factorial(r - 1) * factorial(n - r)) * sum(term))
}, vectorize.args = "r")




#' Moments of Order Statistics from the Complementary Beta Distribution
#'
#' This function computes the moments of order statistics from the complementary beta (CB) distribution.
#' For small values of \code{k} and integer \code{b}, a closed-form formula is used; otherwise,
#' Monte Carlo simulation is applied.
#'
#' @param r rank of the desired order statistic (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param k order of the moment to compute (default is \code{1}).
#' @param a,b positive parameters of the complementary beta distribution.
#' @param rep number of simulations (used when \code{b} is non-integer, default is \code{1e5}).
#' @param seed optional seed for random number generation to ensure reproducibility (used when \code{b} is non-integer, default is \code{42}).
#' @param verbose logical; if \code{TRUE}, prints a message when Monte Carlo simulation is used.
#'
#' @details
#' The computation method varies depending on \code{b} and \code{k}:
#' \itemize{
#'   \item \strong{For integer \code{b} and \code{k = 1, 2}}: The function calculates the moments using the closed-form expression derived in Makouei et al. (2021):
#'   \deqn{\text{E}[X_{r:n}^s] = \frac{1}{B(r, n - r + 1)} \sum_{j=0}^{n-r} \binom{n - r}{j} (-1)^j \mathcal{M}^{(s)}(a, b, r + j - 1),}
#'   Here
#'   \deqn{\mathcal{M}^{(s)}(a, b, k) = \frac{1}{k + 1} \left[1 - \frac{s}{B(a, b)} \sum_{j=0}^{\infty} \binom{b-1}{j}
#'   (-1)^j \mathcal{M}^{(s-1)}(a, b, a + k + j) \right], \quad s \geq 1,}
#'   with the starting point
#'   \deqn{\mathcal{M}^{(1)}(a, b, k) = \frac{B(a + k + 1, b + 1)}{a B(a, b)} \cdot
#'   {_3F_2}\left(a + b, 1, a + k + 1; a + 1, a + b + k + 2; 1\right),}
#'   where \eqn{B(a, b)} is the beta function, \eqn{_3F_2} is the generalized hypergeometric function,
#'   and the upper limit of the summation stops at \eqn{j = b - 1} if \code{b} is an integer.
#'
#'   \item \strong{For non-integer \code{b} or \code{k > 2}}: When \code{b} is non-integer or \code{k} is
#'   greater than 2 the function employs Monte Carlo simulation using the following formula:
#'   \deqn{\text{E}[X^s] \approx \frac{1}{\mathrm{rep}} \sum_{i=1}^{\mathrm{rep}} X_i^s,}
#'   where \eqn{X_i} are the simulated order statistics obtained from the complementary beta distribution.
#'   The method relies on the \code{ros()} function to generate order statistics.
#' }
#'
#' When \code{verbose = TRUE}, the function prints a message only if Monte Carlo simulation is used
#' (i.e., when \code{k > 2} or \code{b} is non-integer).
#'
#' @return
#' The estimated or exact \eqn{k}th moment of the \eqn{r}th order statistic from a complementary beta distribution.
#'
#' @note
#' The closed-form formula is only available for small values of \code{k} and integer \code{b}.
#' Monte Carlo simulation is used otherwise, and results may vary slightly depending on \code{rep}.
#'
#' @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.
#'
#' @seealso
#' \link{ros} for generating random samples of order statistics.
#'
#' @examples
#' # Exact moment when k = 1
#' mo_compbeta(r = 2, n = 15, k = 1, a = 0.5, b = 2)
#'
#' # Simulation when k > 2 or b is non-integer
#' mo_compbeta(r = 2, n = 15, k = 3, a = 2.5, b = 3.7, rep = 1e4)
#'
#' @export
mo_compbeta <- function(r, n, k = 1, a, b, rep = 1e5, seed = 42, verbose = TRUE) {
  if (!all(.is.natural(c(r, n, k)))) stop("invalid arguments")
  if (r > n) stop('"r" must be less than or equal to "n"')
  if (!(a > 0 && b > 0)) stop('both parameters "a" and "b" must be positive')
  if (k == 1) {
    j <- 0:(r - 1)
    term <- choose(n, j) * beta(a + j, b + n - j)
    return(sum(term) / beta(a, b))
  }
  if (k == 2 && b %% 1 == 0) {
    j <- 0:(n - r)
    term <- choose(n - r, j) * (-1)^j * .Ms_recursive(a, b, r + j - 1, k)
    return(sum(term) / beta(r, n - r + 1))
  }
  if (verbose) {
    message("Using Monte Carlo simulation (either b is not integer or k > 2)")
  }
  if (!is.null(seed)) set.seed(seed)
  x <- ros(rep, r, n, qf = ".qcompbeta", a = a, b = b)
  return(mean(x^k))
}




#' Moments of Order Statistics from the Normal Distribution (Simulated)
#'
#' This function computes the moments of order statistics from the normal distribution
#' using simulation.
#'
#' @param r rank of the desired order statistic (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param k order of the moment to compute (default is \code{1}).
#' @param mean mean of the normal distribution (default is \code{0}).
#' @param sd standard deviation of the normal distribution (default is \code{1}).
#' @param rep number of simulations (default is \code{1e5}).
#' @param seed optional seed for random number generation to ensure reproducibility (default is \code{42}).
#'
#' @details
#' This function estimates the \eqn{k}th moment of the \eqn{r}th order statistic in a sample of size \eqn{n}
#' drawn from a normal distribution with the specified mean and standard deviation.
#' The estimation is done via Monte Carlo simulation using the formula:
#'
#' \deqn{\text{E}[X^k] \approx \frac{1}{\mathrm{rep}} \sum_{i=1}^{\mathrm{rep}} X_i^k,}
#' where \eqn{X_i} are the simulated order statistics obtained from the normal distribution.
#'
#' The function relies on the \code{ros()} function to generate order statistics.
#'
#' @return
#' The estimated \eqn{k}th moment of the \eqn{r}th order statistic from a normal distribution.
#'
#' @note
#' The accuracy of the estimated moment depends on the number of simulations (\code{rep}).
#' The default value \code{rep = 1e5} provides a reasonable trade-off between speed and accuracy
#' for most practical cases. For higher order moments or when greater precision is required,
#' users are encouraged to increase \code{rep} (e.g. \code{1e6}).
#'
#' @examples
#' # Compute the first moment (mean) of the 3rd order statistic from a sample of size 10
#' mo_norm(r = 3, n = 10, k = 1, mean = 0, sd = 1)
#'
#' # Compute the second moment of the 2nd order statistic with 1 million simulations
#' mo_norm(r = 2, n = 10, k = 2, rep = 1e6)
#'
#' @seealso \code{\link{ros}}
#'
#' @export
mo_norm <- function(r, n, k = 1, mean = 0, sd = 1, rep = 1e5, seed = 42) {
  if (!all(.is.natural(c(r, n, k)))) stop("invalid arguments")
  if (r > n) stop('"r" must be less than or equal to "n"')
  if (!is.null(seed)) set.seed(seed)
  x <- ros(rep, r, n, dist = "norm", mean = mean, sd = sd)
  return(mean(x^k))
}




#' Moments of Order Statistics from the Student's t-Distribution (Simulated)
#'
#' This function computes the moments of order statistics from the student's t-distribution
#' using simulation.
#'
#' @param r rank of the desired order statistic (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param k order of the moment to compute (default is \code{1}).
#' @param df degrees of freedom for the student's t-distribution.
#' @param rep number of simulations (default is \code{1e5}).
#' @param seed optional seed for random number generation to ensure reproducibility (default is \code{42}).
#'
#' @details
#' This function estimates the \eqn{k}th moment of the \eqn{r}th order statistic in a sample of size \eqn{n}
#' drawn from a student's t-distribution with the specified degrees of freedom (\code{df}).
#'  The estimation is done via Monte Carlo simulation using the formula:
#'
#' \deqn{\text{E}[X^k] \approx \frac{1}{\mathrm{rep}} \sum_{i=1}^{\mathrm{rep}} X_i^k,}
#' where \eqn{X_i} are the simulated order statistics obtained from the student's t-distribution.
#'
#' The function relies on the \code{ros()} function to generate order statistics.
#'
#' @return
#' The estimated \eqn{k}th moment of the \eqn{r}th order statistic from a student's t-distribution.
#'
#' @note
#' The accuracy of the estimated moment depends on the number of simulations (\code{rep}).
#' The default value \code{rep = 1e5} provides a reasonable trade-off between speed and accuracy
#' for most practical cases. For higher order moments or when greater precision is required,
#' users are encouraged to increase \code{rep} (e.g. \code{1e6}).
#'
#' @examples
#' # Compute the first moment (mean) of the 3rd order statistic from a sample of size 10
#' mo_t(r = 3, n = 10, df = 5, k = 1)
#'
#' # Compute the second moment of the 2nd order statistic with 10000 simulations
#' mo_t(r = 2, n = 10, df = 10, k = 2, rep = 1e4)
#'
#' @seealso \code{\link{ros}}
#'
#' @export
mo_t <- function(r, n, k = 1, df, rep = 1e5, seed = 42) {
  if (!all(.is.natural(c(r, n, k)))) stop("invalid arguments")
  if (r > n) stop('"r" must be less than or equal to "n"')
  if (!is.null(seed)) set.seed(seed)
  x <- ros(rep, r, n, dist = "t", df = df)
  return(mean(x^k))
}



#' Moments of Order Statistics from the Beta Distribution (Simulated)
#'
#' This function computes the moments of order statistics from the beta distribution using simulation.
#'
#' @param r rank of the desired order statistic (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param k order of the moment to compute (default is \code{1}).
#' @param a,b non-negative parameters of the beta distribution.
#' @param rep number of simulations (default is \code{1e5}).
#' @param seed optional seed for random number generation to ensure reproducibility (default is \code{42}).
#'
#' @details
#' This function estimates the \eqn{k}th moment of the \eqn{r}th order statistic in a sample of size \eqn{n}
#' drawn from a beta distribution with specified shape parameters. The estimation is done via
#' Monte Carlo simulation using the formula:
#'
#' \deqn{\text{E}[X^k] \approx \frac{1}{\mathrm{rep}} \sum_{i=1}^{\mathrm{rep}} X_i^k,}
#' where \eqn{X_i} are the simulated order statistics from the beta distribution.
#'
#' The function relies on the \code{ros()} function to generate order statistics.
#'
#' @return
#' The estimated \eqn{k}th moment of the \eqn{r}th order statistic from a beta distribution.
#'
#' @note
#' The accuracy of the estimated moment depends on the number of simulations (\code{rep}).
#' The default value \code{rep = 1e5} provides a reasonable trade-off between speed and accuracy
#' for most practical cases. For higher order moments or when greater precision is required,
#' users are encouraged to increase \code{rep} (e.g. \code{1e6}).
#'
#' @seealso
#' \link{ros} for generating random samples of order statistics.
#'
#' @examples
#' # Compute the first moment of the 2nd order statistic from Beta(3, 4) with sample size 5
#' mo_beta(r = 2, n = 5, k = 1, a = 3, b = 4)
#'
#' # Compute the second moment with 10000 simulations
#' mo_beta(r = 2, n = 5, k = 2, a = 2, b = 2.5, rep = 1e4)
#'
#' @export
mo_beta <- function(r, n, k = 1, a, b, rep = 1e5, seed = 42) {
  if (!all(.is.natural(c(r, n, k)))) stop("invalid arguments")
  if (r > n) stop('"r" must be less than or equal to "n"')
  if (!(a > 0 && b > 0)) stop('both parameters "a" and "b" must be positive')
  if (!is.null(seed)) set.seed(seed)
  x <- ros(rep, r, n, dist = "beta", shape1 = a, shape2 = b)
  return(mean(x^k))
}




#' Moments of Order Statistics from the Gamma Distribution (Simulated)
#'
#' This function computes the moments of order statistics from the gamma distribution using simulation.
#'
#' @param r rank of the desired order statistic (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param k order of the moment to compute (default is \code{1}).
#' @param shape shape parameter of the gamma distribution.
#' @param rate rate parameter of the gamma distribution.
#' @param rep number of simulations (default is \code{1e5}).
#' @param seed optional seed for random number generation to ensure reproducibility (default is \code{42}).
#'
#' @details
#' This function estimates the \eqn{k}th moment of the \eqn{r}th order statistic in a sample of size \eqn{n}
#' drawn from a gamma distribution with specified shape and rate parameters. The estimation is done via
#' Monte Carlo simulation using the formula:
#'
#' \deqn{\text{E}[X^k] \approx \frac{1}{\mathrm{rep}} \sum_{i=1}^{\mathrm{rep}} X_i^k,}
#' where \eqn{X_i} are the simulated order statistics from the gamma distribution.
#'
#' The function relies on the \code{ros()} function to generate order statistics.
#'
#' @return
#' The estimated \eqn{k}th moment of the \eqn{r}th order statistic from a gamma distribution.
#'
#' @note
#' The accuracy of the estimated moment depends on the number of simulations (\code{rep}).
#' The default value \code{rep = 1e5} provides a reasonable trade-off between speed and accuracy
#' for most practical cases. For higher order moments or when greater precision is required,
#' users are encouraged to increase \code{rep} (e.g. \code{1e6}).
#'
#' @examples
#' # Compute the first moment (mean) of the 3rd order statistic from a sample of size 10
#' mo_gamma(r = 3, n = 10, shape = 2, rate = 1, k = 1)
#'
#' # Compute the second moment with 10000 simulations
#' mo_gamma(r = 2, n = 10, shape = 2, rate = 0.5, k = 2, rep = 1e4)
#'
#' @seealso \code{\link{ros}}
#'
#' @export
mo_gamma <- function(r, n, k = 1, shape, rate, rep = 1e5, seed = 42) {
  if (!all(.is.natural(c(r, n, k)))) stop("invalid arguments")
  if (r > n) stop('"r" must be less than or equal to "n"')
  if (!is.null(seed)) set.seed(seed)
  x <- ros(rep, r, n, dist = "gamma", shape = shape, rate = rate)
  return(mean(x^k))
}




#' Moments of Order Statistics from the Pareto Distribution (Simulated)
#'
#' This function computes the \eqn{k}th moment of order statistics from the pareto distribution using simulation.
#'
#' @param r rank of the desired order statistic (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param k order of the moment to compute (default is \code{1}).
#' @param scale,shape non-negative parameters of the pareto distribution.
#' @param rep Number of simulations (default is \code{1e5}).
#' @param seed Optional seed for random number generation to ensure reproducibility (default is \code{42}).
#'
#' @details
#' This function estimates the \eqn{k}th moment of the \eqn{r}th order statistic in a sample of size \eqn{n}
#' drawn from a pareto distribution with specified scale and shape parameters. The estimation is done via
#' Monte Carlo simulation using the formula:
#'
#' \deqn{\text{E}[X^k] \approx \frac{1}{\mathrm{rep}} \sum_{i=1}^{\mathrm{rep}} X_i^k,}
#' where \eqn{X_i} are the simulated order statistics from the pareto distribution.
#'
#' The function relies on the \code{ros()} function to generate order statistics.
#'
#' @return
#' The estimated \eqn{k}th moment of the \eqn{r}th order statistic from a pareto distribution.
#'
#' @note
#' The accuracy of the estimated moment depends on the number of simulations (\code{rep}).
#' The default value \code{rep = 1e5} provides a reasonable trade-off between speed and accuracy
#' for most practical cases. For higher order moments or when greater precision is required,
#' users are encouraged to increase \code{rep} (e.g. \code{1e6}).
#'
#' @examples
#' # Compute the first moment (mean) of the 3rd order statistic from a sample of size 10
#' mo_pareto(r = 3, n = 10, scale = 2, shape = 3, k = 1)
#'
#' # Compute the second moment with 1 million simulations
#' mo_pareto(r = 2, n = 10, scale = 1, shape = 2, k = 2, rep = 1e6)
#'
#' @seealso \code{\link{ros}}
#'
#' @export
mo_pareto <- function(r, n, k = 1, scale, shape, rep = 1e5, seed = 42) {
  if (!all(.is.natural(c(r, n, k)))) stop("invalid arguments")
  if (r > n) stop('"r" must be less than or equal to "n"')
  if (!is.null(seed)) set.seed(seed)
  x <- ros(rep, r, n, qf = ".qpareto", scale = scale, shape = shape)
  return(mean(x^k))
}



#' Moments of Order Statistics from the Kumaraswamy Distribution (Simulated)
#'
#' This function computes the moments of order statistics from the kumaraswamy distribution using simulation.
#'
#' @param r rank of the desired order statistic (e.g., \code{1} for the smallest order statistic).
#' @param n sample size from which the order statistic is derived.
#' @param k order of the moment to compute (default is \code{1}).
#' @param a,b positive parameters of the kumaraswamy distribution.
#' @param rep number of simulations (default is \code{1e5}).
#' @param seed optional seed for random number generation to ensure reproducibility (default is \code{42}).
#'
#' @details
#' This function estimates the \eqn{k}th moment of the \eqn{r}th order statistic in a sample of size \eqn{n}
#' drawn from a kumaraswamy distribution with specified shape parameters. The estimation is done via
#' Monte Carlo simulation using the formula:
#'
#' \deqn{\text{E}[X^k] \approx \frac{1}{\mathrm{rep}} \sum_{i=1}^{\mathrm{rep}} X_i^k,}
#' where \eqn{X_i} are the simulated order statistics from the kumaraswamy distribution.
#'
#' The function relies on the \code{ros()} function to generate order statistics.
#'
#' @note
#' The accuracy of the estimated moment depends on the number of simulations (\code{rep}).
#' The default value \code{rep = 1e5} provides a reasonable trade-off between speed and accuracy
#' for most practical cases. For higher order moments or when greater precision is required,
#' users are encouraged to increase \code{rep} (e.g. \code{1e6}).
#'
#' @return
#' The estimated \eqn{k}th moment of the \eqn{r}th order statistic from a kumaraswamy distribution.
#'
#' @seealso
#' \link{ros} for generating random samples of order statistics.
#'
#' @examples
#' # Compute the 2nd moment of the 3rd order statistic from Kumaraswamy(2, 3) with sample size 10
#' mo_kumar(r = 3, n = 10, k = 2, a = 2, b = 3)
#'
#' # Compute the first moment with 10000 simulations
#' mo_kumar(r = 2, n = 5, k = 1, a = 2, b = 2.5, rep = 1e4)
#'
#' @export
mo_kumar <- function(r, n, k = 1, a, b, rep = 1e5, seed = 42) {
  if (!all(.is.natural(c(r, n, k)))) stop("invalid arguments")
  if (r > n) stop('"r" must be less than or equal to "n"')
  if (!(a > 0 && b > 0)) stop('both parameters "a" and "b" must be positive')
  if (!is.null(seed)) set.seed(seed)
  x <- ros(rep, r, n, qf = ".qkumar", a = a, b = b)
  return(mean(x^k))
}







#' Generate Censored Samples (Type I or Type II)
#'
#' This function generates censored samples from a specified distribution, using
#' Type I (time-based) or Type II (failure-based) censoring schemes.
#'
#' @param n total number of items in the sample
#' @param r number of uncensored observations (only for Type II censoring).
#' @param dist a character string specifying the name of the distribution
#'   (e.g., \code{"norm"} for the normal distribution, \code{"exp"} for the exponential distribution).
#' @param type type of censoring: \code{"I"} for Type I (time-based) or \code{"II"} for Type II (failure-based).
#' @param cens.time censoring time for Type I censoring.
#' @param ... further arguments to be passed to \code{dist}.
#'
#' @details
#' This function implements two types of censoring schemes:
#' \enumerate{
#'   \item \strong{Type I censoring:} Observations are censored if they exceed a specified \code{cens.time}.
#'   The function returns all observations less than \code{cens.time}.
#'   \item \strong{Type II censoring:} The smallest \code{r} observations are returned, simulating
#'   a situation where the experiment stops after \code{r} failures.
#' }
#'
#' @return A numeric vector of censored samples.
#'
#' @examples
#' # Type I censoring: Exponential distribution with rate = 1, censored at time 2
#' rcens(n = 10, dist = "exp", type = "I", cens.time = 2, rate = 1)
#'
#' # Type II censoring: Normal distribution, smallest 5 values
#' rcens(n = 10, r = 5, dist = "norm", type = "II", mean = 0, sd = 1)
#'
#' @seealso \code{\link{rpcens2}}
#'
#' @export
rcens <- function(n, r = NULL, dist, type = c("I", "II"), cens.time = NULL, ...) {
  if (!.is.natural(n)) stop("invalid arguments")
  if (!exists(paste0("r", dist))) stop("invalid distribution name, use standard R distributions")
  type <- match.arg(type)

  if (type == "I") {
    if (is.null(cens.time)) stop('argument "cens.time" is missing, with no default')
    s <- do.call(paste0("r", dist), list(n, ...))
    return(s[which(s < cens.time)])
  } else if (type == "II") {
    if (is.null(r)) stop('argument "r" is missing, with no default')
    if (!.is.natural(r)) stop("invalid arguments")
    s <- do.call(paste0("r", dist), list(n, ...))
    return(sort(s)[1:r])
  }
}




#' Generate Progressive Type-II Censored Samples
#'
#' This function generates progressive Type-II censored samples based on the
#' algorithm provided by Balakrishnan and Sandhu (1995).
#'
#' @param n total number of items in the sample.
#' @param R a vector of non-negative integers representing the number of items
#'   censored at each stage of the experiment. The length of \code{R} determines
#'   the number of failure stages (\code{m}), and the sum of \code{R} plus \code{m}
#'   must equal \code{n}.
#' @param dist a character string specifying the name of the distribution
#'   (e.g., \code{"norm"} for the normal distribution, \code{"exp"} for the exponential distribution).
#' @param ... further arguments to be passed to \code{dist}.
#'
#' @details
#' This function implements the algorithm described by Balakrishnan and Sandhu (1995)
#' to simulate progressive Type-II censored samples. Progressive Type-II censoring
#' is a common scheme in reliability and life-testing experiments, where items are
#' progressively removed (censored) during the testing process.
#'
#' @return A numeric vector of size \code{m}, representing the failure times of the observed items.
#'
#' @references
#' Balakrishnan, N., & Sandhu, R. A. (1995). \emph{A simple simulational algorithm for generating progressive
#' Type-II censored samples.} The American Statistician, 49(2), 229-230.
#'
#' @examples
#' # Generate a progressive Type-II censored sample from the normal distribution
#' n <- 10
#' R <- c(2, 1, 2, 0, 0)
#' rpcens2(n, R, dist = "norm", mean = 0, sd = 1)
#'
#' # Generate a progressive Type-II censored sample from the exponential distribution
#' rpcens2(n = 10, R = c(2, 2, 1, 0, 0), dist = "exp", rate = 1)
#'
#' @seealso \code{\link{rcens}}
#'
#' @export
rpcens2 <- function(n, R, dist, ...) {
  if (!.is.natural(n)) stop("invalid arguments")
  if (!all(R >= 0 & floor(R) == R)) stop('"R" must contain non-negative integers')
  m <- length(R)
  if (m + sum(R) != n) stop("invalid arguments")
  if (!exists(paste0("q", dist))) stop("invalid distribution name, use standard R distributions")
  w <- runif(m)
  v <- rep(NA, m)
  u <- v
  for (i in 1:m) v[i] <- w[i]^(1 / (i + sum(R[m:(m - i + 1)])))
  for (i in 1:m) u[i] <- 1 - prod(v[m:(m - i + 1)])
  x <- do.call(paste0("q", dist), list(u, ...))
  return(x)
}




#' Generate Upper and Lower k-Records from Continuous Distributions
#'
#' This function generates \eqn{k}-records (upper or lower) from a specified continuous distribution.
#'
#' Note: Setting \code{k = 1} generates standard (1-)records.
#'
#' @param size number of k-records to generate.
#' @param k the rank of the record to generate (\eqn{k}-record).
#' @param record the type of record to generate: \code{"upper"} for upper k-records,
#'   \code{"lower"} for lower k-records. Default is \code{"upper"}.
#' @param dist a character string specifying the name of the continuous distribution
#'   (e.g., \code{"norm"}, \code{"exp"}, \code{"gamma"}).
#' @param ... further arguments to be passed to \code{dist}.
#'
#' @return A numeric vector of size \code{size}, representing the simulated k-records.
#'
#' @examples
#' # Generate 5 upper 2-records from the normal distribution
#' rkrec(size = 5, k = 2, record = "upper", dist = "norm", mean = 0, sd = 1)
#'
#' # Generate 5 lower 3-records from the exponential distribution
#' rkrec(size = 5, k = 3, record = "lower", dist = "exp", rate = 1)
#'
#' @seealso \code{\link{ros}}
#'
#' @export
rkrec <- function(size, k, record = c("upper", "lower"), dist, ...) {
  if (!all(.is.natural(c(size, k)))) {
    stop('"size" and "k" must be natural values')
  }
  record <- match.arg(record)
  if (!exists(paste0("r", dist))) stop("invalid distribution name, use standard R continuous distributions")
  if (record == "upper") {
    return(ros(size, 1, k, dist, ...))
  } else if (record == "lower") {
    return(ros(size, k, k, dist, ...))
  }
}