Numerical computation of quantiles

Recently I had to define a R function that would compute the {q}-th quantile of a continuous random variable based on an user-defined density function. Since the main objective is to have a general function that computes the quantiles for any user-defined density function it needs be done numerically.

Problem statement

Suppose we are interested in numerically computing the {q}-th quantile of a continuous random variable {X}. Assume for now that {X} is defined on {\Re = (-\infty, \infty)}. So, the problem is to find {x^*} such that {x^* = \{x : F_X(x) = q\}}, where {F_X(x)} is the cumulative distribution function of {X},

\displaystyle F_X(x) = \int _{-\infty} ^ {x} f_X(x) dx

and {f_X(x)} is its density function.


This problem can be cast into an optimization problem,

\displaystyle x^* = \underset{x}{\text{arg min}} (F_X(x) - q)^2

That said, all we need is a numerical optimization routine to find {x} that minimize {(F_X(x) - q)^2} and possibly a numerical integration routine to compute {F_X(x)} in case it is not available in closed form.

Below is one possible implementation of this strategy in R, using the integrate function for numerical integration and nlminb for the numerical optimization.

CDF <- function(x, dist){
objective <- function(x, quantile, dist){
  (CDF(x, dist) - quantile)^2
find_quantile <- function(dist, quantile){
  result = nlminb(start=0, objective=objective,
                  quantile = quantile,
                  dist = dist)$par

and we can test if everything is working as it should using the following
simple task of computing the 95th-quantile of the standard Gaussian distribution.

std_gaussian <- function(x){
  dnorm(x, mean = 0, sd = 1)
find_quantile(dist = std_gaussian,
              quantile = 0.95)
[1] 1.644854
[1] 1.644854

Random variables with bounded domain

If {X} has a bounded domain, the optimization problem above becomes more tricky and unstable. To facilitate the numerical computations we can transform {X} using a monotonic and well behaved function {g}, so that {Y = g(X)} is defined on the real line. Then we apply the numerical solution given above to compute the {q}-th quantile of {Y}, say {y^*}, and then transform the quantile back to the original scale {x^* = g^{-1}(y^*)}. Finally, {x^*} is the {q}-th quantile of {X} that we were looking for.

For example, if {X = [0, \infty)} we could use {Y = \log(X)}.

Fast, simple and useful numerical integration methods

A summary of quick and useful numerical methods to compute a definite integral of the form:

\displaystyle I = \int_{a}^{b} f(x) dx

When implementing some tasks we often use unnecessary complex integration methods just because they are available in most scientific computing languages. They are usually more accurate but at the expense of speed. When computational time is important it is worth to know these faster and easy to implement integration methods.

Trapezoidal rule of integration

The trapezoidal rule of integration approximates the function {f(x)} to be integrated by a first order polynomial {f_1(x) = a_0 + a_1 x}. If we chose {(a, f(a))} and {(b, f(b))} as the two points to approximate {f(x)} by {f_1(x)} we can solve for {a_0} and {a_1}, leading to

\displaystyle \begin{array}{rcl} I & \approx & \int_{a}^{b} (a_0 + a_1x)dx\\ & = & a_0(b-a) + a_1 \bigg(\frac{b^2 - a^2}{2}\bigg)\\ & = & \frac{(b-a)}{2}(f(a) + f(b)) \end{array}

To be more accurate, we usually divide the interval {[a, b]} into {n} sub-intervals, such that the width of each sub-interval is {h = (b - a)/n}. In this case, the trapezoidal rule of integration is given by

\displaystyle I \approx \frac{(b-a)}{2n}\bigg[f(a) + 2\bigg(\sum_{i=1}^{n-1} f(a+ih)\bigg) + f(b)\bigg]

The total error on the multi-segment trapezoidal rule is given by

\displaystyle Error = -\frac{(b-a)^3}{12n^3} \sum_{i = 1}^{n} f''(\zeta_i),

where {\zeta_i} is some point in the interval {[a + (i-1)h, a + ih]}.

Simpson {1/3} rule of integration

The Simpson {1/3} rule of integration approximates the function {f(x)} to be integrated by a second order polynomial {f_2(x) = a_0 + a_1 x + a_2 x^2}. If we chose {(a, f(a))}, {((a+b)/2, f((a+b)/2))} and {(b, f(b))} as the three points to approximate {f(x)} by {f_2(x)} we can solve for {a_0}, {a_1} and {a_2}, leading to

\displaystyle \begin{array}{rcl} I & \approx & \int_{a}^{b} (a_0 + a_1x + a_2 x^2)dx\\ & = & \frac{h}{3}\bigg[f(a) + 4f\bigg(\frac{a+b}{2}\bigg) + f(b)\bigg] \end{array}

Again, if we divide the interval {[a, b]} into {n} equidistant sub-intervals each of length {h}, such that {a = x_0 < x_1 < x_2 < ... < x_n = b}, and apply the Simpson {1/3} rule of integration in each sub-interval we get

\displaystyle I \approx \frac{(b-a)}{3n}\bigg[f(x_0) + 4\sum_{i=1,\ i=\text{odd}}^{n-1} f(x_i) + 2\sum_{i=2,\ i=\text{even}}^{n-2} f(x_i) + f(x_n)\bigg]

The total error on the multi-segment Simpson {1/3} rule is given by

\displaystyle Error = -\frac{(b-a)^5}{90n^5} \sum_{i = 1}^{n/2} f^{(4)}(\zeta_i),

where {\zeta_i} is some point in the interval {[x_{i-1}, x_{i+1}]}.

Gaussian quadrature rule of integration

A quadrature rule of integration [2] approximates the integral {I} by a weighted sum of function values at specified points within the domain of integration,

\displaystyle I \approx \sum_{i=1}^{n} w_i f(x_i).

The Gaussian quadrature rule is constructed to yield an exact result for polynomials of degree {2n - 1} or less by a suitable choice of points {x_i} and weights {w_i}, {i=1,...,n}.

By convention the domain of integration for such a rule is taken as {[-1, 1]},

\displaystyle \int_{-1}^{1} f(x) dx \approx w_i f(x_i),

which means that applying the Gaussian quadrature to solve an integral over {[a, b]} takes the following form:

\displaystyle I \approx \frac{b-a}{2} \sum_{i=1}^{n} w_i f\bigg(\frac{b-a}{2}x_i + \frac{a+b}{2}\bigg)

It can be shown that the evaluation points are just the roots of a polynomial belonging to a class of orthogonal polynomials. The approximation above will give accurate results if the function {f(x)} is well approximated by a polynomial of degree {2n-1} or less. In that case the associated polynomials are known as Legendre polynomials and the method is known as Gauss-Legendre quadrature.

Here we can find the Gauss-Legendre’s quadrature weights {w_i} and Abscissae values {x_i}, {i=1,...,n} for different values of {n}.

The total error of the Gauss-Legendre quadrature rule is given by

\displaystyle \frac{(b-a)^{2n+1}(n!)^4}{(2n+1)[(2n)!]^3} f^{(2n)}(\zeta), \quad a<\zeta<b

Alternative Gaussian quadrature rules

If the integrated function can be written as {f(x) = w(x)g(x)}, where {g(x)} is approximately polynomial, and {w(x)} is known, then there are alternative weights {w_i'} such that

\displaystyle \int_{-1}^{1} f(x) dx = \int_{-1}^{1}w(x)g(x)dx \approx \sum_{i=1}^{n} w_i'g(x_i)

For example, the Gauss-Hermite quadrature rule is used when {w(x) = e^{-x^2}}.


[1] Integration section of the Holistic Numerical Methods website.
[2] Wikipedia page on Gaussian quadrature.
[3] Gauss-Legendre’s quadrature weights and Abscissae values for different values of n.