# How to sample from an exponential distribution using rejection sampling in PHP

I wrote PHP-code that gets me samples of a (truncated) exponential distribution between 0 and 1 with mean 1 ($X\sim \mathrm{Exp}(1)$).

I'm trying to use acceptance-rejection method: I don't know if I get it right, but basically I use one sample of a $\mathrm{Uniform}(0,1)$ distribution to get the x-axis and another to get the y-axis and then check if the y-axis sample is below (valid sample) or above (reject sample) the exponential distribution function curve.

How do I get it to give me samples of a generic $X\sim \mathrm{Exp}(\lambda)$?

$num_samples=1000;$samples=array();$sample=0;$count=0;$counts=0;$time=microtime(true);
for($i=0;$i<$num_samples;$i++){
while (exp(-($sample = mt_rand(0,PHP_INT_MAX)*(1/PHP_INT_MAX))) < mt_rand(0,PHP_INT_MAX)*(1/PHP_INT_MAX)){$count++;
}
$samples[] =$sample;
$counts+=$count;
$count=0; } echo "samples: ".$num_samples." time: ".(microtime(true)-$time)." efficiency: ".$num_samples/$counts."\n";  • An impressively efficient and general answer to this question appears at stats.stackexchange.com/a/13767. – whuber Aug 25 '13 at 20:06 • @whuber Wow, that's much better, I got the best method for an exponential, but this could be very useful sometime. Thank you! Aug 26 '13 at 8:25 ## 2 Answers Your exponential random variable$X$truncated between$0$and$\tau>0$has distribution function $$F_X(x) = \begin{cases} 0 & \textrm{if} & x\leq 0 \,;\\ \frac{1-e^{-\lambda x}}{1-e^{-\lambda\tau}} & \textrm{if} & 0<x<\tau \, ;\\ 1 & \textrm{if} & x \geq \tau \, .\\ \end{cases}$$ Using this result, we can generate realizations of$X$with the following PHP code. <?php$lambda = 1;
$tau = 10;$u = mt_rand(0, PHP_INT_MAX) / PHP_INT_MAX;
$x = - log(1 - (1 - exp(-$lambda * $tau)) *$u) / $lambda; echo$x;
?> samples: 1000000 time: 3.1613190174103 efficiency: 1

• I agree with this advice, but it looks like neither your code nor the formula appear to generate a truncated exponential.
– whuber
Aug 25 '13 at 19:16
• Thanks, Bill. It seems that he added the truncation to the question later, or I'm really tired...
– Zen
Aug 25 '13 at 19:42
• I did. Also my code had some things wrong (sorry it was a first try, should have think it more before posting) I posted an answer correcting them and solving the problem I had. I'll check out your code too. Thanks! Aug 25 '13 at 19:46
• Which code is faster?
– Zen
Aug 25 '13 at 19:48
• This method is much faster. Like 2.5 times faster than mine truncating at 10 the normalized exponential. I think some improvements on the implementation can be made, though. Next time I'm at my parents I will dust off my class notes from this subject on 2008. This is the first time I apply something from this course since then, but I remember there was something about it there. Thanks a lot :) Aug 26 '13 at 7:41

Ok, I just had a moment to take another look at it and following @Glen_b directions I could figure it out. It's terribly inefficient, but in case anybody is interested, here is the code:

$num_samples=1000;$avg = 2; /* ($lambda = 1/$avg) */
$trunc = 10; /* ($tau = $trunc*$avg) */
$start_time = microtime(true); for($i=0;$i<$num_samples;$i++){ do {$sample_exp_lambda_1_tau_10 = mt_rand(0,PHP_INT_MAX)*($trunc/PHP_INT_MAX);$p_exp_lambda_1_tau_10 = exp(-$sample_exp_lambda_1_tau_10);$y_axis = mt_rand(0,PHP_INT_MAX)*(1/PHP_INT_MAX);
$count++; } while ($p_exp_lambda_1_tau_10 < $y_axis);$samples[] = $avg*$sample_exp_lambda_1_tau_10; /* ($x =$avg*$sample_exp_lambda_1_tau_10) */$total_count+=$count;$count=0;
}
echo "samples: ".$num_samples." time: ".(microtime(true)-$start_time)." efficiency: ".$num_samples/$total_count."\n";


And that's all. With this input I generated 1000 samples of a X~exp(2) truncated at 20.

Here is the output:

samples: 1000 time: 0.018975019454956 efficiency: 0.1012555690563


And here the histogram: (Close enough.)

• How could you possibly get values above 10 in your histogram when you're trying to truncate the distribution at 10?
– whuber
Aug 25 '13 at 19:44
• wow, I have to check that out. Aug 25 '13 at 19:48
• Yeah, sorry. I am truncating at $avg*$trunc. It's 20. (I am truncating at 10 the normalized function, then multiplying by lambda to get exp(lambda) as suggested) Aug 25 '13 at 19:49