R/distrMode.R

Defines functions distrMode betaMode cauchyMode chisqMode dagumMode expMode fMode fiskMode frechetMode gammaMode normMode gevMode ghMode ghtMode gldMode gompertzMode gpdMode gumbelMode hypMode koenkerMode kumarMode laplaceMode logisMode lnormMode lomaxMode maxwellMode mvnormMode nakaMode nigMode paralogisticMode paretoMode rayleighMode stableMode stableMode2 tMode unifMode weibullMode yulesMode bernMode binomMode geomMode hyperMode nbinomMode poisMode

Documented in bernMode betaMode binomMode cauchyMode chisqMode dagumMode distrMode expMode fiskMode fMode frechetMode gammaMode geomMode gevMode ghMode ghtMode gldMode gompertzMode gpdMode gumbelMode hyperMode hypMode koenkerMode kumarMode laplaceMode lnormMode logisMode lomaxMode maxwellMode mvnormMode nakaMode nbinomMode nigMode normMode paralogisticMode paretoMode poisMode rayleighMode stableMode stableMode2 tMode unifMode weibullMode yulesMode

#' @title 
#' Mode of some continuous and discrete distributions
#' 
#' @description 
#' These functions return the mode of the main probability 
#' distributions implemented in R.
#' 
#' @note 
#' Some functions like \code{normMode} or \code{cauchyMode}, which relate 
#' to symmetric distributions, are trivial, but are implemented for the sake of 
#' exhaustivity.
#' 
#' @param x
#' character. The name of the distribution to consider. 
#' 
#' @param ...
#' Additional parameters. 
#' 
#' @return 
#' A numeric value is returned, the (true) mode of the distribution.
#' 
#' @author 
#' \code{\link[fBasics]{ghMode}} and \code{\link[fBasics]{ghtMode}} are from 
#' package \pkg{fBasics}; 
#' \code{\link[fBasics]{hypMode}} was written by David Scott; 
#' \code{\link[fBasics]{gldMode}}, \code{\link[fBasics]{nigMode}} and 
#' \code{\link[stabledist]{stableMode}} were written by Diethelm Wuertz.
#' 
#' @seealso 
#' \code{\link[modeest]{mlv}} for the estimation of the mode; 
#' the documentation of the related distributions 
#' \code{\link[stats]{Beta}}, \code{\link[stats]{GammaDist}}, etc.
#' 
#' @importFrom statip name2distr
#' @export
#' 
distrMode <- 
function(x, 
         ...)
{
  stopifnot(is.character(x))
  x <- match.arg(statip::name2distr(x), .distributionsList())
  do.call(paste0(x, "Mode"), list(...))
}


# Beta distribution
#' @inheritParams stats::dbeta
#' 
#' @export
#' @rdname distrMode
#' 
#' @examples 
#' ## Beta distribution
#' curve(dbeta(x, shape1 = 2, shape2 = 3.1), 
#'       xlim = c(0,1), ylab = "Beta density")
#' M <- betaMode(shape1 = 2, shape2 = 3.1)
#' abline(v = M, col = 2)
#' mlv("beta", shape1 = 2, shape2 = 3.1)
#' 
betaMode <-
function(shape1,
         shape2,
         ncp = 0)
{
  if (ncp == 0) {
    M <- (shape1-1)/(shape1+shape2-2)
  } else {
    warning("still to be done. 'NA' is returned", 
            call. = FALSE)
    M <- NA
  }
  M
}


# Cauchy distribution
#' @inheritParams stats::dcauchy
#' 
#' @export
#' @rdname distrMode
#' 
cauchyMode <-
function(location = 0,
         ...)
{
  location
}


# Chi-square distribution
#' @inheritParams stats::dchisq
#' 
#' @export
#' @rdname distrMode
#' 
chisqMode <-
function(df,
         ncp = 0)
{
  if (ncp == 0) {
    M <- max(df-2, 0)
  } else {
    warning("still to be done. 'NA' is returned", 
            call. = FALSE)
    M <- NA
  }
  M
}


# Dagum distribution
#' @inheritParams VGAM::ddagum
#' 
#' @export
#' @rdname distrMode
#' 
dagumMode <- 
function(scale = 1, shape1.a, shape2.p)
{
  scale*((shape1.a*shape2.p-1)/(shape1.a+1))^(1/shape1.a)
}


# Exponential distribution
#' @inheritParams stats::dexp
#' 
#' @export
#' @rdname distrMode
#' 
expMode <-
function(...)
{
  0
}


# F distribution
#' @inheritParams stats::df
#' 
#' @export
#' @rdname distrMode
#' 
fMode <-
function(df1,
         df2)
{ 
  if (df1 > 2) {
    M <- (1-2/df1)*(df2/(2+df2))
  } else {
    warning("still to be done. 'NA' is returned", 
            call. = FALSE)
    M <- NA
  }
  M
}


# Fisk distribution
#' @inheritParams VGAM::dfisk
#' 
#' @export
#' @rdname distrMode
#' 
fiskMode <- 
function(scale = 1, 
         shape1.a)
{
  scale*((shape1.a-1)/(shape1.a+1))^(1/shape1.a)
}


# Frechet distribution
#' @inheritParams VGAM::dfrechet
#' 
#' @export
#' @rdname distrMode
#' 
frechetMode <-
function(location = 0,
         scale = 1,
         shape = 1,
         ...)
{
  location + scale*(shape/(1+shape))^(1/shape)
}


# Gamma distribution
#' @inheritParams stats::dgamma
#' 
#' @export
#' @rdname distrMode
#' 
gammaMode <-
function(shape,
         rate = 1,
         scale = 1/rate)
{
  scale*(shape-1)
}


# Gaussian (normal) distribution
#' @inheritParams stats::dnorm
#' 
#' @export
#' @rdname distrMode
#' 
normMode <-
function(mean = 0,
         ...)
{
  mean
}


# Generalized extreme value distribution
#' @inheritParams VGAM::dgev
#' 
#' @export
#' @rdname distrMode
#' 
gevMode <-
function(location = 0,
         scale = 1,
         shape = 0,
         ...)
{
  k <- pmax(0,(1+shape))^(-shape)-1
  shape[shape==0] <- Inf
  location + (scale/shape)*k
}


# Generalized hyperbolic distribution
#' @inheritParams fBasics::ghMode
#' 
#' @importFrom fBasics ghMode
#' @export
#' @rdname distrMode
#' 
ghMode <-
function(alpha = 1,
         beta = 0,
         delta = 1,
         mu = 0,
         lambda = -1/2)
{
  fBasics::ghMode(alpha, beta, delta, mu, lambda)
}


# Generalized Hyperbolic Student-t
#' @inheritParams fBasics::ghtMode
#' 
#' @importFrom fBasics ghtMode
#' @export
#' @rdname distrMode
#' 
ghtMode <-
function(beta = 0.1, 
         delta = 1, 
         mu = 0, 
         nu = 10) 
{
  fBasics::ghtMode(beta, delta, mu, nu)
}


# Generalized Lambda Distribution
#' @inheritParams fBasics::gldMode
#' 
#' @importFrom fBasics gldMode
#' @export
#' @rdname distrMode
#' 
gldMode <-
function(lambda1 = 0, 
         lambda2 = -1, 
         lambda3 = -1/8, 
         lambda4 = -1/8) 
{
  fBasics::gldMode(lambda1, lambda2, lambda3, lambda4)
}


# Gompertz distribution
#' @inheritParams VGAM::dgompertz
#' 
#' @export
#' @rdname distrMode
#' 
gompertzMode <-
function(scale = 1, 
         shape)
{
  if (shape < scale) {
    M <- log(scale/shape)/scale
  } else {
    M <- 0
  }
  M
}


# Generalized Pareto distribution
#' @inheritParams VGAM::dgpd
#' 
#' @export
#' @rdname distrMode
#' 
gpdMode <-
function(location = 0,
         scale = 1,
         shape = 0)
{
  if (shape == -1) {
    warning("all values between 'loc' and 'loc+scale' are modes, only the mean value is returned", 
            call. = FALSE)
    M <- location + scale/2  
  } else if (-2-1/shape > 0) {
    M <- location - scale/shape
  } else {
    warning("the density is not continuous at the mode",
            call. = FALSE)
    M <- location
  }
  M
}


# Gumbel distribution
#' @inheritParams VGAM::dgumbel
#' 
#' @export
#' @rdname distrMode
#' 
gumbelMode <-
function(location = 0,
         ...)
{
  location
}


# Hyperbolic distribution
#' @inheritParams fBasics::hypMode
#' 
#' @importFrom fBasics hypMode
#' @export
#' @rdname distrMode
#' 
hypMode <-
function(alpha = 1,
         beta = 0,
         delta = 1,
         mu = 0,
         pm = c(1, 2, 3, 4)) 
{
  fBasics::hypMode(alpha, beta, delta, mu, pm)
}


# Koenker distribution 
#' @inheritParams stats::dcauchy
#'
#' @export
#' @rdname distrMode
#' 
koenkerMode <- 
function(location = 0, 
         ...)
{
  location
}


