jroth@allvax.dec.com (Jim Roth) (05/03/90)
In article <49791@lanl.gov>, matthias@mpx1.lampf.lanl.gov (Matthias, Bjorn E.) writes... >This topic has certainly been beaten to death recently. Still, I believe >that there are some more useful things to say. .. stuff deleted ... >in Monte Carlo simulations to generate quantities with certain distributions. >Often it is not possible to invert the distribution in question and then one >resorts to rejection methods (e.g. von Neumann rejection) or to random walk >traversal algorithms (e.g. Metropolis algorithm). But when it is possible, >it is usually faster to do things directly as outlined above for this case. In this case, it is considerably faster to do a random walk on the sphere, because this avoids relatively expensive transcendental function calls; in fact the only library call besides the random number generator is a square root which could be done faster as an inverse square root if such a routine were available (or better if a normalize to unit vector routine were available.) Try the following code, which cyclically permutes the coordinates and performs a random plane rotation; the distribution of theta is not too important, so rejection is not needed here. The limit distribution is the same as long as it is not pathologically bad. x = sqrt(1.0/3.0) ! initial unit vector y = sqrt(1.0/3.0) z = sqrt(1.0/3.0) do i=1,succ cs = 1.0-2.0*ran(ix) sn = 1.0-2.0*ran(ix) t = 1.0/sqrt(cs**2+sn**2) cs = cs*t sn = sn*t t = x x = cs*y-sn*z y = sn*y+cs*z z = t ax(i)=x ay(i)=y az(i)=z enddo Needless to say, the rejection method could be tweaked by randomly selecting a sub-cube not completely outside the unit sphere and uniformly choosing a point within it, or some other table driven technique to improve the succes rate. It would approach the speed of the random walk then. However, for high dimensions the random walk is very attractive, since a plane rotation in one of n choose 2 planes can be done just as cheaply as in 2-D (as long as successive vectors need not be completely independant.) - Jim
matthias@mpx1.lampf.lanl.gov (Matthias, Bjorn E.) (05/04/90)
In a recent posting of mine, I was too hasty in claiming that random walk
sampling of a distribution in usually slower than a direct inversion when
the latter is possible. Sorry about that! A followup posting from
Jim Roth presented what I think is a very elegant random walk of sorts on
the surface of a sphere. He says:
>In this case, it is considerably faster to do a random walk on the sphere,
..
(his code is available in his posting)
..
I timed his random walk and checked the first and second moments of the
PHI and THETA distributions and thought I should post the results together
with those of the other two methods:
rejection method cpu = 31620 msec
failures = 476793 ; fraction = 47.68 %
successes = 523207 ; fraction = 52.32 %
diagnostics on distribution:
number of points = 523207
first moment of phi = 3.1417506
second moment of phi = 13.1581593
first moment of theta = 1.5707290
second moment of theta = 2.9352241
direct method cpu = 26050 msec
diagnostics on distribution:
number of points = 523207
first moment of phi = 3.1412365
second moment of phi = 13.1574106
first moment of theta = 1.5701435
second moment of theta = 2.9331174
random walk method cpu = 17780 msec
diagnostics on distribution:
number of points = 523207
first moment of phi = 3.1411185
second moment of phi = 13.1555815
first moment of theta = 1.5704634
second moment of theta = 2.9331996
Clearly, the random walk is the superior performer here.
Thanks, Jim!
Regards,
Bjorn
---
Bjorn E. Matthias MATTHIAS@LAMPF (BITNET)
matthias@mpx0.lampf.lanl.gov or matthias@mpx1.lampf.lanl.gov (internet)dvb@Apple.COM (David Van Brink) (05/07/90)
A mathemetician friend of mine worked out an absurdly simple means of finding a random point on the surface of a unit radius sphere. 1. Pick a random angle, 0 <= theta < 360. This is in the x-y plane. 2. Pick a random z-value from -1 to 1 ( "<=" or "<" ? I'm not sure). Turns out the distribution of area is linear along an axis. I was somewhat skeptical of this, but at least one non-disproof is that the surface area of the sphere is the same as the surface area of the curved part of a bounding cylinder. (4*pi*r^2).
rick@hanauma.stanford.edu (Richard Ottolini) (05/07/90)
In article <40768@apple.Apple.COM+ dvb@Apple.COM (David Van Brink) writes:
+
+A mathemetician friend of mine worked out an absurdly
+simple means of finding a random point on the surface
+of a unit radius sphere.
+
+1. Pick a random angle, 0 <= theta < 360. This is in
+the x-y plane.
+
+2. Pick a random z-value from -1 to 1 ( "<=" or "<" ?
+I'm not sure).
Huh? Each latitude will have the same fraction of points and things
will appear more crowded at the poles.gbrown@tybalt.caltech.edu (Glenn C. Brown) (05/07/90)
rick@hanauma.stanford.edu (Richard Ottolini) writes: >+1. Pick a random angle, 0 <= theta < 360. This is in >+the x-y plane. >+ >+2. Pick a random z-value from -1 to 1 ( "<=" or "<" ? >+I'm not sure). >Huh? Each latitude will have the same fraction of points and things ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ >will appear more crowded at the poles. Latitude is the angle of elevation (usually refered to as phi) in polar coordinates. Choosing a random z is not equivalent to choosing a random phi. Therefore, each latitude will not have the same fraction of points. I believe the posted method may be right. If it is, the proof should be straightforward (I just completed a course in surface integrals), and I'll post it in a few days if noone else beats me to it. (That's assuming the above is true.) --Glenn
kneller@cgl.ucsf.edu (Don Kneller) (05/09/90)
gbrown@tybalt.caltech.edu (Glenn C. Brown) writes: >rick@hanauma.stanford.edu (Richard Ottolini) writes: >>+1. Pick a random angle, 0 <= theta < 360. This is in >>+the x-y plane. >>+ >>+2. Pick a random z-value from -1 to 1 ( "<=" or "<" ? >>+I'm not sure). >I believe the posted method may be right. I believe so too. I had previously posted that you should pick phi such that: phi = arccos(random number between -1 and 1) and since z is R * cos(phi), this becomes: z = random number between -R and R, which is what Richard proposed. To see this is indeed correct, look at the area of a strip of height "dz" at angle phi on a sphere of radius R. area = circumference * width of strip on sphere = 2 pi R sin(phi) * R d(phi) as: z = R cos(phi), dz = - R sin(phi) d(phi) so area = 2 pi R dz (don't worry about the sign, dz and d(phi) are in opposite directions) That is, the area on the sphere of a strip of height "dz" is *independent* of z. If the linear density of points along z is rho, then N = rho * dz will fall between z and z + dz. The density of points on the sphere will be: rho2 = N / area = rho * dz / (2 pi R dz) = rho / (2 pi R) The bottom line?: 1) pick theta randomly between -pi and pi 2) pick z randomly between -R and R ----- Don Kneller UUCP: ...ucbvax!ucsfcgl!kneller INTERNET: kneller@cgl.ucsf.edu
kneller@socrates.ucsf.edu (Don Kneller) (05/09/90)
rick@hanauma.stanford.edu (Richard Ottolini) writes: >+1. Pick a random angle, 0 <= theta < 360. This is in >+the x-y plane. >+ >+2. Pick a random z-value from -1 to 1 ( "<=" or "<" ? >+I'm not sure). I believe this is correct. I had previously posted that you should pick phi such that: phi = arccos(random number between -1 and 1) and since z is R * cos(phi), this becomes: z = random number between -R and R, which is what Richard proposed. To see this is indeed correct, look at the area of a strip of height "dz" at angle phi on a sphere of radius R. area = circumference * width of strip on sphere = 2 pi R sin(phi) * R d(phi) as: z = R cos(phi), dz = - R sin(phi) d(phi) so area = 2 pi R dz (don't worry about the sign, dz and d(phi) are in opposite directions) That is, the area on the sphere of a strip of height "dz" is *independent* of z. If the linear density of points along z is rho, then N = rho * dz will fall between z and z + dz. The density of points on the sphere will be: rho2 = N / area = rho * dz / (2 pi R dz) = rho / (2 pi R) , independent of z. The bottom line?: 1) pick theta randomly between -pi and pi 2) pick z randomly between -R and R ----- Don Kneller UUCP: ...ucbvax!ucsfcgl!kneller INTERNET: kneller@cgl.ucsf.edu BITNET: kneller@ucsfcgl.BITNET
cca@newton.physics.purdue.edu (Charles C. Allen) (05/09/90)
> 1) pick theta randomly between -pi and pi > 2) pick z randomly between -R and R The above method (azimuthal angle between -pi and pi, z between -R and R) works fine, but it does not save you a trig function call unless that's all you need. You have one arccos call if you want the polar angle, and at least one more plus some square roots if you want the (x,y,z) triple. Which method is "fastest" depends not only on what hardware you have, but also what you want for a final result. [Minor personal note - every single math & physics text I have uses theta for the polar angle and phi for the azimuthal angle when discussing spherical coordinates. Perhaps my brain has calcified, but it *is* a little disconcerting to realize that a posting isn't making sense because someone is using a different notation.] Charles Allen Internet: cca@newton.physics.purdue.edu Department of Physics HEPnet: purdnu::allen, fnal::cca Purdue University talknet: 317/494-9776 West Lafayette, IN 47907