An R function for finding parameters of normal, beta, and gamma distributions that match a given pair of quantiles

params-from-quantiles.R

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")