# Kumaraswamy distribution
#' @inheritParams VGAM::dkumar
#' 
#' @export
#' @rdname distrMode
#' 
kumarMode <-
function(shape1,
         shape2)
{
  (shape1-1)/(shape1*shape2 - 1)^(1/shape1)
}


# Laplace distribution
#' @inheritParams VGAM::dlaplace
#' 
#' @export
#' @rdname distrMode
#' 
laplaceMode <-
function(location = 0,
         ...)
{
  location
}


# Logistic distribution
#' @inheritParams stats::dlogis
#' 
#' @export
#' @rdname distrMode
#' 
logisMode <-
function(location = 0,
         ...)
{
  location
}


# Lognormal distribution
#' @inheritParams stats::dlnorm
#' 
#' @export
#' @rdname distrMode
#' 
#' @examples 
#' ## Lognormal distribution
#' curve(stats::dlnorm(x, meanlog = 3, sdlog = 1.1), 
#'       xlim = c(0, 10), ylab = "Lognormal density")
#' M <- lnormMode(meanlog = 3, sdlog = 1.1)
#' abline(v = M, col = 2)
#' mlv("lnorm", meanlog = 3, sdlog = 1.1)
#' 
lnormMode <-
function(meanlog = 0,
         sdlog = 1)
{
  exp(meanlog - sdlog^2)
}


# Lomax distribution
#' @inheritParams VGAM::dlomax
#' 
#' @export
#' @rdname distrMode
#' 
lomaxMode <- 
function(...)
{
  0
}


# Maxwell-Boltzmann distribution
#' @inheritParams VGAM::dmaxwell
#' 
#' @export
#' @rdname distrMode
#' 
maxwellMode <-
function(rate)
{
  sqrt(2/rate)
}


# Multivariate normal distribution
#' @inheritParams mvtnorm::dmvnorm
#' 
#' @export
#' @rdname distrMode
#' 
mvnormMode <- 
function(mean, 
         ...)
{
  mean
}


# Nakagami distribution
#' @inheritParams VGAM::dnaka
#' 
#' @export
#' @rdname distrMode
#'
nakaMode <- 
function(scale = 1, 
         shape)
{
  sqrt(2*(2*shape - 1)*scale/(4*shape))
}


# Normal Inverse Gaussian distribution
#' @inheritParams fBasics::nigMode
#' 
#' @importFrom fBasics nigMode
#' @export
#' @rdname distrMode
#' 
nigMode <-
function(alpha = 1,
         beta = 0,
         delta = 1,
         mu = 0)
{
  fBasics::nigMode(alpha, beta, delta, mu)
}


# Paralogistic distribution
#' @inheritParams VGAM::dparalogistic
#' 
#' @export
#' @rdname distrMode
#' 
paralogisticMode <- 
function(scale = 1, 
         shape1.a)
{
  if (shape1.a <= 1) {
    warning("the density is not continuous at the mode", 
            call. = FALSE)
    0
  } else {
    scale*((shape1.a-1)/(shape1.a^2+1))^(1/shape1.a)
  }
}


# Pareto distribution
#' @inheritParams VGAM::dpareto
#' 
#' @export
#' @rdname distrMode
#' 
#' @examples 
#' curve(VGAM::dpareto(x, scale = 1, shape = 1), xlim = c(0, 10))
#' abline(v = paretoMode(scale = 1), col = 2)
#' 
paretoMode <-
function(scale = 1,
         ...)
{
  warning("the density is not continuous at the mode", 
          call. = FALSE)  
  scale
}


# Rayleigh distribution       
#' @inheritParams VGAM::drayleigh
#' 
#' @export
#' @rdname distrMode
#' 
rayleighMode <-
function(scale = 1)
{
  scale
}


# Stable distribution
#' @inheritParams stabledist::dstable
#' 
#' @importFrom stabledist dstable
#' @importFrom stats optimize
#' @export
#' @rdname distrMode
#' 
stableMode <-
function(alpha,
         beta,
         gamma = 1,
         delta = 0,
         pm = 0, 
         ...) 
{
  beta.max <- 1 - 1e-11
  tol <- .Machine$double.eps^0.25
  if (gamma == 1 & delta == 0 & pm == 0) {
    stopifnot(0 < alpha, alpha <= 2, length(alpha) == 1, -1 <= beta, beta <= 1, length(beta) == 1, length(beta.max) == 1)
    if (alpha * beta == 0) {
      M <- 0
    }
    if (beta > beta.max) {
      beta <- beta.max
    }
    M <- stats::optimize(stabledist::dstable, interval = c(-0.7, 0) * sign(beta), alpha = alpha, 
                         beta = beta, pm = 0, maximum = TRUE, tol = tol)$maximum    
    return(M)
  } else {
    warning("still to be done. 'NA' is returned", 
            call. = FALSE)
    return(NA)
  }
}


