7.2 Transformation Method:Exponential and Normal Deviates 289 An important example of the use of (7.2.8)is the Box-Muller method for generating random deviates with a normal(Gaussian)distribution, 1 p(y)dy _e-v'/2dy V2 (7.2.9) Consider the transformation between two uniform deviates on(0,1),z1,z2,and two quantities y1,2, y1 =V-21n 1 cos 2Tx2 (7.2.10) 2=√-2lnx1sin2xx2 Equivalently we can write 3 (7.2.11) RECIPES x2= 2 arctan 1 Now the Jacobian determinant can readily be calculated (try it!) 气 Press. 1 1 0(x1,x2】 /2 OT. (7.2.12) 0(y1,y2 Programs Since this is the product of a function of y2 alone and a function of y alone,we see IENTIFIC that each y is independently distributed according to the normal distribution(7.2.9). 6 One further trick is useful in applying(7.2.10).Suppose that,instead of picking uniform deviates and 2 in the unit square,we instead pick v and v2 as the ordinate and abscissa of a random point inside the unit circle around the origin.Then the sum of their squares,R2=v is a uniform deviate,which can be used for1, 1920 COMPUTING (ISBN while the angle that(v1,v2)defines with respect to the v axis can serve as the random angle 22.What's the advantage?It's that the cosine and sine in(7.2.10)can now Numerical Recipes 10-621 be written as v/vR2 and v2/vR2,obviating the trigonometric function calls! We thus have 43106 #include <math.h> (outside float gasdev(long *idum) Software. Returns a normally distributed deviate with zero mean and unit variance,using ran1(idum) North as the source of uniform deviates. float ran1(long *idum); static int iset=0; static float gset; float fac,rsq,v1,v2; 1f(*1dum<0)1set=0: Reinitialize. if (iset ==0){ We don't have an extra deviate handy,so do v1=2.0*ran1(1dum)-1.0; pick two uniform numbers in the square ex- v2=2.0*ran1(1dum)-1.0; tending from-1 to +1 in each direction, rsq=v1*v1+v2*v2; see if they are in the unit circle,7.2 Transformation Method: Exponential and Normal Deviates 289 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). An important example of the use of (7.2.8) is the Box-Muller method for generating random deviates with a normal (Gaussian) distribution, p(y)dy = 1 √ 2π e−y2/2dy (7.2.9) Consider the transformation between two uniform deviates on (0,1), x 1, x2, and two quantities y1, y2, y1 = −2 ln x1 cos 2πx2 y2 = −2 ln x1 sin 2πx2 (7.2.10) Equivalently we can write x1 = exp  −1 2 (y2 1 + y2 2)  x2 = 1 2π arctany2 y1 (7.2.11) Now the Jacobian determinant can readily be calculated (try it!): ∂(x1, x2) ∂(y1, y2) = ∂x1 ∂y1 ∂x1 ∂y2 ∂x2 ∂y1 ∂x2 ∂y2 = −  1 √ 2π e−y2 1/2   1 √ 2π e−y2 2/2  (7.2.12) Since this is the product of a function of y2 alone and a function of y1 alone, we see that each y is independently distributed according to the normal distribution (7.2.9). One further trick is useful in applying (7.2.10). Suppose that, instead of picking uniform deviates x1 and x2 in the unit square, we instead pick v1 and v2 as the ordinate and abscissa of a random point inside the unit circle around the origin. Then the sum of their squares, R2 ≡ v2 1 +v2 2 is a uniform deviate, which can be used for x1, while the angle that(v1, v2) defines with respect to the v1 axis can serve as the random angle 2πx2. What’s the advantage? It’s that the cosine and sine in (7.2.10) can now be written as v1/ √ R2 and v2/ √ R2, obviating the trigonometric function calls! We thus have #include <math.h> float gasdev(long *idum) Returns a normally distributed deviate with zero mean and unit variance, using ran1(idum) as the source of uniform deviates. { float ran1(long *idum); static int iset=0; static float gset; float fac,rsq,v1,v2; if (*idum < 0) iset=0; Reinitialize. if (iset == 0) { We don’t have an extra deviate handy, so do { v1=2.0*ran1(idum)-1.0; pick two uniform numbers in the square ex￾v2=2.0*ran1(idum)-1.0; tending from -1 to +1 in each direction, rsq=v1*v1+v2*v2; see if they are in the unit circle,