From: ags@seaman.cc.purdue.edu (Dave Seaman) Newsgroups: sci.math.num-analysis Subject: Re: N-dim spherical random number drawing Date: 20 Sep 1996 13:04:58 -0500 In article , Shing-Te Li wrote: >Does anyone know a good algorithm that can draw a vector of random >numbers uniformly from a N dimensional sphere? Here are four methods for generating uniformly-distributed points on a unit sphere: (1) The normal-deviate method. Choose x, y, and z, each normally-distributed with mean 0 and variance 1. (Actually, the variance does not matter, provided it remains fixed.) Normalize the result to obtain a uniform density on the unit sphere. This method also generalizes nicely to n dimensions. It works because the vector chosen (before normalization) has a density that depends only on the distance from the origin. The method is mentioned in Knuth, _Seminumerical Algorithms_, 2nd. ed., pp. 130-131. (2) The hypercube rejection method. Choose x, y, and z, uniformly distributed on [-1,1]. Compute s = x^2 + y^2 + z^2. If s > 1, reject the triplet and start over. Once you have an acceptable vector, normalize it to obtain a point on the unit sphere. This method also generalizes to n dimensions, but as the dimension increases the probability of rejection becomes high and many random numbers are wasted. (3) The trig method. This method works only in 3-space, but it is very fast. It depends on the slightly counterintuitive fact (see proof below) that each of the three coordinates of a uniformly distributed point on S^2 is uniformly distributed on [-1,1] (but the three are not independent, obviously). Therefore, it suffices to choose one axis (Z, say) and generate a uniformly distributed value on that axis. This constrains the chosen point to lie on a circle parallel to the X-Y plane, and the obvious trig method may be used to obtain the remaining coordinates. (a) Choose z uniformly distributed in [-1,1]. (b) Choose t uniformly distributed on [0, 2*pi). (c) Let r = sqrt(1-z^2). (d) Let x = r * cos(t). (e) Let y = r * sin(t). (4) A variation of method (3) that doesn't use trig may be found in Knuth, _loc. cit._. The method is equivalent to the following: (a) Choose u and v uniformly distributed on [-1,1]. (b) Let s = u^2 + v^2. (c) If s > 1, return to step (a). (d) Let a = 2 * sqrt(1-s). (e) The desired point is (a*u, a*v, 2*s-1). This method uses two-dimensional rejection, which gives a higher probability of acceptance compared to algorithm (2) in three dimensions. I group this with the trig method because both depend on the fact that the projection on any axis is uniformly distributed on [-1,1], and therefore both methods are limited to use on S^2. I have found the trig method to be fastest in tests on an IBM RS/6000 model 590, using the XLF 3.2 Fortran compiler, with the IMSL subroutines DRNUN and DRNNOA generating the random numbers. The trig routines were accelerated by using the MASS library, and sqrt computations were performed by the hardware sqrt instruction available on the 590. Timings for 1 million random points on S^2: (1) normal-deviate method 9.60 sec. (2) 3-D rejection 8.84 sec. (3) trig method 2.71 sec. (4) polar with 2-D rejection 5.08 sec. The fact that algorithm (3) does not involve rejection turns out to be a huge advantage because it becomes possible to generate all the points in a single loop with no conditional statements, which optimizes very effectively and more than offsets the cost of using trig functions. The alternative with algorithms (2) and (4) is to generate the points one at a time and test each one for acceptance before proceeding to the next, which requires nested loops and many more subroutine calls to generate the random numbers. Algorithm (1) does not explicitly use rejection, but the IMSL subroutine DRNNOA uses rejection in generating the normally-distributed numbers that are required. In higher dimensions, the normal-deviate method is preferred because the probability of rejection in the hypercube method grows very rapidly with the dimension. -------------------------------------------------------------------- Here is a proof of correctness for the fact that the z-coordinate is uniformly distributed. The proof also works for the x- and y- coordinates, but note that this works only in 3-space. The area of a surface of revolution in 3-space is given by A = 2 * pi * int_a^b f(x) * sqrt(1 + [f'(x}]^2) dx Consider a zone of a sphere of radius R. Since we are integrating in the z direction, we have f(z) = sqrt(R^2 - z^2) f'(z) = -z / sqrt(R^2-z^2) 1 + [f'(z)]^2 = r^2 / (R^2-z^2) A = 2 * pi * int_a^b sqrt(R^2-z^2) * R/sqrt(R^2-z^2) dz = 2 * pi * R int_a^b dz = 2 * pi * R * (b-a) = 2 * pi * R * h, where h = b-a is the vertical height of the zone. Notice how the integrand reduces to a constant. The density is therefore uniform. -- Dave Seaman dseaman@purdue.edu ++++ stop the execution of Mumia Abu-Jamal ++++ ++++ if you agree copy these lines to your sig ++++ ++++ see http://www.xs4all.nl/~tank/spg-l/sigaction.htm ++++ ============================================================================== From: George Marsaglia Subject: Random points on a sphere. Date: Tue, 15 Jun 1999 18:03:07 -0400 Newsgroups: sci.math.num-analysis,sci.math,sci.stat.math Questions on methods for random sampling from the surface of a sphere keep coming up, and seemingly get lost as new generations of computer users encounter the problem. I described the methods that I had used for years in "Sampling from the surface of a sphere", Ann Math Stat, v43, 1972, recommending the two that seemed the fastest: 1) Choose standard normal variates x1,x2,x3, put R=1/sqrt(x1^2+x2^2+x3^2) and return the point (x1*R,x2*R,x3*R). 2) Generate u and v, uniform in [-1,1] until S=u^2+v^2 < 1, then return (1-2*S,u*r,v*r), with r=2*sqrt(1-S). The second method is fast and easiest to program---unless one already has a fast normal generator in his library. Method 1) seems by far the fastest (and easily applies to higher dimensions) if my latest normal generator is used. It produces normal variates at the rate of seven million per second in a 300MHz Pc. George Marsaglia