# from package 'stable'
#' @inheritParams stable::stable.mode
#' 
#' @importFrom stable stable.mode
#' @export
#' @rdname distrMode
#' 
stableMode2 <- 
function(loc, 
         disp, 
         skew, 
         tail)
{
  stable::stable.mode(loc, disp, skew, tail)
}


# Weibull distribution (in the context of extreme value theory)
# 
# #' @export
# #' @rdname distrMode
# #' 
# nweibullMode <-
# function(loc = 0,
#          scale = 1,
#          shape = 1)
# {
#   #return(loc + (scale/-abs(shape))*((1-abs(shape))^(abs(shape))-1))
#   if (shape < 1) {
#     M <- loc
#   } else {
#     M <- loc - scale*((shape-1)/shape)^(1/shape)
#   }
#   M
# }


# Student distribution
#' @inheritParams stats::dt
#' 
#' @export
#' @rdname distrMode
#' 
tMode <-
function(df,
         ncp)
{
  if (ncp == 0) {
    M <- 0
  } else {
    warning("still to be done. 'NA' is returned", 
            call. = FALSE)
    M <- NA
  }
  M
}


# Uniform distribution
#' @inheritParams stats::dunif
#' 
#' @export
#' @rdname distrMode
#' 
unifMode <-
function(min = 0,
         max = 1)
{
  warning("all values between 'min' and 'max' are modes, only the mean value is returned", 
          call. = FALSE)
  (min+max)/2
}


# Weibull distribution
#' @inheritParams stats::dweibull
#' 
#' @export
#' @rdname distrMode
#' 
weibullMode <-
function(shape,
         scale = 1)
{
  scale*(1-1/shape)^(1/shape)
}


# Yules-Simon distribution
#' @inheritParams VGAM::dyules
#' 
#' @export
#' @rdname distrMode
#' 
yulesMode <- 
function(...)
{
  1
}


## Discrete distributions
#------------------------------------------------------

# Bernoulli distribution
#' @inheritParams statip::dbern
#' 
#' @export
#' @rdname distrMode
#' 
bernMode <-
function(prob)
{
  if (prob > 1 || prob < 0) return(NaN)
  q <- 1 - prob
  if (q > prob) return(0)
  if (q == prob) {
    c(0,1)
  } else {
    1
  }
}


# Binomial distribution
#' @inheritParams stats::dbinom
#' 
#' @export
#' @rdname distrMode
#' 
binomMode <-
function(size,
         prob)
{
  if (prob > 1 || prob < 0) return(NaN)
  if (prob == 0) {
    return(0)
  } else {
    if (prob == 1) {
      return(size)
    } else {
      x <- ceiling((size+1)*prob - 1)
      if (x == (size+1)*prob - 1) {
        return(c(x, x+1))
      } else {
        return(x)
      }      
    }
  }     
}


# Geometric distribution
#' @inheritParams stats::dgeom
#' 
#' @export
#' @rdname distrMode
#' 
geomMode <-
function(...)
{
  1
}


# Hypergeometric distribution
#' @inheritParams stats::dhyper
#' 
#' @export
#' @rdname distrMode
#' 
hyperMode <-
function(m,
         n,
         k,
         ...)
{
  lambda <- (m+1)*(k+1)/(m+n+1)  
  if (lambda == 0) return(0)
  x <- floor(lambda)
  if (lambda == x) {
    return(c(lambda-1,lambda))
  } else {
    return(x)
  }
}


# Negative binomial distribution
#' @inheritParams stats::dnbinom
#' 
#' @export
#' @rdname distrMode
#' 
nbinomMode <-
function(size,
         prob,
         mu)
{
  if (prob > 1 || prob < 0) return(NaN)
  if (!missing(mu)) {
    prob <- size/(size+mu)    
  }
  if (size <= 1) {
    0
  } else {
    floor((size-1)*(1-prob)/prob)
  }
}


# Poisson distribution
#' @inheritParams stats::dpois
#' 
#' @export
#' @rdname distrMode
#' 
#' @examples 
#' ## Poisson distribution
#' poisMode(lambda = 6)
#' poisMode(lambda = 6.1)
#' mlv("poisson", lambda = 6.1)
#' 
poisMode <-
function(lambda)
{
  if (lambda < 0) return(NaN)
  if (lambda == 0) return(0)
  x <- floor(lambda)
  if (lambda == x) {
    c(lambda-1,lambda)
  } else {
    x
  }
}
paulponcet/modeest documentation built on Nov. 19, 2019, 9:16 p.m.