The exponential mixture is a semi-parametric distribution used in Loss Modeling to fit collections of similiar losses. The exponential mixture has a decreasing hazard rate, increasing mean residual life and has a thicker tail than the usual single parameter exponential. The exponential mixture is considered semi-parametric because no parametric assumption is made about the form of the mixing distribution.

The R library Renext exposes the MixExp2 function, which makes available the standard d, p, q and r-prefixed statistical functions to determine the density, distribution, quantile or generate random draws from 2-component exponential mixtures, but for exponential mixtures of more than two terms there’s little available.

In this post, we’ll demonstrate a method for calculating quantiles from exponential mixture models with any number of terms using inverse transform sampling. We’ll leverage the Newton-Raphson method to obtain quantiles, assuming the mean and weight for each component have already been determined. In a future post, we’ll demonstrate how to use maximum likelihood and MCMC to arrive at the optimal weight distribution across terms.

Background

The distribution function of an exponential mixture model can be represented as

$$ F(x) = 1 - w_{1}e^{-\lambda_{1}x} - w_{2}e^{-\lambda_{2}x} - \cdots - w_{n}e^{-\lambda_{n}x}, $$

where \(\sum_{i}^{n}w_{i} = 1\) and \(\lambda_{i}\) is the inverse of the expected value for the \(i^{th}\) component in the mixture. In order to use Newton-Raphson to find roots for an exponential mixture model with an arbitrary number of components, we need to arrange the CDF such that all terms on one side are set identically to zero. Rearranging yields

$$ w_{1}e^{-\lambda_{1}x} + w_{2}e^{-\lambda_{2}x} + \cdots + w_{n}e^{-\lambda_{n}x} - (1 - F(x)) = 0 $$

The Newton-Raphson update expression can be represented as

$$ x_{n+1} = x_{n} - \frac{f(x)}{f^{\prime}(x)}, $$

with successive iterations providing better approximations of the root until a sufficiently accurate value is obtained.

For our purposes, the rearranged CDF expression will serve as \(f(x)\). \(f^{\prime}(x)\) will be the derivative of the rearranged CDF, specifically:

$$ f^{\prime}(x) = -\lambda_{1}w_{1}e^{-\lambda_{1}x} -\lambda_{2}w_{2}e^{-\lambda_{2}x} - \cdots -\lambda_{n}w_{n}e^{-\lambda_{n}x}. $$

We can commence iteration at \(x_0 = 0\) in all cases, and don’t need to concern ourselves with searching for a optimal initialization point since the mixture is monotonic.

Implementation

What follows is an implementation of the Newton-Raphson algorithm written in R. It takes as arguments a starting point, a function f and its derivative f_prime. epsilon and maxiters represent the threshold below which a root is assumed to have been located and and the maximum number of iterations to carry out as a safeguard against divergence:

newtonRaphson = function(
    x_init=0, f, f_prime, epsilon=.000001, 
    maxiters=1000, returniters=FALSE) {

    root = NULL
    x_i = x_init
    tolerance = 1
    itervals = vector()
    itercntr = 0    # ensure `maxiters` is not exceeded

    while (TRUE) {
        itercntr = itercntr + 1
        itervals[[itercntr]] = x_i
        if ((itercntr>=maxiters) || (tolerance<=epsilon))  {
            root = itervals[[length(itervals)]]
            break
        } else {
            x_n = (x_i - (f(x_i) / f_prime(x_i)))
            tolerance = abs(x_n - x_i)
            x_i = x_n 
        }
    }

    resList = list(iterations=itercntr, root=root, approx_seq=itervals)
    if (returniters==FALSE) resList$approx_seq = NULL
    return(resList)
}

The function returns a list of three elements:

  • iterations is the number of Newton-Raphson approximations.
  • root is the value at which the minimum of the function is obtained.
  • approx_seq is a vector of length equal to iterations of successive root approximations.

Numerical Example

Assume we’re interested in determining quantiles for an exponential mixture model with weights = (.60, .30, .10) and means = (10, 50, 100). The corresponding distribution function is

$$ F(x) = 1 - .60e^{-x/10} - .30e^{-x/50} - .10e^{-x/100}, $$

where the hazard rate for each term is the inverse of the term’s mean.

