1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
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
|
