R/densityExpressions.R

In OBRE: Optimal B-Robust Estimator Tools

#' Distributions formulas for OBRE
#'
#' Function containing expressions of density and cumulative functions, plus the first and second derivatives.
#'
#' @importFrom stats D
#' @importFrom methods is
#' @param strDistribution Distribution input between "normal" (Normal distribution), "logNormal" (logNormal distribution),
#' "weibull" (Weibull distribution), "logLogistic" (logLogistic distribution), "gpd2" (Generalized Pareto
#' Distribution with two parameters) or "custom" if the distribution is written by the user.
#' @param eDensityFun The density of a two parameters distribution. This should be an expression object, the
#' two parameters should be called "nTheta1" and "nTheta2", the data "nvData" and its formula should be derivable
#'
#' @return Returns list containing all the symbolic functions.
#'
#' @export
#' @examples
#' # Generates the Normal distribution input for OBRE
#' distrForOBRE <- densityExpressions(strDistribution = "normal")
#' # The same result can be generated by inserting manually the formula
#' distrForOBRE <- densityExpressions(strDistribution = "custom",
#' eDensityFun = expression((exp( -((nvData - nTheta1)^2) / (2 * nTheta2^2)) /
#' (sqrt(2 * pi) * nTheta2))))
#'

densityExpressions = function(strDistribution = "normal", eDensityFun = NA) {

# Preliminary controls on the distribution name ####
  if (strDistribution == "custom") {
    if (!is(eDensityFun, "expression")) {
      cat("The input is not an expression object.")
      return("Input error.")
    } else {
      nNumberOfParameters = 0
      nControl = TRUE
      while (nControl == TRUE) {
        nControl = grepl(paste0("nTheta", (nNumberOfParameters + 1)), as.character(x = eDensityFun))
        nNumberOfParameters = nNumberOfParameters + 1
      }
      nNumberOfParameters = nNumberOfParameters - 1
      nControlData = grepl("nvData", as.character(x = eDensityFun))
      if (nNumberOfParameters != 2 ||  nControlData == FALSE) {
        cat("The input is not an admissible expression object.")
        return("Input error.")
      }
    }
  } else {
    if (length(which(strDistribution %in% c("normal", "logNormal", "weibull", "logLogistic", "gpd2"))) == 0) {
      cat("The input is not an admissible distribution.")
      return("Input error.")
    }
  }

# Density expression depending on the current distribution ####
  eDerivDensityFunTheta = expression()
  switch(strDistribution,
         "normal" = {
           eDensityFun = expression((exp( -((nvData - nTheta1)^2) / (2 * nTheta2^2)) / (sqrt(2 * pi) * nTheta2)))
           nNumberOfParameters = 2
         },
         "logNormal" = {
           eDensityFun = expression((exp( -((log(nvData) - nTheta1)^2) / (2 * nTheta2^2)) /
                                       (sqrt(2 * pi) * nTheta2 * nvData)))
           nNumberOfParameters = 2
         },
         "weibull" = {
           eDensityFun = expression(((nTheta1 / nTheta2) * ((nvData / nTheta2)^(nTheta1 - 1)) *
                                       exp(-((nvData / nTheta2)^nTheta1))))
           nNumberOfParameters = 2
         },
         "logLogistic" = {
           eDensityFun = expression(((nTheta1 * ((nvData / nTheta2)^nTheta1)) / (nvData*((1 + ((nvData / nTheta2)^nTheta1))^2))))
           nNumberOfParameters = 2
         },
         "gpd2" = {
           eDensityFun = expression(((1 / nTheta2) * ((1 + (nTheta1 * (nvData)) / nTheta2)^(-1 / nTheta1 - 1))))
           nNumberOfParameters = 2
         })

# First derivatives ####
  eDerivDensityFunTheta[1] = as.expression(D(eDensityFun, "nTheta1"))
  eDerivDensityFunTheta[2] = as.expression(D(eDensityFun, "nTheta2"))

# Second derivatives ####
  eDeriv2DensityFunTheta1Theta1 = as.expression(D(eDerivDensityFunTheta[1], "nTheta1"))
  eDeriv2DensityFunTheta2Theta2 = as.expression(D(eDerivDensityFunTheta[2], "nTheta2"))
  eDeriv2DensityFunTheta1Theta2 = as.expression(D(eDerivDensityFunTheta[1], "nTheta2"))
  eDeriv2DensityFunTheta2Theta1 = as.expression(D(eDerivDensityFunTheta[2], "nTheta1"))

# Output ####
  lExpressionsDist = list(eDensityFun = eDensityFun,
                          eDerivDensityFunTheta = eDerivDensityFunTheta,
                          eDeriv2DensityFunTheta1Theta1 = eDeriv2DensityFunTheta1Theta1,
                          eDeriv2DensityFunTheta2Theta2 = eDeriv2DensityFunTheta2Theta2,
                          eDeriv2DensityFunTheta1Theta2 = eDeriv2DensityFunTheta1Theta2,
                          eDeriv2DensityFunTheta2Theta1 = eDeriv2DensityFunTheta2Theta1,
                          nNumberOfParameters = nNumberOfParameters)
  class(lExpressionsDist) = "OBREdist"
  return(lExpressionsDist)
}