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.
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.
