R/qcnb.R

In cnbdistr: Conditional Negative Binomial Distribution

#'
#' @title
#' Quantile Function of Conditional Negative Binomial
#'
#'
#' @description
#' Quantile function of the conditional distribution of X given X + Y = D,
#' where X ~ NB(r1, p1) and Y ~ NB(r2, p2) are drawn from two negative binomials,
#' independent of each other,
#' and assuming p1/p2 = lambda.
#'
#' @details
#' Need to specify full list of arguments, as default values have not been set.
#'
#' @author Xiaotian Zhu, \email{xiaotian.zhu.psualum@@gmail.com}
#'
#' @param p a nonempty vector of probabilities (0 <= p[i] <= 1 for all i).
#' @param D a positive integer.
#' @param r1 a positive value.
#' @param r2 a positive value.
#' @param lambda a positive value.
#'
#' @return A vector x such that x[i] = Inf\{x in 0:D, p[i] <= Pr(X <= x | X + Y = D)\} for all i.
#'
#' @examples
#' qcnb(0.035193, 2000, 120, 90, 0.994)
#' qcnb(seq(0, 1, 0.05), 7, 2, 0.4, 0.6)
#'
#' @seealso
#' \code{\link{dcnb}, \link{pcnb}, \link{rcnb}.}
#'
#'
#' @export
"qcnb" <- function(p, D, r1, r2, lambda){

  if (missing(p) || missing(D) || missing(r1) || missing(r2) || missing(lambda))
    stop("Need to specify a full set of arguments: p, D, r1, r2, lambda.")

  if (!is.numeric(D) || length(D)!=1 || !D%%1==0 || D<=0)
    stop("D needs to be a positive integer.")

  if (!is.numeric(p) || length(p)<1 || any(p<0) || any(p>1))
    stop("p needs to be a nonempty vector of probabilities (0 <= p[i] <= 1 for all i).")

  if (!is.numeric(r1) || length(r1)!=1 || r1<=0)
    stop("r1 needs to be a positive value.")

  if (!is.numeric(r2) || length(r2)!=1 || r2<=0)
    stop("r2 needs to be a positive value.")

  if (!is.numeric(lambda) || length(lambda)!=1 || lambda<=0)
    stop("lambda needs to be a positive value.")


  return_value <- rep(0, length(p))

  for (i in 1:length(p)) {

    if (p[i] <= dcnb(0, D, r1, r2, lambda)){
      return_value[i] <- 0
    } else if (p[i] == 1) {
      return_value[i] <- D
    } else {

      tmplower <- 0
      tmpupper <- D

      while ((tmpupper - tmplower) > 1) {

        tmptest <- floor((tmplower + tmpupper)/2)

        if (pcnb(tmptest, D, r1, r2, lambda) < p[i]) {
          tmplower <- tmptest
        } else {
          tmpupper <- tmptest
        }

      }

      return_value[i] <- tmpupper

    }

  }


  return_value

}