Sample exactly from an exponential distribution. More...
|template<class Random >|
|RandomNumber< bits >||operator() (Random &r) const|
Sample exactly from an exponential distribution.
Sample x ≥ 0 from exp(−x). See:
The following code illustrates the basic method given by von Neumann:
This returns a result consuming exp(1)/(1 − exp(-1)) = 4.30 random numbers on average. (Von Neumann incorrectly states that the method takes (1 + exp(1))/(1 − exp(-1)) = 5.88 random numbers on average.) Because of the finite precision of Random::Fixed(), the code snippet above only approximates exp(−x). Instead, it returns x with probability h(1 − h)x/h for x = ih, h = 2−53, and integer i ≥ 0.
The above is precisely von Neumann's method. Abramowitz and Stegun consider a variant: sample uniform variants until the first is less than the sum of the rest. Forsythe converts the < ranking for the runs to ≤ which makes the analysis of the discrete case more difficult. He also drops the "trick" by which the integer part of the deviate is given by the number of rejections when finding the fractional part.
Von Neumann says of his method: "The machine has in effect computed a logarithm by performing only discriminations on the relative magnitude of numbers in (0,1). It is a sad fact of life, however, that under the particular conditions of the Eniac it was slightly quicker to use a truncated power series for log(1−T) than to carry out all the discriminations."
Here the code is modified to make it more efficient:
In addition, the method is extended to use infinite precision uniform deviates implemented by RandomNumber and returning exact results for the exponential distribution. This is possible because only comparisons are done in the method. The template parameter bits specifies the number of bits in the base used for RandomNumber (i.e., base = 2bits).
For example the following samples from an exponential distribution and prints various representations of the result.
Here's a possible result:
0.0111... = (0.4375,0.5) = 0.474126 10.000... = (2,2.125) = 2.05196 1.00... = (1,1.25) = 1.05766 0.010... = (0.25,0.375) = 0.318289 10.1... = (2.5,3) = 2.8732 0.0... = (0,0.5) = 0.30753 0.101... = (0.625,0.75) = 0.697654 0.00... = (0,0.25) = 0.0969214 0.0... = (0,0.5) = 0.194053 0.11... = (0.75,1) = 0.867946
First number is in binary with ... indicating an infinite sequence of random bits. Second number gives the corresponding interval. Third number is the result of filling in the missing bits and rounding exactly to the nearest representable double.
This class uses some mutable RandomNumber objects. So a single ExactExponential object cannot safely be used by multiple threads. In a multi-processing environment, each thread should use a thread-specific ExactExponential object. In addition, these should be invoked with thread-specific random generator objects.
|bits||the number of bits in each digit.|
|RandomNumber< bits > RandomLib::ExactExponential< bits >::operator()||(||Random &||r||)||const|
Return a random deviate with an exponential distribution, exp(−x) for x ≥ 0. Requires 7.232 bits per invocation for bits = 1. The average number of bits in the fraction = 1.743. The relative frequency of the results for the fractional part with bits = 1 is shown by the histogram
The base of each rectangle gives the range represented by the corresponding binary number and the area is proportional to its frequency. A PDF version of this figure is given here. This allows the figure to be magnified to show the rectangles for all binary numbers up to 9 bits. Note that this histogram was generated using an earlier version of ExactExponential (thru version 1.3) that implements von Neumann's original method. The histogram generated with the current version of ExactExponential is the same as this figure for u in [0, 1/2]. The histogram for u in [1/2, 1] is obtained by shifting and scaling the part for u in [0, 1/2] to fit under the exponential curve.
Another way of assessing the efficiency of the algorithm is thru the mean value of the balance = (number of random bits consumed) − (number of bits in the result). If we code the result in mixed Knuth and Yao's unary-binary notation (integer is given in unary, followed by "0" as a separator, followed by the fraction in binary), then the mean balance is 3.906. (Flajolet and Saheb analyzed the algorithm based on the original von Neumann method and showed that the balance is 5.680 in that case.)
For bits large, the mean number of random digits is exp(1/2)/(1 − exp(−1/2)) = 4.19 (versus 4.30 digits for the original method).
|Random||the type of the random generator.|
|[in,out]||r||a random generator.|