14.3 Normally Distributed Random Numbers

   

 
Java Number Cruncher: The Java Programmer's Guide to Numerical Computing
By Ronald  Mak

Table of Contents
Chapter  14.   Generating Random Numbers

14.3 Normally Distributed Random Numbers

Suppose we used a random number generator that produces uniformly distributed float values between 0 and 1 to produce n random values at a time, and then we averaged each group of values. How would the averages be distributed?

If we look at a large number of random values that are uniformly distributed between 0 and 1, their overall average should be 0.5. But if we average small groups of these values, say n = 12 at a time, then these group averages should "bunch up" around 0.5 some of the averages should be less than 0.5 and others greater, but most of them should be close to 0.5. In fact, the Central Limit Theorem of statistics states that these group averages will be normally distributed with a mean of 0.5. In other words, their distribution is the familiar bell curve with its peak at 0.5. So we've managed to derive our first algorithm for generating normally distributed random values.

Other algorithms for generating normally distributed random values are more obscure. We'll investigate two more that also start with uniformly distributed random values.

Method java.util.Random.nextGaussian() uses the polar algorithm. We begin this algorithm by generating two random values v 1 and v 2 that are uniformly distributed between -1 and +1. We then compute the value s = v 1 2 + v 2 2 . If s > 1, we reject the values of v 1 and v 2 and generate new ones. In other words, if we consider v 1 and v 2 to be the x and y coordinates of a point, we want only uniformly distributed random points that are inside the unit circle.

With two good values of v 1 and v 2 , we then compute

graphics/14equ02.gif


and

graphics/14equ03.gif


We won't prove it here the references at the end of this chapter contain the proofs but the values of x 1 and x 2 are normally distributed. This algorithm cranks them out two at a time.

Another algorithm is the ratio algorithm. Again, we start with two uniformly distributed random values, this time between 0 and 1, and which we'll call u and v. The value of u cannot be 0. We compute the value

graphics/14equ04.gif


Now we subject the value of x to up to three tests:

  • Quick acceptance test : If graphics/14inl02.gif then terminate the algorithm and deliver the value of x.

  • Quick rejection test : If graphics/14inl03.gif then reject this value of x and compute new values for u and v.

  • Final acceptance test : If x 2 -4 ln u then terminate the algorithm and deliver the value of x. Otherwise, reject this value of x and compute new values for u and v.

The values of x that this algorithm delivers are normally distributed. Again, see the references at the end of this chapter for proofs.

Listing 14-1a shows class RandomNormal in package numbercruncher.mathutils , which implements these three algorithms. Each of the three methods nextCentral(), nextPolar() , and nextRatio() returns the next random value from a normal distribution with a given mean and standard deviation.

