Beefy Boxes and Bandwidth Generously Provided by pair Networks
Don't ask to ask, just ask
 
PerlMonks  

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

by Roy Johnson (Monsignor)
on Feb 17, 2005 at 15:45 UTC ( #431978=categorized answer: print w/replies, xml ) Need Help??

Q&A > math > How do I get random numbers that follow standard distribution? - Answer contributed by Roy Johnson

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; }
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 Answer: 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.

      $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

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.
Log In?
Username:
Password:

What's my password?
Create A New User
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others scrutinizing the Monastery: (3)
As of 2020-10-22 21:27 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?
    My favourite web site is:












    Results (230 votes). Check out past polls.

    Notices?