cj@hpsad.HP.COM (Chris Johnson) (09/21/89)
Hi, I'm looking for a random number routine which will generate numbers with a Gaussian distribution. Rand(3C) and other such programs produce (semi) even distributions. I suppose some sort of conversion routine would also work. If anyone out there has such a beast, please let me know. -cj cj@hpsad.hp.com
dann@hpsad.HP.COM (Dan Needham) (09/22/89)
>>>> cj@hpsad.HP.COM (Chris Johnson) writes ...... >>>> I'm looking for a random number routine which will generate >>>> numbers with a Gaussian distribution. Rand(3C) and other such >>>> programs produce (semi) even distributions. I suppose some sort >>>> of conversion routine would also work. If anyone out there has >>>> such a beast, please let me know. -cj cj@hpsad.hp.com ---------- I had to do this once in school and still have the paper I wrote it down on. This part of the program was handed out by the professor. function [a,b] = grv(sigma) % Generates a pair of Gaussian random variables with mean zero, variance sigma rand('uniform'); % Generates #'s between [0,1] pi = 3.1415927; p2 = 2*pi*rand(1); p1 = (-2) * log(rand(1)); a = sigma * sqrt(p1) * cos(p2); b = sigma * sqrt(p1) * sin(p2); ------------------------------------------------------------ Dan Needham
cj@hpsad.HP.COM (Chris Johnson) (09/23/89)
Thanks to all of you who sent replies (both to notes and directly) I have what I need. If anyone wants some of these, send me a note and I'll ship then off to you. -cj cj@hpsad.hp.com
ken@aiai.ed.ac.uk (Ken Johnson) (09/24/89)
> I'm looking for a random number routine which will generate > numbers with a Gaussian distribution. Rand(3C) and other such > programs produce (semi) even distributions. Assuming that rand( ) returns a rectangular distribution: [this is untested and based on what I remember of the old SPSS library] gauss(mean,sd) /* sd: standard deviation */ float mean; float sd; { float t; float rand( ); int i; t = 0.0; for (i = 0; i < 12; ++i) { t += rand( ); } return(sd * ((t - 6.0) / 12.0)); } The 12 is a magic number, so don't change it. Summing 12 rectangular distributions gives a gaussian one. If anyone can explain why clearly enough for me to understand, I will send them a Beer Token. -- Ken Johnson, AI Applications Institute, 80 South Bridge, Edinburgh EH1 1HN E-mail ken@aiai.ed.ac.uk, phone 031-225 4464 extension 212 `I have read your article, Mr. Johnson, and I am no wiser than when I started.' -- `Possibly not, sir, but far better informed.'
karl@haddock.ima.isc.com (Karl Heuer) (09/25/89)
In article <927@skye.ed.ac.uk> ken@aiai.UUCP (Ken Johnson) writes: >>[wants random numbers with Gaussian distribution] >[includes code to sum 12 uniform-distribution random numbers] First, the standard rand() function is uniform on a range of *integers*, so you'd have to replace rand() with rand()/(RAND_MAX+1.0) to get what you want. Second, the resulting distribution is just an approximation, not a true Gaussian distribution. (All its values are in the range [-6,6], whereas a true Gaussian goes all the way to infinity.) It may be close enough for your application, but you should be aware of the differences. Third, a benchmark once revealed that it can be *slower* than the (already posted) mathematically exact Gaussian that you can get from two uniform variables and sin/cos/log. (Your results may vary.) >The 12 is a magic number... If anyone can explain why clearly enough for me >to understand, I will send them a Beer Token. It's because 12 is the number in a dozen. For bakery applications, e.g. if this code is to be used in a toast-flipping simulation of the Murphy's Law phenomenon, one should use 13. (<-- Culture-dependent joke.) But seriously, the reason is simply that the uniform distribution on [0,1) has a variance of exactly 1/12. Karl W. Z. Heuer (ima!haddock!karl or karl@haddock.isc.com), The Walking Lint
mouse@mcgill-vision.UUCP (der Mouse) (10/01/89)
In article <1820002@hpsad.HP.COM>, cj@hpsad.HP.COM (Chris Johnson) writes: > I'm looking for a random number routine which will generate numbers > with a Gaussian distribution. Rand(3C) and other such programs > produce (semi) even distributions. They also produce integers. If you have a routine that produces floating-point random numbers in the range [0,1), Knuth has an algorithm which produces nice Gaussian numbers. Here's an extract from the comment header for the Gaussian random number routine from a local library: * Algorithm is the polar method, from Knuth _Art_Of_Computer_Programming_ * Volume 2 (Seminumerical Algorithms), chapter 3 (random numbers). * This is Algorithm P in section 3.4.1. * * Algorithm P (Polar method for normal deviates) * This algorithm calculates two independent normally distributed * variables X1 and X2, given two independent uniformly distributed * variables U1 and U2. * * P1. [Get uniform variables.] * Generate U1 and U2. * Set V1 <- 2U1 - 1, V2 <- 2U2-1. * P2. [Compute S] * Set S <- (V1*V1) + (V2*V2) * P3. [Is S >= 1] * If S >= 1, return to step P1. * P4. [Compute X1, X2] * Set X1 <- V1 * t, X2 <- V2 * t, where t is sqrt((-2 ln(S))/S) The code itself is in assembly, which is why I'm not posting it. There are two points that perhaps need comment. One is that the number are generated in pairs. Knuth includes a proof that if the underlying uniform generator produces independent numbers when called to generate pairs this way, the resulting Gaussian numbers will be independent. The other point is that the above algorithm is a theoretical algorithm. In particular, S==0 can occur with non-zero probability in an implementation, while the above algorithm appears to assume that since the probability of S being 0 is mathematically zero, the issue can be ignored. (Note that the computation of t involves dividing by S.) It also may be unacceptably expensive; a square root and a natural logarithm are required for each pair of numbers. Depending on what you want them for, you may be able to cheat and just add together a handful of uniform numbers. This produces near-Gaussian distributions relatively cheaply. (Adding N independent uniform [0,1] numbers together produces a pseudo-Gaussian of mean N/2 and variance N/12. A reasonable choice of N is 12, thus automatically producing a variance of 1.) These pseudo-Gaussians are missing the long tails, but for some applications (eg, adding Gaussian noise to a picture) the tails don't matter much. der Mouse old: mcgill-vision!mouse new: mouse@larry.mcrcim.mcgill.edu