quantileFunctions.R

### QUANTILE FUNCTIONS ###

# Load required packages.
library(MCMCpack)
library(rootSolve)

# Define desired quantile.
p = 0.9

### EMPIRICAL QUANTILE ESTIMATION ###

# Helper function for computing empirical quantile estimates.
# Input: a numeric vector x.
# Output: a numeric quantile estimate.
empiricalQuantileEstimate = function(x) {
  quantile(x, c(p))[1]
}

### INVERSE-GAMMA QUANTILE ESTIMATION ###

# Helper function for determining whether critical point is a local maximum.
maximumCheck = function(n, alpha, beta) {
  D = (n^2 / beta^2) * (alpha * trigamma(alpha) - 1)
  secDer = -n * trigamma(alpha)
  
  return(D > 0 & secDer < 0)
}

# Helper function for computing MLE estimates for an inverse-gamma 
# distribution.
igMLEEstimate = function(x) {
  # Define function for which we wish to find a root for alpha.
  igConstant = -(log(mean(1 / x)) + mean(log(x)))
  igMLEFunc = Vectorize(function(x) {
    igConstant + log(x) - digamma(x)
  })
  # Solve equation for roots.
  alphaRoots = uniroot.all(igMLEFunc,c(0.000001, 1000))
  betaRoots = alphaRoots / mean(1 / x)
  mleMatrix = matrix(c(alphaRoots, betaRoots), byrow = T)
  maximumIndices = which(
    apply(mleMatrix, 2, function(z) {maximumCheck(length(x), z[1], z[2])}))
  
  return(mleMatrix[, maximumIndices])
}

# Helper function for creating inverse-gamma density function.
# Input: parameters alpha and beta.
# Output: resulting density function.
getInvGammaDensity = function(alpha, beta) {
  d = Vectorize(function(x) {
    dinvgamma(x, shape = alpha, scale = beta)
  })
  return(d)
}

# Helper function for computing quantiles of arbitrary density functions.
# Input: Density function func.
# Output: Quantile equation function to be solved.
quantileFunc = function(p, func) {
  output = Vectorize(function(x) {
    integrate(func, 0, x)$value - p
  })
  
  return(output)
}

# Helper function for computing inverse-gamma quantile estimate of an input 
# numeric vector x. 
igQuantileEstimate = function(x) {
  mleEstimates = igMLEEstimate(x)
  densityFunc = getInvGammaDensity(mleEstimates[1], mleEstimates[2])
  
  uniroot.all(quantileFunc(p, densityFunc), c(0.000001, 5000))
}