Algorithm for generating hyperexponential distribution

Question

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.

Answer 1

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

Answer 2

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.