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