Listing 14-1a Class RandomNormal , which implements three algorithms for generating normally distributed random values.
 package numbercruncher.mathutils; import java.util.Random; /**  * Utility class that generates normally-distributed  * random values using several algorithms.  */ public class RandomNormal {     /** mean */                 private float mean;     /** standard deviation */   private float stddev;     /** next random value from         the polar algorithm  */ private float   nextPolar;     /** true if the next polar         value is available  */  private boolean haveNextPolar = false;     /** generator of uniformly-distributed random values */     private static Random gen = new Random();     /**      * Set the mean and standard deviation.      * @param mean the mean      * @param stddev the standard deviation      */     public void setParameters(float mean, float stddev)     {         this.mean   = mean;         this.stddev = stddev;     }     /**      * Compute the next random value using the Central Limit Theorem,      * which states that the averages of sets of uniformly-distributed      * random values are normally distributed.      */     public float nextCentral()     {         // Average 12 uniformly-distributed random values.         float sum = 0.0f;         for (int j = 0; j < 12; ++j) sum += gen.nextFloat();         // Subtract 6 to center about 0.         return stddev*(sum - 6) + mean;     }     /**      * Compute the next randomn value using the polar algorithm.      * Requires two uniformly-distributed random values in [-1, +1).      * Actually computes two random values and saves the second one      * for the next invokation.      */     public float nextPolar()     {         // If there's a saved value, return it.         if (haveNextPolar) {             haveNextPolar = false;             return nextPolar;         }         float u1, u2, r;    // point coordinates and their radius         do {             // u1 and u2 will be uniformly-distributed             // random values in [-1, +1).             u1 = 2*gen.nextFloat() - 1;             u2 = 2*gen.nextFloat() - 1;             // Want radius r inside the unit circle.             r = u1*u1 + u2*u2;         } while (r >= 1);         // Factor incorporates the standard deviation.         float factor = (float) (stddev*Math.sqrt(-2*Math.log(r)/r));         // v1 and v2 are normally-distributed random values.         float v1 = factor*u1 + mean;         float v2 = factor*u2 + mean;         // Save v1 for next time.         nextPolar     = v1;         haveNextPolar = true;         return v2;     }     // Constants for the ratio algorithm.     private static final float C1 = (float) Math.sqrt(8/Math.E);     private static final float C2 = (float) (4*Math.exp(0.25));     private static final float C3 = (float) (4*Math.exp(-1.35));     /**      * Compute the next random value using the ratio algorithm.      * Requires two uniformly-distributed random values in [0, 1).      */     public float nextRatio()     {         float u, v, x, xx;         do {             // u and v are two uniformly-distributed random values             // in [0, 1), and u != 0.             while ((u = gen.nextFloat()) == 0);  // try again if 0             v = gen.nextFloat();             float y = C1*(v - 0.5f);    // y coord of point (u, y)             x = y/u;                    // ratio of point's coords             xx = x*x;         } while (               (xx > 5f - C2*u)                  // quick acceptance                     &&             ( (xx >= C3/u + 1.4f)             // quick rejection               (xx > (float) (-4*Math.log(u))) ) // final test         );         return stddev*x + mean;     } } 

Program 14 C1 demonstrates these algorithms. It generates 100,000 random values with each algorithm and counts how many values fall into each of 32 equal-width intervals. The program prints the results of each algorithm as a bar chart, in which the length of each bar for an interval is determined by that interval's count. The bars form bell curves. See Listing 14-1b.

Listing 14-1b A demonstration of three algorithms for generating normally distributed random values.
 package numbercruncher.program14_1; import numbercruncher.mathutils.RandomNormal; import numbercruncher.randomutils.Buckets; /**  * PROGRAM 14-1: Normally-Distributed Random Numbers  *  * Demonstrate algorithms for generating normally-distributed  * random numbers.  */ public class GenerateRandomNormal {     private static final int BUCKET_COUNT = 32;     private static final int NUMBER_COUNT = 100000;    // 100K     /** counters of random values that fall within each interval */     private Buckets buckets = new Buckets(BUCKET_COUNT);     /** generator of normally-distributed random values */     private RandomNormal normal;     /**      * Test the algorithms with a given mean and standard deviation.      * @param mean the mean      * @param stddev the standard deviation      */     private void run(float mean, float stddev)     {         long startTime;     // starting time of each algorithm         // Initialize the random number generator         // and the interval buckets.         normal = new RandomNormal();         normal.setParameters(mean, stddev);         buckets.setLimits(0, BUCKET_COUNT-1);         // Central limit theorem algorithm.         startTime = System.currentTimeMillis();         buckets.clear();         central();         print("Central", startTime);         // Polar algorithm.         startTime = System.currentTimeMillis();         buckets.clear();         polar();         print("Polar", startTime);         // Ratio algorithm.         startTime = System.currentTimeMillis();         buckets.clear();         ratio();         print("Ratio", startTime);     }     /**      * Print the results of an algorithm with its elapsed time.      * @param label the algorithm label      * @param startTime the starting time      */     private void print(String label, long startTime)     {         long elapsedTime = System.currentTimeMillis() - startTime;         System.out.println("\n" + label + " (" + elapsedTime +                            " ms):\n");         buckets.print();     }     /**      * Invoke the Central Limit Theorem algorithm.      */     private void central()     {         for (int i = 0; i < NUMBER_COUNT; ++i) {             buckets.put(normal.nextCentral());         }     }     /**      * Invoke the polar algorithm.      */     private void polar()     {         for (int i = 0; i < NUMBER_COUNT; ++i) {             buckets.put(normal.nextPolar());         }     }     /**      * Invoke the ratio algorithm.      */     private void ratio()     {         for (int i = 0; i < NUMBER_COUNT; ++i ) {             buckets.put(normal.nextRatio());         }     }     /**      * Main.      * @param args the array of program arguments      */     public static void main(String args[])     {         float mean   = BUCKET_COUNT/2f;         float stddev = BUCKET_COUNT/6f;         GenerateRandomNormal test = new GenerateRandomNormal();         test.run(mean, stddev);     } } 

Output:

 Central (710 ms):  0     87: *  1    164: *  2    313: **  3    413: ***  4    691: *****  5    994: *******  6   1430: **********  7   2005: **************  8   2572: ******************  9   3158: ********************** 10   3968: **************************** 11   4804: ********************************** 12   5508: *************************************** 13   6241: ******************************************** 14   6784: ************************************************ 15   7014: ************************************************** 16   7068: ************************************************** 17   7028: ************************************************** 18   6791: ************************************************ 19   6246: ******************************************** 20   5535: *************************************** 21   4840: ********************************** 22   3968: **************************** 23   3349: ************************ 24   2656: ******************* 25   1966: ************** 26   1472: ********** 27   1042: ******* 28    702: ***** 29    419: *** 30    283: ** 31    172: * Polar (490 ms):  0    117: *  1    190: *  2    271: **  3    431: ***  4    665: *****  5    929: ******  6   1432: **********  7   1989: **************  8   2478: *****************  9   3269: *********************** 10   3925: *************************** 11   4886: ********************************** 12   5509: ************************************** 13   6212: ******************************************* 14   6628: ********************************************** 15   7106: ************************************************* 16   7261: ************************************************** 17   7101: ************************************************* 18   6743: ********************************************** 19   6339: ******************************************** 20   5573: ************************************** 21   4816: ********************************* 22   3994: **************************** 23   3301: *********************** 24   2514: ***************** 25   1932: ************* 26   1391: ********** 27    960: ******* 28    690: ***** 29    443: *** 30    301: ** 31    192: * Ratio (550 ms):  0     98: *  1    183: *  2    283: **  3    454: ***  4    685: *****  5    995: *******  6   1428: **********  7   1933: *************  8   2497: *****************  9   3209: ********************** 10   4053: **************************** 11   4762: ********************************* 12   5502: ************************************** 13   6204: ******************************************* 14   6758: ********************************************** 15   7013: ************************************************ 16   7293: ************************************************** 17   7121: ************************************************* 18   6838: *********************************************** 19   6237: ******************************************* 20   5453: ************************************* 21   4802: ********************************* 22   3950: *************************** 23   3389: *********************** 24   2533: ***************** 25   2009: ************** 26   1342: ********* 27   1016: ******* 28    673: ***** 29    440: *** 30    277: ** 31    216: * 

The program also prints the elapsed time for each algorithm. The algorithm that uses the Central Limit Theorem is the simplest to program and understand, but it's the slowest because it requires the most uniformly distributed random values. The polar algorithm is the fastest of the three, and it is the algorithm used by the method java.util.Random.nextGaussian() .

The interval counters used by Program 14 C1 are implemented by class Buckets , which is shown in Listing 14-1c.

Listing 14-1c Class Buckets , which implements the interval counters.
 package numbercruncher.randomutils; import numbercruncher.mathutils.AlignRight; /**  * Counters of random values that fall within each interval.  */ public class Buckets {     private static final int MAX_BAR_SIZE = 50;     private AlignRight ar = new AlignRight();     /** number of intervals */      private int n;     /** counters per interval */    private int counters[];     /** minimum random value */     private float rMin;     /** maximum random value */     private float rMax;     /** from min to max */          private float width;     /**      * Constructor.      * @param n the number of intervals      */     public Buckets(int n)     {         this.n = n;         this.counters = new int[n];         clear();     }     /**      * Return the counter value for interval i.      * @param i the value of i      * @return the counter value      */     public int get(int i) { return counters[i]; }     /**      * Set the minimum and maximum random values.      * @param rMin the minimum value      * @param rMax the maximum value      */     public void setLimits(float rMin, float rMax)     {         this.rMin  = rMin;         this.rMax  = rMax;         this.width = (rMax - rMin)/n;     }     /**      * Determine a random value's interval and count it.      * @param r the random value      */     public void put(float r)     {         // Ignore the value if it's out of range.         if ((r < rMin)  (r > rMax)) return;         // Determine its interval and count it.         int i = (int) ((r - rMin)/width);         ++counters[i];     }     /**      * Clear all the interval counters.      */     public void clear()     {         for (int i = 0; i < counters.length; ++1) counters (i) =0;     }     /**      * Print the counter values as a horizontal bar chart. Scale the      * chart so that the longest bar is MAX_BAR_SIZE.      */     public void print()     {         // Get the longest bar's length.         int maxCount = 0;         for (int i = 0; i < n; ++i) {             maxCount = Math.max(maxCount, counters[i]);         }         // Compute the scaling factor.         float factor = ((float) MAX_BAR_SIZE)/maxCount;         // Loop to print each bar.         for (int i = 0; i < n; ++i) {             int b = counters[i];             // Interval number.             ar.print(i, 2);             ar.print(b, 7);             System.out.print(": ");             // Bar.             int length = Math.round(factor*b);             for (int j = 0; j < length; ++j) System.out.print("*");             System.out.println();         }     } } 

The screen shot in Screen 14-1 shows the interactive [1] version of Program 14 C1. With far fewer random values, the bell curve formed by the interval counts is much rougher. The interactive program also plots the bell curve defined by the normal distribution function

[1] Each interactive program in this book has two versions, an applet and a standalone application. You can download all the Java source code. See the downloading instructions in the preface of this book.

graphics/14equ05.gif


Screen 14-1. A screen shot of the interactive version of Program 14 C1.

graphics/14scr01.jpg

where m is the mean of the distribution and s is its standard deviation.


   
Top
 


Java Number Cruncher. The Java Programmer's Guide to Numerical Computing
Java Number Cruncher: The Java Programmers Guide to Numerical Computing
ISBN: 0130460419
EAN: 2147483647
Year: 2001
Pages: 141
Authors: Ronald Mak

flylib.com © 2008-2017.
If you may any questions please contact us: flylib@qtcs.net