Skip to content

Instantly share code, notes, and snippets.

# marcovivero/quantileFunctions.R Created Jan 22, 2016

 ### 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)) }
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.