Create a gist now

Instantly share code, notes, and snippets.

### 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))
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment