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.

  • 1
    $\begingroup$ 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. $\endgroup$ – COOLSerdash Jun 17 '13 at 10:21
  • $\begingroup$ 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? $\endgroup$ – Anirudh Vemula Jun 17 '13 at 10:34
  • $\begingroup$ 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. $\endgroup$ – COOLSerdash Jun 17 '13 at 10:39
  • $\begingroup$ 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. $\endgroup$ – Glen_b Nov 10 '18 at 4:49
  • 2
    $\begingroup$ @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. ] $\endgroup$ – Glen_b 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)
%        p: mixing probability
%    Rates: rates for the two exponential distributions (2 x 1 vector)
%        N: number of samples to generate
%        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);

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.

Wiki for Hyperexponential Distribution

| cite | improve this answer | |

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

    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.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.