Dear Booster, I have a problem in understanding how the boost::random::discrete_distribution works. As far as I know, given the following events: x1, x2 and x2 and weights: 1, 4, and 5, * the associated PDF of the discrete distribution would be: \Pr{X=x1}=1/10 \Pr{X=x2}=4/10 \Pr{X=x3}=5/10 * and the CDF would be: \Pr{X<=x1}=1/10 \Pr{X<=x2}=5/10 \Pr{X<=x3}=1 So, if my random generator generates the value 0.814724, the discrete random variate would return the event x3 (if I'm not wrong). Instead, it seems that boost::random::discrete_distribution does not work in this way. To test it I have implemented a small program (see attachment). Moreover, I have also implemented a naive version of discrete distribution (called my::discrete_distribution), which simply draws a uniform [0,1] random number and select the first event whose CDF is greater than or equal to the drawn value. I have also added some debug prints in other boost classes (just to get the value of the generated random numbers). I run it in Linux with Boost devel version + GCC 4.6.3 (with flags: -Wall -ansi -pedantic) Here's the output: Testing Boost... [boost::random::mersenne_twister] 3499211612 [boost::random::mersenne_twister] 581869302 [boost::random::discrete_distribution] test=0.814724 - result=0 Random Event: 1 [boost::random::mersenne_twister] 3890346734 [boost::random::mersenne_twister] 3586334585 [boost::random::discrete_distribution] test=0.905792 - result=2 Random Event: 2 [boost::random::mersenne_twister] 545404204 [boost::random::mersenne_twister] 4161255391 [boost::random::discrete_distribution] test=0.126987 - result=2 Random Event: 2 Testing My own version... [boost::random::mersenne_twister] 3499211612 [my::discrete_distribution::rand] RNG: 0.814724 Random Event: 2 [boost::random::mersenne_twister] 581869302 [my::discrete_distribution::rand] RNG: 0.135477 Random Event: 1 [boost::random::mersenne_twister] 3890346734 [my::discrete_distribution::rand] RNG: 0.905792 Random Event: 2 As you can see, the sequence of Random Event is different. Looking at the output of my::distribution, it seems it is working OK. Can you help me in understanding what is going on? Thank you very much for the help!! Best, -- Marco
AMDG On 06/04/2012 01:53 PM, sguazt wrote:
Dear Booster,
I have a problem in understanding how the boost::random::discrete_distribution works.
As far as I know, given the following events: x1, x2 and x2 and weights: 1, 4, and 5, * the associated PDF of the discrete distribution would be: \Pr{X=x1}=1/10 \Pr{X=x2}=4/10 \Pr{X=x3}=5/10
* and the CDF would be: \Pr{X<=x1}=1/10 \Pr{X<=x2}=5/10 \Pr{X<=x3}=1
So, if my random generator generates the value 0.814724, the discrete random variate would return the event x3 (if I'm not wrong).
Wrong. The implementation does not use the inverse cdf. It uses Walker's alias algorithm which generated variates in constant time. The algorithm generates a table that looks like this: std::vector<std::pair<double, int> > table = { { 0.3, 1 }, { 0.5, 2 }, { 1.0, 0 } }; Then the generation algorithm is double x = uniform_01<>(engine); int result = uniform_int<>(0, 2)(engine); if (x < table[result].first) return result; else return table[result].second; If you work out the probabilities, you'll find that they add up correctly. In Christ, Steven Watanabe
On Tue, Jun 5, 2012 at 7:57 AM, Steven Watanabe <watanabesj@gmail.com> wrote:
AMDG
Thanks Steven! [cut]
The algorithm generates a table that looks like this:
std::vector<std::pair<double, int> > table = { { 0.3, 1 }, { 0.5, 2 }, { 1.0, 0 } };
Then the generation algorithm is
double x = uniform_01<>(engine); int result = uniform_int<>(0, 2)(engine); if (x < table[result].first) return result; else return table[result].second;
Sorry, I didn't have enough background knowledge on this topic. My error was to relate the U(0,1) x number to a CDF value. Anyway I performed a Monte-Carlo simulation and the results (the empirical PDF) match with the original PDF. Since I still don't understand how event weights are preserved and I don't find a way to get the Walker's paper http://dx.doi.org/10.1145/355744.355749 could you point me to some reference or URL so that I can get acquainted with Walker's method? Thank you so much! Best, -- Marco
On Tue, Jun 5, 2012 at 9:51 AM, sguazt <marco.guazzone@gmail.com> wrote:
On Tue, Jun 5, 2012 at 7:57 AM, Steven Watanabe <watanabesj@gmail.com> wrote:
AMDG
Thanks Steven!
[cut]
The algorithm generates a table that looks like this:
std::vector<std::pair<double, int> > table = { { 0.3, 1 }, { 0.5, 2 }, { 1.0, 0 } };
Then the generation algorithm is
double x = uniform_01<>(engine); int result = uniform_int<>(0, 2)(engine); if (x < table[result].first) return result; else return table[result].second;
Sorry, I didn't have enough background knowledge on this topic.
My error was to relate the U(0,1) x number to a CDF value. Anyway I performed a Monte-Carlo simulation and the results (the empirical PDF) match with the original PDF.
Since I still don't understand how event weights are preserved and I don't find a way to get the Walker's paper http://dx.doi.org/10.1145/355744.355749 could you point me to some reference or URL so that I can get acquainted with Walker's method?
Got it! I found Walker's method in Knuth "The Art of Computer Programming, Vol 2". Thank you -- Marco
participants (3)
-
Leo Cacciari -
sguazt -
Steven Watanabe