Skip to content

Instantly share code, notes, and snippets.

@schnellp
Created April 28, 2015 13:09
Show Gist options
  • Save schnellp/d774c127d09d4b5b691e to your computer and use it in GitHub Desktop.
Save schnellp/d774c127d09d4b5b691e to your computer and use it in GitHub Desktop.
An R function for finding parameters of normal, beta, and gamma distributions that match a given pair of quantiles
params.from.quantiles <- function(q, p, family="normal", start=c(0.5, 0.5)) {
# Finds parameter values of distributions that closely match
# the specified quantiles.
#
# Args:
# q: A vector of quantiles (length 2, strictly increasing)
# p: A vector of percentiles (length 2, strictly increasing)
# family: One of "normal", "beta", or "gamma" specifying
# which family of distributions to use
# start: starting values for parameter search
#
# Returns:
# A list containing a named vector of parameter values (par),
# a named vector containing the fitted quantiles (cdf),
# and the family used (family)
# basic input validation
stopifnot(length(p) == 2,
length(q) == 2,
p[1] < p[2],
q[1] < q[2],
0 < p[1],
p[2] < 1)
# select distribution function
g <- switch(family,
normal=pnorm,
beta=pbeta,
gamma=pgamma)
# a wrapper of the cdf for optim
f <-function(theta) {
# if parameter values are invalid, return NA
tryCatch({
sum((g(q, theta[1], theta[2]) - p)^2)
}, warning=function(w) {
NA
})
}
# find best parameter values
sol <- optim(start, f)
# prepare output
out <- list()
out$par <- sol$par
names(out$par) <- switch(family,
normal=c("mean", "sd"),
beta=c("a", "b"),
gamma=c("shape", "rate"))
out$cdf <- g(q, out$par[1], out$par[2])
names(out$cdf) <- q
out$family <- family
out
}
# example usage
# params.from.quantiles(q=c(0.01, 0.05), p=c(0.5, 0.99), family="beta")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment