R/dcomper.R

In dsfa: Distributional Stochastic Frontier Analysis

#' Composed-Error distribution
#'
#' Probablitiy density function, distribution, quantile function and random number generation for the composed-error distribution
#'
#' @return \code{dcomper()} gives the density, \code{pcomper()} give the distribution function, \code{qcomper()} gives the quantile function, and \code{rcomper()} generates random numbers, with given parameters.
#' \code{dcomper()} and \code{pcomper()} returns a \code{derivs} object.
#'
#' @details  This is wrapper function for the normal-halfnormal and normal-exponential distribution. A random variable \eqn{X} follows a composed error distribution if \eqn{X = V + s \cdot U }, where \eqn{V \sim N(\mu, \sigma_V^2)} and \eqn{U \sim HN(0,\sigma_U^2)} or \eqn{U \sim Exp(\sigma_U^2)}.
#' For more details see \code{\link{dnormhnorm}} and \code{\link{dnormexp}}. Here, \eqn{s=-1} for production and \eqn{s=1} for cost function.
#' 
#'
#' @param x numeric vector of quantiles.
#' @param mu numeric vector of \eqn{\mu}.
#' @param sigma_v numeric vector of \eqn{\sigma_V}. Must be positive.
#' @param sigma_u numeric vector of \eqn{\sigma_U}. Must be positive.
#' @param s integer; \eqn{s=-1} for production and \eqn{s=1} for cost function.
#' @param tri optional; index matrix for upper triangular, generated by \code{\link{trind_generator}}.
#' @param log.p logical; if TRUE, probabilities p are given as log(p).
#' @param distr string; determines the distribution:\cr
#' `normhnorm`, Normal-halfnormal distribution \cr
#' `normexp`, Normal-exponential distribution \cr
#' @inheritParams list2derivs
#' 
#' @examples
#' pdf <- dcomper(x=5, mu=1, sigma_v=2, sigma_u=3, s=-1, distr="normhnorm")
#' cdf <- pcomper(q=5, mu=1, sigma_v=2, sigma_u=3, s=-1, distr="normhnorm")
#' q <- qcomper(p=seq(0.1, 0.9, by=0.1), mu=1, sigma_v=2, sigma_u=3, s=-1, distr="normhnorm")
#' r <- rcomper(n=10, mu=1, sigma_v=2, sigma_u=3, s=-1, distr="normhnorm")
#'
#' @references
#' \itemize{
#' \item \insertRef{aigner1977formulation}{dsfa}
#' \item \insertRef{kumbhakar2015practitioner}{dsfa}
#' \item \insertRef{schmidt2020analytic}{dsfa}
#' \item \insertRef{gradshteyn2014table}{dsfa}
#' \item \insertRef{azzalini2013skew}{dsfa}
#' }
#' 
#' @family distribution
#' 
#' @export
dcomper <- function(x, mu=0, sigma_v=1, sigma_u=1, s=-1, distr="normhnorm", deriv_order=0, tri=NULL, log.p = FALSE){
  #Density function of the comper distribution
  if (any(sigma_v <= 0)|any(sigma_u <= 0)){
    stop(paste("sigma_v and sigma_u must be positive", "\n", ""))
  }
  
  if (any(!s%in%c(-1,1))){
    stop(paste("s must be {-1, 1}", "\n", ""))
  }
  
  distr<-match.arg(distr,c("normhnorm","normexp"))
  
  X<-tryCatch(cbind(x, mu, sigma_v, sigma_u), warning=function(w) {
    stop("Input vectors have incompatible lengths")})
  
  if(is.null(tri)){
    tri=trind_generator(3)
  }
  
  out<-dcomper_cpp(x=X[,1, drop=T], m=X[,2, drop=T], v=X[,3, drop=T], u=X[,4, drop=T], s=s, distr=distr, deriv_order=deriv_order, tri=tri, logp=log.p)
  
  #Return ouptut
  return(out)
}

#' @describeIn dcomper distribution function for the composed-error distribution.
#' @param p numeric vector of probabilities.
#' @export
pcomper <- function(q, mu=0, sigma_v=1, sigma_u=1, s=-1, distr="normhnorm", deriv_order=0, tri=NULL, log.p = FALSE){
  #Probability function of the comper distribution
  if (any(sigma_v <= 0)|any(sigma_u <= 0)){
    stop(paste("sigma_v and sigma_u must be positive", "\n", ""))
  }
  
  if (any(!s%in%c(-1,1))){
    stop(paste("s must be {-1, 1}", "\n", ""))
  }
  
  distr<-match.arg(distr,c("normhnorm","normexp"))
  
  X<-tryCatch(cbind(q, mu, sigma_v, sigma_u), warning=function(w) {
    stop("Input vectors have incompatible lengths")})
  
  if(is.null(tri)){
    tri=trind_generator(3)
  }
  
  out<-pcomper_cpp(q=X[,1, drop=T], m=X[,2, drop=T], v=X[,3, drop=T], u=X[,4, drop=T], s=s, distr=distr, deriv_order=deriv_order, tri=tri, logp=log.p)
  
  #Return ouptut
  return(out)
}

#' @describeIn dcomper quantile function for the composed-error distribution.
#' @param q numeric vector of quantiles.
#' @export
qcomper <- function(p, mu=0, sigma_v=1, sigma_u=1, s=-1, distr="normhnorm", log.p = FALSE){
  #Quantile function of the comper distribution
  distr<-match.arg(distr,c("normhnorm","normexp"))
  
  if(distr=="normhnorm"){
    out<-qnormhnorm(p, mu, sigma_v, sigma_u, s)
  }
  
  if(distr=="normexp"){
    out<-qnormexp(p, mu, sigma_v, sigma_u, s)
  }
  
  #Return output
  return(out)
}

#' @describeIn dcomper random number generation for the composed-error distribution.
#' @param n positive integer; number of observations.
#' @export
rcomper <- function(n, mu=0, sigma_v=1, sigma_u=1, s=-1, distr="normhnorm"){
  #Function to generate n random numbers from the comper distribution
  distr<-match.arg(distr,c("normhnorm","normexp"))
  
  if(distr=="normhnorm"){
    out<-rnormhnorm(n, mu, sigma_v, sigma_u, s)
  }
  
  if(distr=="normexp"){
    out<-rnormexp(n, mu, sigma_v, sigma_u, s)
  }
  
  #Return output
  return(out)
}