Re: How do I get random numbers that follow standard distribution?

The Perl Cookbook provides this code:

sub gaussian_rand {
    my ($u1, $u2);  # uniformly distributed random numbers
    my $w;          # variance, then a weight
    my ($g1, $g2);  # gaussian-distributed numbers

    do {
        $u1 = 2 * rand() - 1;
        $u2 = 2 * rand() - 1;
        $w = $u1*$u1 + $u2*$u2;
    } while ( $w >= 1 );

    $w = sqrt( (-2 * log($w))  / $w );
    $g2 = $u1 * $w;
    $g1 = $u2 * $w;
    # return both if wanted, else just one
    return wantarray ? ($g1, $g2) : $g1;
}
[download]

and notes:

The gaussian_rand function [generates] two numbers with a mean of 0 and a standard deviation of 1 (i.e., a Gaussian distribution). To generate numbers with a different mean and standard deviation, multiply the output of gaussian_rand by the new standard deviation, and then add the new mean

Comment on Re: How do I get random numbers that follow standard distribution? Download Code

Replies are listed 'Best First'.
Re: Answer: How do I get random numbers that follow standard distribution? by Roy Johnson (Monsignor) on Feb 17, 2005 at 16:04 UTC
This page claims that the while condition should be `} while ($w>=1 or $w==0);` Caution: Contents may have been coded under pressure.	[reply] [d/l]
Re^2: Answer: How do I get random numbers that follow standard distribution? by hv (Prior) on Feb 18, 2005 at 11:56 UTC
`$w == 0` will occur only when `$u1 == $u2 == 1/2`, which will happen rarely enough that I can't see it skewing the results noticeably - even if you only have 32-bit random numbers, you'll hit this case only 1.27 in 2^64 times. (That 1.27 is based on my rusty attempts to calculate how many times we loop until we get `$w < 1`. I think the probability of a success is `arcsin(1)/2`, which is about 0.78.) Hugo	[reply] [d/l] [select]
Re: Answer: How do I get random numbers that follow standard distribution? by polettix (Vicar) on Mar 24, 2005 at 17:38 UTC
I do support the previous answer: the code should include the test for `$w` being zero, otherwise you risk a division by zero in the very next line. The proposed solution roots in Numerical Recipes in C, chapter 7, section 2; probably it' better to update the Perl Cookbook and the answer in Q&A as well. Flavio Don't fool yourself.	[reply] [d/l]