diff options
Diffstat (limited to 'reference/sphere.rand')
| -rw-r--r-- | reference/sphere.rand | 157 |
1 files changed, 157 insertions, 0 deletions
diff --git a/reference/sphere.rand b/reference/sphere.rand new file mode 100644 index 0000000..bf56705 --- /dev/null +++ b/reference/sphere.rand @@ -0,0 +1,157 @@ +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 <wmEgVDG00VUtMEhUwg@andrew.cmu.edu>, +Shing-Te Li <sl2x+@andrew.cmu.edu> 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 <geo@stat.fsu.edu> +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 |