5
$\begingroup$

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.

|cite|improve this question
$\endgroup$
  • 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 -Reinstate Monica 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 -Reinstate Monica Nov 10 '18 at 4:51
4
$\begingroup$

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

|cite|improve this answer
$\endgroup$

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.