ken@ut-emx.UUCP (kenneth moore) (02/14/89)
Hi., I am looking for a program to generate -Normal distributed random numbers (for any mean and standard deviation).I would be very thankful if any one can suggest a good source code or send me one. Thanks., Kenneth Moore e-mail: (ken@emx.utexas.edu)
toms@ncifcrf.gov (Tom Schneider) (02/14/89)
In article <10471@ut-emx.UUCP> ken@ut-emx.UUCP (kenneth Moore) writes: >I am looking for a program to generate -Normal distributed > random numbers (for any mean and standard deviation).I would > be very thankful if any one can suggest a good source code > or send me one. > e-mail: (ken@emx.utexas.edu) The simplest way to generate such numbers is to add together several flat distribution numbers. This is very costly, and there is a better way. A friend showed me how to do this last week, so I put the code together today. It works fine. The program is in Pascal, and should work on any machine with a good compiler. Send me a note if you want a copy. It's not hard to do. If U is a member of the set [0..1] and Un and Un+1 are two members, then define theta = Un 2 pi r = sqrt(-2 ln(Un+1)) then when these polar coordinates are converted to Cartesian coordinates, one gets two independent Normally distributed numbers with mean 0 and standard deviation 1. To get other standard deviations multiply by a constant, and to get other means, add a constant. I don't know a reference, but I'm sure there are many. Tom Schneider National Cancer Institute Laboratory of Mathematical Biology Frederick, Maryland toms@ncifcrf.gov
bs@linus.UUCP (Robert D. Silverman) (02/14/89)
In article <734@fcs280s.ncifcrf.gov> toms@ncifcrf.gov (Tom Schneider) writes: >In article <10471@ut-emx.UUCP> ken@ut-emx.UUCP (kenneth Moore) writes: >>I am looking for a program to generate -Normal distributed >> random numbers (for any mean and standard deviation).I would >The simplest way to generate such numbers is to add together several >flat distribution numbers. This is very costly, and there is a better way. >A friend showed me how to do this last week, so I put the code together Your friend misinformed you about the cost of the method. See below. >today. It works fine. The program is in Pascal, and should work >on any machine with a good compiler. Send me a note if you want a copy. >It's not hard to do. If U is a member of the set [0..1] and Un and Un+1 >are two members, then define > theta = Un 2 pi > r = sqrt(-2 ln(Un+1)) stuff deleted.... This is essentially the polar method. It does work. Your assertion that it is less costly than using the central limit theorem is incorrect however. It is far faster to do 6 floating point additions than to compute Un followed by a LOG and a SQUARE ROOT. The latter computation is MUCH slower. Especially on vector machines where a vectorized add is very quick. I might also suggest looking at rejection methods and Marsaglia's wedge tail method. Check out Knuth V. 2 Bob Silverman
pokey@well.UUCP (Jef Poskanzer) (02/14/89)
In the referenced message, toms@ncifcrf.gov (Tom Schneider) wrote: }In article <10471@ut-emx.UUCP> ken@ut-emx.UUCP (kenneth Moore) writes: }>I am looking for a program to generate -Normal distributed }> random numbers } }The simplest way to generate such numbers is to add together several }flat distribution numbers. This is very costly, and there is a better way. } theta = Un 2 pi } r = sqrt(-2 ln(Un+1)) }then when these polar coordinates are converted to Cartesian }coordinates, one gets two independent Normally distributed numbers Excuse me, but for what meaning of the word "costly" is generating several uniform numbers even SLIGHTLY more costly than generating two uniform numbers and then doing a log, a square root, and a sine or cosine?? Not to mention all those multiplies. I prefer the polar method too (see below for an implementation in C, straight out of Knuth), but not because of "cost". I use it because it's *correct*. The "add up N uniforms" method is only an approximation. --- Jef Jef Poskanzer jef@helios.ee.lbl.gov ...well!pokey "No wonder white people are going to Hell..." -- Andrew Purshottam float rndnor( mean, stddev ) float mean, stddev; { float v1, v2, z; do { v1 = -log( 1.0 - rnd() ); v2 = -log( 1.0 - rnd() ); } while ( 2.0 * v1 < ( (v2-1.0)*(v2-1.0) ) ); if ( rnd() > 0.5 ) z = 1.0; else z = -1.0; return stddev * z * v2 + mean; }
rjfrey@kepler1.UUCP (Robert J Frey) (02/18/89)
In article <734@fcs280s.ncifcrf.gov> toms@ncifcrf.gov (Tom Schneider) writes: >In article <10471@ut-emx.UUCP> ken@ut-emx.UUCP (kenneth Moore) writes: >>I am looking for a program to generate -Normal distributed > >...I don't know a reference, but I'm sure there are many. > A good reference for this sort of thing is George Fishman, Principles of Discrete Event Simulation. This is an excellent text and covers many other issues relevant here, such as appropriate statistical test a RNG should pass. There are FORTRAN-like psuedo-code descriptions of many RNG's; most of the algorithms are quite simple. ============================================================================== |Dr. Robert J. Frey | {acsm, sbcs, polyof}!kepler1!rjfrey | |Kepler Financial Management, Ltd.|------------------------------------------| |100 North Country Rd., Bldg. B | The views expressed are wholly my own and| |Setauket, NY 11766 | and do not reflect those of the Indepen- | |(516) 689-6300 x.16 | dent Republic of Latvia. | ==============================================================================
news@pyr.gatech.EDU (News Administrator) (02/26/89)
The algorithm gains its speed by having a fast part 99% of the time (a multiplication and a comparison) and a slow part 1% of the time. This algorithm is not suited to vector machines though because of the condition checks. In fact any rejection method suffers from the same problem. I have coded this algorithm bot in Fortran and 8086 assembler and I'll mail them upon request. One gains a very good knowledge about the cost of functions like Log etc, even when you have a co-processor by using the probability
integral transformation to generate an exponetial and then generating the same by the above method.
B. Narasimhan e-mail: naras@stat.fsu.edu on internet.
From: naras@stat.uucp (B. Narasimhan)
Path: stat!naras
--
News Administrator -- Don Deal
Georgia Institute of Technology, Atlanta GA 30332
...!{akgua,allegra,amd,hplabs,ihnp4,masscomp,ut-ngp}!gatech!gitpyr!news