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 ++++