As a sanity check, let’s evaluate \(F(x)\) at a few points to determine the CDF, which will serve as a stand-in for the uniform random variates that will be sampled when the algorithm is used in practice. Let’s find the associated percentiles for inputs 10, 25, 75 and 200:

mixtureCDF = function(x) {

    return(1 - .60 * exp(-x / 10) - .30 * exp(-x / 50) - .10 * exp(-x / 100))

}

mixtureCDF(10)    #  0.4431693676
mixtureCDF(25)    #  0.6909097246
mixtureCDF(75)    #  0.8854924461
mixtureCDF(200)   #  0.9809717788

These percentiles represent the uniform random numbers that will be passed to the simulation framework in practice. The goal is upon being passed a random number such as 0.6909097246 (corresponding to the mixture model evaluated at 25), the function will return 25.

We need one additional piece of functionality to complete the framework: A function that wraps newtonRaphson and takes for arguments a vector of means, a vector of weights and returns another function which can then be passed a uniform random variate, which ultimately returns the corresponding quantile. This logic is encapsulated in invertCDF:

invertCDF = function(weights, means) {

    rates = 1.0 / means

    function(urand) {
        # Initialize numerator and denominator functional
        # components for Newton-Raphson iteration.
        f = function(x) {
            func = 0
            for (i in 1:length(rates)) {
                iterrate = rates[i]
                iterweight = weights[i]
                iter_func =  iterweight * exp(-iterrate * x)
                func = func + iter_func
            }
            return(func - (1 - urand))
        }

        f_prime = function(x) {
            func_prime = 0
            for (i in 1:length(rates)) {
                iterrate = rates[i]
                iterweight = weights[i]
                iter_func_prime =  -((iterweight * iterrate) * exp(-iterrate * x))
                func_prime = func_prime + iter_func_prime
            }
            return(func_prime)
        }

        # call newtonRaphson to back-out quantile from mixed exponential
        # distribution corresponding to CDF value of urand.
        result = newtonRaphson(x_init=0, f=f, f_prime=f_prime)
        return(result$root)
    }
}

Finally, we demonstrate using invertCDF to find quantiles from the sample exponential mixture CDF \(F(x) = 1 - .60e^{-x/10} - .30e^{-x/50} - .10e^{-x/100}\). Let’s assume our uniform random samples are 0.4431693676, 0.6909097246, 0.8854924461 and 0.9809717788, which correspond to values 10, 25, 75 and 200. Put another way, if we pass 0.4431693676, 0.6909097246, 0.8854924461 and 0.9809717788 into our simulation framework, we’d expect the values 10, 25, 75 and 200 to be returned. This serves as verification that the framework is functioning as expected:

weights = c(.60, .30, .10)
means = c(10, 50, 100)
random_draws = c(0.4431693676, 0.6909097246, 0.8854924461, 0.9809717788)  

# Call invertCDF, which returns a function.
inverter = invertCDF(weights=weights, means=means)

# Pass random_draws to inverter, and obtain estimated percentiles.
result = sapply(random_draws, inverter, simplify=TRUE)

print(result)
# [1]  10.00000000  25.00000000  75.00000002 200.00000011

Below are plots of the exponential mixture model, along with the shaded area under the curve that corresponds to the quantiles found above:

em44

em69

em86

em98

The plots were created with ggplot2 as follows:

library(ggplot2)

f = function(x) {

    return(.06 * exp(-x / 10.) + .006 * exp(-x / 50.)+.001 * exp(-x / 100.))

}

x = seq(from=0,to=100,by=.005)
y = sapply(X=x,FUN=f,simplify=TRUE,USE.NAMES=FALSE)
df = data.frame(X=x,Y=y)

ggplot(data=df, mapping=aes(x=x, y=y)) +
    geom_line() +
    geom_area(mapping=aes(x=ifelse(x>0 & x<10, x, 0)), fill="#9898fb", alpha=1.) +
    xlim(0, 100) +
    ylim(0, .1) + 
    ylab("Exponential Mixture Density") +
    ggtitle("AUC for Exponential Mixture: ~44th Percentile") +
    geom_vline(xintercept=0) +
    geom_hline(yintercept=0)