# R/gbeta.R In gbeta: Generalized Beta and Beta Prime Distributions

#### Documented in dgbetapgbetaqgbetargbeta

#' Generalized Beta distribution
#' @description Density, distribution function, quantile function, and random
#'   generation for the generalized Beta distribution.
#'
#' @param u numeric vector
#' @param q numeric vector of quantiles
#' @param p numeric vector of probabilities
#' @param n positive integer, the desired number of simulations
#' @param c,d,kappa,tau parameters; they must be strictly positive numbers,
#'   except \code{kappa} which can take any value
#' @param log logical, whether to return the log-density
#' @param method the method of random generation, \code{"mixture"} or
#'   \code{"arou"}; only a positive \code{kappa} is allowed for the
#'   \code{"mixture"} method, but this method is faster
#'
#' @references
#' \itemize{
#'   \item Marwa Hamza & Pierre Vallois.
#'     \emph{On Kummerâ€™s distributions of type two and generalized Beta
#'           distributions}.
#'     Statistics & Probability Letters 118 (2016), pp. 60-69.
#'     <doi:10.1016/j.spl.2016.03.014>
#'   \item James J. Chen & Melvin R. Novick.
#'     \emph{Bayesian Analysis for Binomial Models with Generalized Beta Prior
#'           Distributions}.
#'     Journal of Educational Statistics 9, No. 2 (1984), pp. 163-175.
#'     <doi:10.3102/10769986009002163>
#' }
#'
#' @examples library(gbeta)
#' curve(dgbeta(x, 4, 12, 10, 0.01), axes = FALSE, lwd = 2)
#' axis(1)
#'
#' @importFrom stats qbeta
#' @importFrom Runuran uq
#'
#' @rdname GBeta
#' @name GBeta
#' @export
dgbeta <- function(u, c, d, kappa, tau, log = FALSE){
stopifnot(c > 0, d > 0, tau > 0)
out <- numeric(length(u))
outside01 <- u < 0 | u > 1
out[outside01] <- ifelse(log, -Inf, 0)
inside01 <- !outside01
if(any(inside01)){
u <- u[inside01]
out[inside01] <- if(log){
-2 * log(1-u) + dgbetap(u/(1-u), c, d, kappa, tau, scale = 1, log = TRUE)
}else{
1/(1-u)^2 * dgbetap(u/(1-u), c, d, kappa, tau, scale = 1)
}
}
out
}

#' @rdname GBeta
#' @export
pgbeta <- function(q, c, d, kappa, tau){
stopifnot(c > 0, d > 0, tau > 0)
out <- numeric(length(q))
negs <- q <= 0
sup1 <- q >= 1
out[sup1] <- 1
in01 <- !(negs || sup1)
if(any(in01)){
q <- q[in01]
out[in01] <- pgbetap(q/(1-q), c, d, kappa, tau, scale = 1)
}
out
}

#' @rdname GBeta
#' @export
rgbeta <- function(n, c, d, kappa, tau, method = "mixture"){
betap <- rgbetap(n, c, d, kappa, tau, scale = 1, method = method)
betap / (1 + betap)
}

#' @rdname GBeta
#' @export
qgbeta <- function(p, c, d, kappa, tau){
stopifnot(all(p >= 0 & p <= 1))
stopifnot(c > 0, d > 0, tau > 0)
if(kappa == 0 || tau == 1){
return(qbeta(p, c, d))
}
uq(unuran_gbeta_pinv(c, d, kappa, tau), p)
}


## Try the gbeta package in your browser

Any scripts or data that you put into this service are public.

gbeta documentation built on Nov. 19, 2020, 9:07 a.m.