# Algorithm for generating hyperexponential distribution

I need to generate a hyperexponential distribution for my project. I have already implemented a poisson generating algorithm given by Donald Knuth, but I couldn't find an algorithm for generating a hyper exponential random variable.

I am provided with the mean and variance required of the distribution and I need an algorithm which can generate a random variable from this distribution when I execute it.

• Hi and welcome to the site, @Anirudh Vemula. This book contains an algorithm for generating RVs from a hyperexponential distribution (page 107). This site also seems to offer some guidance. Here is an algorithm implemented in C. Jun 17, 2013 at 10:21
• Thanks a lot! There is an alias method which is used in the algorithm(in book) to select lambda, do you have any idea what that is? Jun 17, 2013 at 10:34
• Unfortunately, I don't know what the alias method is. But a quick Google search reveiled several interesting documents (especially the third): first, second, third, fourth. Jun 17, 2013 at 10:39
• The alias method is just a method for generating discrete random variates; since the hyperexponential is a finite mixture of exponentials you need to choose which exponential component to generate first. The alias method is particularly convenient if the discrete variable has a finite range (as here). I describe the basic approach here along with some other methods. Nov 10, 2018 at 4:49
• @Anirudh just having a mean and variance is not sufficient to determine a hyperexponential. Even a two-component exponential has 3 parameters ($\lambda_1,\lambda_2, p_1$,) while specifying mean and variance will only fix two parameters. [Once you've chosen your $p$ and $\lambda$ vectors, it's very easy to generate. e.g. it's a single line in R: rexp(n,lambda[sample(p,n,TRUE,p)]) generates n values from the hyperexponential if you've already specified p, lambda and n. ] Nov 10, 2018 at 4:51

If you have a Hyperexponential-2 (H2), then with probability $$p$$ you sample from $$F_1$$ and with $$1-p$$ sample from $$F_2$$. Obviously $$p$$ must be on the interval $$[0,1]$$, $$F_1$$ ~ Exponential($$\lambda_1$$), and $$F_2$$ ~ Exponential($$\lambda_2$$).

You can mix a larger number of exponentials by extending this idea.

If you have a target mean and variance, select $$p,\lambda_1,\lambda_2$$ so that the resulting hyperexponential will have the target mean and variance.

To choose the parameters, solve the equations for the mean and variance. If $$X$$ ~ $$H2(p,\lambda_1,\lambda_2)$$ then $$\text{E}[X] = p/\lambda_1 + (1-p)/\lambda_2$$.
$$\text{Var}[X] = \text{E}[X^2]-\text{E}[X]^2 = 2p/\lambda_1 ^2 + 2(1-p)/\lambda_2^2 - \text{E}[X]^2$$

Set these equal to your targets and solve for $$p,\lambda_1,\lambda_2$$; it will be underdetermined so a number of hyperexponentials will have your target mean and variance. Just pick one solution.

MATLAB Code to generate from $$H2(p,\lambda_1,\lambda_2)$$:

function [ X ] = h2rnd( p,Rates,N )
%H2RAND Samples from specified Hyperexponential distrbution
%   [ X ] = h2rnd(p,Rates,N)
%   INPUTS
%        p: mixing probability
%    Rates: rates for the two exponential distributions (2 x 1 vector)
%        N: number of samples to generate
%  OUTPUTS
%        X ~ H2(p,Rates(1),Rates(2))
% %%%%%%%%%%%%% BEGIN ERROR CHECKING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if p < 0 | p > 1, error('p must be on interval [0,1]'), end
if min(Rates)<=0, error('Rates must be positive'), end
% %%%%%%%%%%%%% END ERROR CHECKING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X = zeros(N,1);
P = rand(N,1);
X(P<=p)= -log(1-rand(sum(P<=p),1))/Rates(1);
X(P>p)= -log(1-rand(sum(P>p),1))/Rates(2);
end


Update: Adding another exponential to have $$Y\sim H3(p_1,p_2,\lambda_1,\lambda_2,\lambda_3)$$ allows one to match more moments but increases the complexity as well.

Reference:
Wiki for Hyperexponential Distribution

Composition method is used in @SecretAgentMan's answer.

Instead, inverse transformation is used in the following Python function.

import random as rd
import math
from scipy.optimize import fsolve

def sim_exp_hyper(pmf:list, expects:list, n:int, whi_seed:int=123) -> list:
"""Simulate realisation of hyper-exponentially distributed random
variables using inverse transformation method.

Keyword Arguments
=================
pmf: probability mass function
expects: expectation of each individual exponential distribution
n: number of realisations
whi_seed: seed
"""

if len(pmf) != len(expects):
raise ValueError("len(pmf) != len(expects).")
elif sum([i < 0 for i in pmf]) > 0:
raise ValueError(f"There are negative values in the probablity mass "
f"function {pmf}.")

def get_eq(u, x):
eq = 1 - sum(pmf[i] * math.exp(- x / expects[i]) for i in range(len(pmf))) - u
return eq

rd.seed(whi_seed)
us = [rd.random() for i in range(n)]
xs = [fsolve(lambda x: get_eq(u, x), 0.1)[0] for u in us]
return xs, us


This method is inefficient because a nonlinear equation is solved to generate one realisation. There is no way to express the inverse of CDF of hyper-exponential distributions. So we can only obtain $$x$$ by solving: $$F(x)=1-\sum_{i=1}^{m} p_{i} e^{-\lambda_{i} x}=\sum_{i=1}^{m} p_{i}\left(1-e^{-\lambda_{i} x}\right) = u$$ where $$u$$ is a simulated random number from uniform distribution over $$[0, 1]$$.

However, one random number (instead of two in composition method) is used for every realisation, which is essential in common random number method for variance reduction.