# 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. – COOLSerdash Jun 17 '13 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? – Anirudh Vemula Jun 17 '13 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. – COOLSerdash Jun 17 '13 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. – Glen_b -Reinstate Monica Nov 10 '18 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. ] – Glen_b -Reinstate Monica Nov 10 '18 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