Implementation of Exponential Distribution
Hi, I've been looking into statistical distributions, and while comparing Wikipedia's mathematical formulation and Boost's implementation of the Exponential Distribution, I have found two points of interest. Before continuing, I just need to add that I'm not a statistician or a mathematician, so any errors either mathematical or in interpretation below are my own. The Wikipedia reference: http://en.wikipedia.org/wiki/Exponential_distribution The Boost implementation (line 54 of exponential_distribution.hpp) -result_type(1) / _lambda * log(result_type(1)-eng()); Where the eng() call generates a uniform random number. The first, less critical point: If eng() is uniform, then so is (1 - eng()), which allows for the simplification of the statement passed to log() from: log(result_type(1)-eng()) to log(eng()) The second, more critical point: The implementation does not seem to be the correct inverse of the original function. Boost's implementation (with above modification) is: -result_type(1) / _lambda * log(eng()); But should actually be: -log(eng()) / _lamda; Could someone perhaps clarify the above? Again, apologies if my interpretations were wrong. Kind regards, Nelis Franken
The second, more critical point: The implementation does not seem to be the correct inverse of the original function. Boost's implementation (with above modification) is: -result_type(1) / _lambda * log(eng());
But should actually be: -log(eng()) / _lamda;
The expressions are the same because of operator associativiy. (i.e. x/y*z = (x/y)*z ) I believe you are right about the first point though. cheers Arnaldur
Nelis Franken wrote:
Hi,
I've been looking into statistical distributions, and while comparing Wikipedia's mathematical formulation and Boost's implementation of the Exponential Distribution, I have found two points of interest. Before continuing, I just need to add that I'm not a statistician or a mathematician, so any errors either mathematical or in interpretation below are my own.
The Wikipedia reference: http://en.wikipedia.org/wiki/Exponential_distribution
The Boost implementation (line 54 of exponential_distribution.hpp) -result_type(1) / _lambda * log(result_type(1)-eng());
Where the eng() call generates a uniform random number.
The first, less critical point: If eng() is uniform, then so is (1 - eng()), which allows for the simplification of the statement passed to log() from: log(result_type(1)-eng()) to log(eng())
The second, more critical point: The implementation does not seem to be the correct inverse of the original function. Boost's implementation (with above modification) is: -result_type(1) / _lambda * log(eng());
But should actually be: -log(eng()) / _lamda;
Could someone perhaps clarify the above? Again, apologies if my interpretations were wrong.
No you are correct on both counts I believe, I'm cc'ing the random lib's author: Jens, any comments? John.
Thanks for the forward. John Maddock wrote:
Nelis Franken wrote:
The Wikipedia reference: http://en.wikipedia.org/wiki/Exponential_distribution
The Boost implementation (line 54 of exponential_distribution.hpp) -result_type(1) / _lambda * log(result_type(1)-eng());
This is -1/lambda * log(1-eng()) This matches what wikipedia shows.
The first, less critical point: If eng() is uniform, then so is (1 - eng()), which allows for the simplification of the statement passed to log() from: log(result_type(1)-eng()) to log(eng())
Not quite. Our engines return values in the range [0,1) (lower bound inclusive, upper bound exclusive). log(0) is undefined, thus we avoid that by using (1-eng()). Wikipedia states "Moreover, if U is uniform on (0;1), then so is 1 − U. " Note the interval that is open on both ends in the precondition. I believe the boost implementation is correct here. Thanks, Jens
On Mon, 23 Oct 2006, Nelis Franken wrote: [...]
The Boost implementation (line 54 of exponential_distribution.hpp) -result_type(1) / _lambda * log(result_type(1)-eng());
Where the eng() call generates a uniform random number.
The first, less critical point: If eng() is uniform, then so is (1 - eng()), which allows for the simplification of the statement passed to log() from: log(result_type(1)-eng()) to log(eng())
If I remember correctly, eng() returns a value in the interval [0,1), which poses a problem for log at 0. But 1-eng() returns a value in the interval (0,1], which is ok for log. -- François Duranleau LIGUM, Université de Montréal
participants (5)
-
Arnaldur Gylfason -
François Duranleau -
Jens Maurer -
John Maddock -
Nelis Franken