Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Generating random numbers on unit sphere

Subject: Generating random numbers on unit sphere

From: kanika

Date: 26 Jun, 2007 15:48:25

Message: 1 of 31

Hi everyone,

I am trying to generate random numbers on a unit sphere and further in
higher dimensions also. Any suggestions to generate these
efficiently? I know that numerical recipes should have a code for
this in C, which I cannot find. Any suggestions or pointers or
mathematical tips for this?

Thank you
Kanika

Subject: Generating random numbers on unit sphere

From: Per S

Date: 26 Jun, 2007 12:21:03

Message: 2 of 31

kanika wrote:
>
>
> Hi everyone,
>
> I am trying to generate random numbers on a unit sphere and further
> in
> higher dimensions also. Any suggestions to generate these
> efficiently? I know that numerical recipes should have a code for
> this in C, which I cannot find. Any suggestions or pointers or
> mathematical tips for this?
>
> Thank you
> Kanika
>
help sph2cart
%note I may have taken angles and radii in wrong order below
n=1000;
theta=pi*rand(1,n);
phi=2*pi*rand(1,n);
[x,y,z]=sph2cart(theta,phi,ones(1,n));
plot3(x,y,z,'.');

For hyper-spheres, you have to put hyper-angles randomly and then use
the appropiate transform to get x,y,z,w,... etc
/Per

Subject: Generating random numbers on unit sphere

From: Peter Boettcher

Date: 26 Jun, 2007 12:49:21

Message: 3 of 31

kanika <kanika.dhyani@gmail.com> writes:

> Hi everyone,
>
> I am trying to generate random numbers on a unit sphere and further in
> higher dimensions also. Any suggestions to generate these
> efficiently? I know that numerical recipes should have a code for
> this in C, which I cannot find. Any suggestions or pointers or
> mathematical tips for this?

The simple solution I like is to use randn to generate a set of
N-dimensional Gaussian numbers, then normalize each one. Since
multidimensional Gaussians are radially symmetric, this gives the
correct density on the unit (hyper-)sphere.

-Peter

Subject: Generating random numbers on unit sphere

From: ellieandrogerxyzzy@mindspring.com.invalid (Roger Stafford)

Date: 26 Jun, 2007 20:30:27

Message: 4 of 31

In article <ef5bc41.0@webcrossing.raydaftYaTP>, "Per S"
<per.sundqvist@uam.es> wrote:

> kanika wrote:
> > ........
> > I am trying to generate random numbers on a unit sphere and further
> > ........
> help sph2cart
> %note I may have taken angles and radii in wrong order below
> n=1000;
> theta=pi*rand(1,n);
> phi=2*pi*rand(1,n);
> [x,y,z]=sph2cart(theta,phi,ones(1,n));
> plot3(x,y,z,'.');
>
> For hyper-spheres, you have to put hyper-angles randomly and then use
> the appropiate transform to get x,y,z,w,... etc
> /Per
-------------------
  The difficulty with that distribution, Per, is that it is not uniformly
distributed area-wise over the unit sphere's surface. It packs in the
points more densely at its two poles.

  Matlab's 'sph2cart' function uses elevation as its second argument, phi,
which ranges from -pi/2 to +pi/2 and azimuth for its first argument,
theta, ranging from 0 to 2*pi. The differential for area would be

 dA = 1^2*cos(phi)*dphi*dtheta,

and as you see, dA becomes vanishingly small at the two poles where
cos(phi) approaches zero. This would cause a theoretically infinite
density there, using straight 'rand' outputs for phi.

  What you need to generate a uniform area distribution using 'rand' is:

 theta = 2*pi*rand(1,n);
 phi = asin(2*rand(1,n)-1);

  In N dimensions life gets increasingly difficult. The method mentioned
by Peter using 'randn' is the one I always use.

Roger Stafford

Subject: Generating random numbers on unit sphere

From: kanika

Date: 27 Jun, 2007 08:53:22

Message: 5 of 31

Thank you Per , Peter and Roger for the suggestions...

Another question that comes to my mind related to this...is depending
on the dimension(d) how many random points would be enough to cover
(approximate)
a (hyper)sphere. For the moment I use something like 10^d points as a
good number of random points.

Kanika

On Jun 26, 10:30 pm, ellieandrogerxy...@mindspring.com.invalid (Roger
Stafford) wrote:
> In article <ef5bc4...@webcrossing.raydaftYaTP>, "Per S"
>
>
>
> <per.sundqv...@uam.es> wrote:
> >kanikawrote:
> > > ........
> > > I am trying to generate random numbers on a unit sphere and further
> > > ........
> > help sph2cart
> > %note I may have taken angles and radii in wrong order below
> > n=1000;
> > theta=pi*rand(1,n);
> > phi=2*pi*rand(1,n);
> > [x,y,z]=sph2cart(theta,phi,ones(1,n));
> > plot3(x,y,z,'.');
>
> > For hyper-spheres, you have to put hyper-angles randomly and then use
> > the appropiate transform to get x,y,z,w,... etc
> > /Per
>
> -------------------
> The difficulty with that distribution, Per, is that it is not uniformly
> distributed area-wise over the unit sphere's surface. It packs in the
> points more densely at its two poles.
>
> Matlab's 'sph2cart' function uses elevation as its second argument, phi,
> which ranges from -pi/2 to +pi/2 and azimuth for its first argument,
> theta, ranging from 0 to 2*pi. The differential for area would be
>
> dA = 1^2*cos(phi)*dphi*dtheta,
>
> and as you see, dA becomes vanishingly small at the two poles where
> cos(phi) approaches zero. This would cause a theoretically infinite
> density there, using straight 'rand' outputs for phi.
>
> What you need to generate a uniform area distribution using 'rand' is:
>
> theta = 2*pi*rand(1,n);
> phi = asin(2*rand(1,n)-1);
>
> In N dimensions life gets increasingly difficult. The method mentioned
> by Peter using 'randn' is the one I always use.
>
> Roger Stafford

Subject: Generating random numbers on unit sphere

From: John D'Errico

Date: 27 Jun, 2007 06:21:22

Message: 6 of 31

kanika wrote:
>
>
> Thank you Per , Peter and Roger for the suggestions...
>
> Another question that comes to my mind related to this...is
> depending
> on the dimension(d) how many random points would be enough to cover
> (approximate)
> a (hyper)sphere. For the moment I use something like 10^d points
> as a
> good number of random points.
>
> Kanika

How much of a cover is "enough"? How good
of an approximation do you need?

You could use the "area" of the hypersphere
as a measure, so the surface area of such a
sphere would be:

  2*pi^(n/2)*r^(n-1)/gamma(n/2)

But only you can decide how many points
you need.

John

Subject: Generating random numbers on unit sphere

From: ellieandrogerxyzzy@mindspring.com.invalid (Roger Stafford)

Date: 27 Jun, 2007 11:16:36

Message: 7 of 31

In article <1182934402.356198.245770@k79g2000hse.googlegroups.com>, kanika
<kanika.dhyani@gmail.com> wrote:

> Thank you Per , Peter and Roger for the suggestions...
>
> Another question that comes to my mind related to this...is depending
> on the dimension(d) how many random points would be enough to cover
> (approximate)
> a (hyper)sphere. For the moment I use something like 10^d points as a
> good number of random points.
>
> Kanika
----------------
  Whether 10^d random points are appropriate or not depends on what you
intend to do with your results. We have no way of judging that for you.

  However, the fact that you have written 10^d rather than 10^(d-1) makes
me wonder whether, when you originally said "random numbers on a unit
sphere", you actually meant random points uniformly distributed in the
volume of the unit sphere rather than on its surface. In n-dimensional
space the surface of the unit sphere is only (n-1)-dimensional, whereas
its volume is n-dimensional. If it is volume, even n-dimensional volume,
you are seeking to fill, this makes it a quite different problem. I wrote
a function called 'randsphere' which fills the interior of an
n-dimensional hypersphere uniformly with random points. You can find it
in the matlab central file exchange at:

<http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=9443&objectType=file>

Roger Stafford

Subject: Generating random numbers on unit sphere

From: kanika

Date: 27 Jun, 2007 15:27:17

Message: 8 of 31


Well I meant points on the surface of a unit sphere and not in the
interior.
But it is clear to me now.

Thank you
Kanika

Subject: Generating random numbers on unit sphere

From: Per S

Date: 27 Jun, 2007 11:48:41

Message: 9 of 31

Roger Stafford wrote:
>
>
> In article <ef5bc41.0@webcrossing.raydaftYaTP>, "Per S"
> <per.sundqvist@uam.es> wrote:
>
>> kanika wrote:
>> > ........
>> > I am trying to generate random numbers on a unit sphere
and
> further
>> > ........
>> help sph2cart
>> %note I may have taken angles and radii in wrong order below
>> n=1000;
>> theta=pi*rand(1,n);
>> phi=2*pi*rand(1,n);
>> [x,y,z]=sph2cart(theta,phi,ones(1,n));
>> plot3(x,y,z,'.');
>>
>> For hyper-spheres, you have to put hyper-angles randomly and
then
> use
>> the appropiate transform to get x,y,z,w,... etc
>> /Per
> -------------------
> The difficulty with that distribution, Per, is that it is not
> uniformly
> distributed area-wise over the unit sphere's surface. It packs in
> the
> points more densely at its two poles.
>
> Matlab's 'sph2cart' function uses elevation as its second
> argument, phi,
> which ranges from -pi/2 to +pi/2 and azimuth for its first
> argument,
> theta, ranging from 0 to 2*pi. The differential for area would be
>
> dA = 1^2*cos(phi)*dphi*dtheta,
>
> and as you see, dA becomes vanishingly small at the two poles where
> cos(phi) approaches zero. This would cause a theoretically
> infinite
> density there, using straight 'rand' outputs for phi.
>
> What you need to generate a uniform area distribution using
> 'rand' is:
>
> theta = 2*pi*rand(1,n);
> phi = asin(2*rand(1,n)-1);
>
> In N dimensions life gets increasingly difficult. The method
> mentioned
> by Peter using 'randn' is the one I always use.
>
> Roger Stafford
>
  Right about that you are. The thing with randn, must be that you
have to define both an inner and an outer radii to get the points on
the surface?

xyz=randn(3,1000);
rad2=xyz(1,:).^2+xyz(2,:).^2+xyz(3,:).^2;
r1=0.99^2;
r2=1.01^2;
ix=find(and(rad2>=r1,rad2<=r2));

if r1->r2 then you need infinte number of points! Or have I
missunderstood you?

/Per

Subject: Generating random numbers on unit sphere

From: Per S

Date: 27 Jun, 2007 12:02:44

Message: 10 of 31

Per S wrote:
>
>
> Roger Stafford wrote:
>>
>>
>> In article <ef5bc41.0@webcrossing.raydaftYaTP>, "Per S"
>> <per.sundqvist@uam.es> wrote:
>>
>>> kanika wrote:
>>> > ........
>>> > I am trying to generate random numbers on a unit
sphere
> and
>> further
>>> > ........
>>> help sph2cart
>>> %note I may have taken angles and radii in wrong order
below
>>> n=1000;
>>> theta=pi*rand(1,n);
>>> phi=2*pi*rand(1,n);
>>> [x,y,z]=sph2cart(theta,phi,ones(1,n));
>>> plot3(x,y,z,'.');
>>>
>>> For hyper-spheres, you have to put hyper-angles randomly
and
> then
>> use
>>> the appropiate transform to get x,y,z,w,... etc
>>> /Per
>> -------------------
>> The difficulty with that distribution, Per, is that it is not
>> uniformly
>> distributed area-wise over the unit sphere's surface. It packs
> in
>> the
>> points more densely at its two poles.
>>
>> Matlab's 'sph2cart' function uses elevation as its second
>> argument, phi,
>> which ranges from -pi/2 to +pi/2 and azimuth for its first
>> argument,
>> theta, ranging from 0 to 2*pi. The differential for area would
> be
>>
>> dA = 1^2*cos(phi)*dphi*dtheta,
>>
>> and as you see, dA becomes vanishingly small at the two poles
> where
>> cos(phi) approaches zero. This would cause a theoretically
>> infinite
>> density there, using straight 'rand' outputs for phi.
>>
>> What you need to generate a uniform area distribution using
>> 'rand' is:
>>
>> theta = 2*pi*rand(1,n);
>> phi = asin(2*rand(1,n)-1);
>>
>> In N dimensions life gets increasingly difficult. The method
>> mentioned
>> by Peter using 'randn' is the one I always use.
>>
>> Roger Stafford
>>
> Right about that you are. The thing with randn, must be that you
> have to define both an inner and an outer radii to get the points
> on
> the surface?
>
> xyz=randn(3,1000);
> rad2=xyz(1,:).^2+xyz(2,:).^2+xyz(3,:).^2;
> r1=0.99^2;
> r2=1.01^2;
> ix=find(and(rad2>=r1,rad2<=r2));
>
> if r1->r2 then you need infinte number of points! Or have I
> missunderstood you?
>
> /Per
  
Oh, I'm stupid! It´s only to scale them all to the unit sphere!
x=xyz(1,:)./sqrt(rad2);

etc, but except for rad2=0, using find.

Nice idea,
Per

Subject: Generating random numbers on unit sphere

From: aklassyguy@hotmail.com

Date: 1 Jul, 2007 08:00:02

Message: 11 of 31

On Tue, 26 Jun 2007 15:48:25 -0000, kanika <kanika.dhyani@gmail.com>
wrote:

>Hi everyone,
>
>I am trying to generate random numbers on a unit sphere and further in
>higher dimensions also. Any suggestions to generate these
>efficiently? I know that numerical recipes should have a code for
>this in C, which I cannot find. Any suggestions or pointers or
>mathematical tips for this?
>
>Thank you
>Kanika

Another method is the "throw-away" method. Simply generate uniformly
distributed points in the interior of a hyper-cube (easy to do with a
uniform random number generator), throw away those points outside the
radius of the hyper-sphere, and normalize the rest. This gives you
uniformly distributed points on the surface of the hyper-sphere. It
costs you some extra uniform random number generator calls because of
the throw aways and a sqrt for the normalization, but you stay away
from trig and exp functions.

James Tursa

Subject: Generating random numbers on unit sphere

From: John D'Errico

Date: 1 Jul, 2007 04:50:19

Message: 12 of 31

aklassyguy wrote:
>
>
> On Tue, 26 Jun 2007 15:48:25 -0000, kanika
> <kanika.dhyani@gmail.com>
> wrote:
>
>>Hi everyone,
>>
>>I am trying to generate random numbers on a unit sphere and
> further in
>>higher dimensions also. Any suggestions to generate these
>>efficiently? I know that numerical recipes should have a code
for
>>this in C, which I cannot find. Any suggestions or pointers or
>>mathematical tips for this?
>>
>>Thank you
>>Kanika
>
> Another method is the "throw-away" method. Simply generate
> uniformly
> distributed points in the interior of a hyper-cube (easy to do with
> a
> uniform random number generator), throw away those points outside
> the
> radius of the hyper-sphere, and normalize the rest. This gives you
> uniformly distributed points on the surface of the hyper-sphere. It
> costs you some extra uniform random number generator calls because
> of
> the throw aways and a sqrt for the normalization, but you stay away
> from trig and exp functions.
>
> James Tursa

As a general rule, rejection techniques
are a good way to generate samples for
the most difficult domains.

You do need to be aware of the rejection
rate however, as this can severely impact
the speed of your code. For example, in
two dimensions, you will reject only a
relatively small fraction.

  1 - pi/4 = 0.2146...

In n dimensions, the fraction rejected
will grow large relatively fast.

  n = (2:10)';
  1 - pi.^(n/2)./gamma(1+n/2)./2.^n

ans =
       0.2146
       0.4764
      0.69157
      0.83551
      0.91925
      0.96309
      0.98415
      0.99356
      0.99751

Thus in 10 dimensions, expect to reject
fully 99.75% of the samples that you
generate. This means that you will need
to oversample your requirements by a
factor of roughly 400 to 1. In 100
dimensions the oversampling requirement
becomes so onerous that you will need to
generate roughly 10^70 samples merely
to gain one success. I think it likely
that you would need to buy a faster
computer than your current one to
generate that large of a random sample.

So when generating points on the surface
of a spherical domain in higher dimensions,
I'd strongly recommend use of the randn
solution.

John

Subject: Generating random numbers on unit sphere

From: ellieandrogerxyzzy@mindspring.com.invalid (Roger Stafford)

Date: 1 Jul, 2007 09:32:31

Message: 13 of 31

In article <qgne83dt3ngurbu5svhsrr3bu5v7hook81@4ax.com>,
aklassyguy@hotmail.com wrote:

> On Tue, 26 Jun 2007 15:48:25 -0000, kanika <kanika.dhyani@gmail.com>
> wrote:
>
> >Hi everyone,
> >
> >I am trying to generate random numbers on a unit sphere and further in
> >higher dimensions also. Any suggestions to generate these
> >efficiently? I know that numerical recipes should have a code for
> >this in C, which I cannot find. Any suggestions or pointers or
> >mathematical tips for this?
> >
> >Thank you
> >Kanika
>
> Another method is the "throw-away" method. Simply generate uniformly
> distributed points in the interior of a hyper-cube (easy to do with a
> uniform random number generator), throw away those points outside the
> radius of the hyper-sphere, and normalize the rest. This gives you
> uniformly distributed points on the surface of the hyper-sphere. It
> costs you some extra uniform random number generator calls because of
> the throw aways and a sqrt for the normalization, but you stay away
> from trig and exp functions.
>
> James Tursa
---------------------
  James, as the number of dimensions gets larger, the fraction of a
hypercube's volume that is occupied by an inscribed hypersphere becomes
smaller and smaller. For n equal to 10, the fraction becomes seriously
small for your rejection scheme to be practical. Only one in every 400
points randomly selected in a ten-dimensional hypercube will fall within
the hypersphere. With n equal to 20, the ratio becomes one in every forty
million! That's a pretty poor return on all the computation involved.

  To be precise, a hypercube with unit side equal to one has an
n-dimensional volume of one. A contained hypersphere within it has a
radius, R, of 1/2, and the formula for the volume of an n-dimensional
hypersphere of radius R is

 Vn(R) = 2*pi^(n/2)*R^n/(n*gamma(n/2))

See <http://people.csail.mit.edu/jrennie/writing/sphereVolume.pdf>

  Using matlab,

 R = 1/2;
 n = 10; V10 = 2*pi^(n/2)*R^n/(n*gamma(n/2))
 n = 20; V20 = 2*pi^(n/2)*R^n/(n*gamma(n/2))
 [V10;V20]

ans =

   0.00249039457019
   0.00000002461137

  Anyway, what is so hard about using the 'randn' function? Normalization
of n randn numbers to define a point on the hypersphere's surface is very
easy to do and there are no rejections at all.

Roger Stafford

Subject: Generating random numbers on unit sphere

From: James Tursa

Date: 1 Jul, 2007 19:51:10

Message: 14 of 31

On Sun, 1 Jul 2007 04:50:19 -0400, "John D'Errico"
<woodchips@rochester.rr.com> wrote:

>aklassyguy wrote:
>>
>>
>> On Tue, 26 Jun 2007 15:48:25 -0000, kanika
>> <kanika.dhyani@gmail.com>
>> wrote:
>>
>>>Hi everyone,
>>>
>>>I am trying to generate random numbers on a unit sphere and
>> further in
>>>higher dimensions also. Any suggestions to generate these
>>>efficiently? I know that numerical recipes should have a code
>for
>>>this in C, which I cannot find. Any suggestions or pointers or
>>>mathematical tips for this?
>>>
>>>Thank you
>>>Kanika
>>
>> Another method is the "throw-away" method. Simply generate
>> uniformly
>> distributed points in the interior of a hyper-cube (easy to do with
>> a
>> uniform random number generator), throw away those points outside
>> the
>> radius of the hyper-sphere, and normalize the rest. This gives you
>> uniformly distributed points on the surface of the hyper-sphere. It
>> costs you some extra uniform random number generator calls because
>> of
>> the throw aways and a sqrt for the normalization, but you stay away
>> from trig and exp functions.
>>
>> James Tursa
>
>As a general rule, rejection techniques
>are a good way to generate samples for
>the most difficult domains.
>
>You do need to be aware of the rejection
>rate however, as this can severely impact
>the speed of your code. For example, in
>two dimensions, you will reject only a
>relatively small fraction.
>
> 1 - pi/4 = 0.2146...
>
>In n dimensions, the fraction rejected
>will grow large relatively fast.
>
> n = (2:10)';
> 1 - pi.^(n/2)./gamma(1+n/2)./2.^n
>
>ans =
> 0.2146
> 0.4764
> 0.69157
> 0.83551
> 0.91925
> 0.96309
> 0.98415
> 0.99356
> 0.99751
>
>Thus in 10 dimensions, expect to reject
>fully 99.75% of the samples that you
>generate. This means that you will need
>to oversample your requirements by a
>factor of roughly 400 to 1. In 100
>dimensions the oversampling requirement
>becomes so onerous that you will need to
>generate roughly 10^70 samples merely
>to gain one success. I think it likely
>that you would need to buy a faster
>computer than your current one to
>generate that large of a random sample.
>
>So when generating points on the surface
>of a spherical domain in higher dimensions,
>I'd strongly recommend use of the randn
>solution.
>
>John

Oops.(Guess I should have done that calculation before posting :)

Subject: Generating random numbers on unit sphere

From: Jeffrey Erlich

Date: 1 Jul, 2007 16:09:52

Message: 15 of 31

This may be naive, but instead of generating n random numbers and
rejecting the ones that don't fall on the sphere, wouldn't it be much
more efficient to generate n-1 random number and then calculate the
last number.
x=rand*2-1;
y=rand*2-1;
z=1-randsign*sqrt(x^2+y^2);

function y=randsign
  y=1;
 if rand>0.5
y=-1;
end

Subject: Generating random numbers on unit sphere

From: ellieandrogerxyzzy@mindspring.com.invalid (Roger Stafford)

Date: 1 Jul, 2007 22:24:59

Message: 16 of 31

In article <ef5bc41.13@webcrossing.raydaftYaTP>, "Jeffrey Erlich"
<jerlich@NOprincetonSPAM.edu> wrote:

> This may be naive, but instead of generating n random numbers and
> rejecting the ones that don't fall on the sphere, wouldn't it be much
> more efficient to generate n-1 random number and then calculate the
> last number.
> x=rand*2-1;
> y=rand*2-1;
> z=1-randsign*sqrt(x^2+y^2);
>
> function y=randsign
> y=1;
> if rand>0.5
> y=-1;
> end
----------------
  Among other difficulties, this surface isn't a sphere, not even in three
dimensions. It's a conical surface, with equation x^2+y^2-(z-1)^2 = 0.

  There are ways of solving the n-dimensional form of this problem using
only n-1 random numbers but they aren't this simple.

Roger Stafford

Subject: Generating random numbers on unit sphere

From: John D'Errico

Date: 1 Jul, 2007 22:04:08

Message: 17 of 31

Jeffrey Erlich wrote:
>
>
> This may be naive, but instead of generating n random numbers and
> rejecting the ones that don't fall on the sphere, wouldn't it be
> much
> more efficient to generate n-1 random number and then calculate the
> last number.
> x=rand*2-1;
> y=rand*2-1;
> z=1-randsign*sqrt(x^2+y^2);
>
> function y=randsign
> y=1;
> if rand>0.5
> y=-1;
> end

Roger mentioned an important flaw.
There are others, as he alluded to.

One other is its not a uniform sampling
on the conic surface that it does
sample.

John

Subject: Generating random numbers on unit sphere

From: ellieandrogerxyzzy@mindspring.com.invalid (Roger Stafford)

Date: 2 Jul, 2007 03:54:20

Message: 18 of 31

In article <ef5bc41.15@webcrossing.raydaftYaTP>, "John D'Errico"
<woodchips@rochester.rr.com> wrote:

> Jeffrey Erlich wrote:
> >
> > This may be naive, but instead of generating n random numbers and
> > rejecting the ones that don't fall on the sphere, wouldn't it be
> > much
> > more efficient to generate n-1 random number and then calculate the
> > last number.
> > x=rand*2-1;
> > y=rand*2-1;
> > z=1-randsign*sqrt(x^2+y^2);
> >
> > function y=randsign
> > y=1;
> > if rand>0.5
> > y=-1;
> > end
>
> Roger mentioned an important flaw.
> There are others, as he alluded to.
>
> One other is its not a uniform sampling
> on the conic surface that it does
> sample.
>
> John
-------------------
  This surface of Jeffrey's is a peculiar one, consisting of two cone
halves meeting at the common vertex (0,0,1) and with vertical sides
defined by four vertical planes x = 1, x = -1, y = 1, and y = -1, each
bounded by two parabolas and two vertical lines. It would make an
interesting picture artistically, but I feel sure it is not the surface he
intended.

  However, strangely enough, if we ignore the fact that no random points
are being placed on the four vertical surfaces, those on the two conic
portions of the surface must, in fact, be uniformly distributed area-wise,
because of the constant angle, pi/4, the rays from the cone vertex make
with respect to the horizontal direction, and because the original
distribution of (x,y) points in the plane is uniform. There is a constant
differential relationship between areas, dAc = sqrt(2)*dAp, where Ac is
cone surface area and Ap is horizontal plane surface area. Of course this
would not be true if he had defined a true spherical surface as he
presumably intended.

Roger Stafford

Subject: Generating random numbers on unit sphere

From: John D'Errico

Date: 2 Jul, 2007 05:45:28

Message: 19 of 31

Roger Stafford wrote:
> -------------------
> This surface of Jeffrey's is a peculiar one, consisting of two
> cone
> halves meeting at the common vertex (0,0,1) and with vertical sides
> defined by four vertical planes x = 1, x = -1, y = 1, and y = -1,
> each
> bounded by two parabolas and two vertical lines. It would make an
> interesting picture artistically, but I feel sure it is not the
> surface he
> intended.
>
> However, strangely enough, if we ignore the fact that no random
> points
> are being placed on the four vertical surfaces, those on the two
> conic
> portions of the surface must, in fact, be uniformly distributed
> area-wise,
> because of the constant angle, pi/4, the rays from the cone vertex
> make
> with respect to the horizontal direction, and because the original
> distribution of (x,y) points in the plane is uniform. There is a
> constant
> differential relationship between areas, dAc = sqrt(2)*dAp, where
> Ac is
> cone surface area and Ap is horizontal plane surface area. Of
> course this
> would not be true if he had defined a true spherical surface as he
> presumably intended.
>
> Roger Stafford
>
  
Yes, I should have looked more carefully.
A plot would have sufficed. While its
clearly not spherical, the sampling will
indeed be uniform.

John

Subject: Generating random numbers on unit sphere

From: dave

Date: 2 Jul, 2007 05:57:50

Message: 20 of 31

What about setting a vector radius of 1 and calculating n-1 random
rotations, one around each axis? does that come out properly
distributed?

Subject: Generating random numbers on unit sphere

From: Per S

Date: 2 Jul, 2007 09:58:43

Message: 21 of 31

dave wrote:
>
>
> What about setting a vector radius of 1 and calculating n-1 random
> rotations, one around each axis? does that come out properly
> distributed?
  
You mean (0,0,1) rotated by Euler matrices? I think you need n
rotations not n-1. Only rotating by Rx and Ry is maby not sufficient,
but I am not sure? In 4d, you have to define the rotations also (I
remember a posting on this subject). Who are interested in n-d
spheres anyway (I'm just curious)? To avoid loops I think that you
will need to form 3x3 (in 3D) blocks in diagonal of a sparse matrix
and then do the rotation (matrix multiplication). Try it and do some
statistics to se if its true or not (a 2d histogram in theta phi
angles)!
/Per

Subject: Generating random numbers on unit sphere

From: dave

Date: 2 Jul, 2007 10:22:00

Message: 22 of 31

Per S wrote:
>
>
> dave wrote:
>>
>>
>> What about setting a vector radius of 1 and calculating n-1
> random
>> rotations, one around each axis? does that come out properly
>> distributed?
>
> You mean (0,0,1) rotated by Euler matrices? I think you need n
> rotations not n-1. Only rotating by Rx and Ry is maby not
> sufficient,
> but I am not sure? In 4d, you have to define the rotations also (I
> remember a posting on this subject). Who are interested in n-d
> spheres anyway (I'm just curious)? To avoid loops I think that you
> will need to form 3x3 (in 3D) blocks in diagonal of a sparse matrix
> and then do the rotation (matrix multiplication). Try it and do
> some
> statistics to se if its true or not (a 2d histogram in theta phi
> angles)!
> /Per
  
i'm more of an idea rat... but i expect it would be n-1. in 2d there
is only 1 angle, 3d only 2 angles, and thats as far as i usually work
in.

I am not big enough on statistics to be sure how to tell if they are
properly distributed on a sphere. i can see it on a circle because
you can unroll it into a line and see the distribution easily enough.
 but on a sphere if you try to map to a square you get distortion
that would make it not look right. also, what is the metric in 3d?
number of points per unit surface area?

Subject: Generating random numbers on unit sphere

From: Jeffrey Erlich

Date: 2 Jul, 2007 21:06:35

Message: 23 of 31

As several people have pointed out, the equation I described was not
a sphere.
i should have written

n=1000;
x=2*rand(n,1)-1; y=2*rand(n,1)-1;
z=(round(rand(n,1))*2-1).*sqrt(1-x.^2-y.^2);

gd=find(real(z)~=0); % reject cases where x^2+y^2>1
plot3(x(gd), y(gd), z(gd));

This produces points on a sphere, but as it was pointed out they are
not uniform in z, only in x and y.

Sorry for the careless typing.

Subject: Generating random numbers on unit sphere

From: Bernie Schattka

Date: 31 Jul, 2007 05:28:02

Message: 24 of 31

 kanika <kanika.dhyani@gmail.com> wrote in message <1182872905.154561.195950@o61g2000hsh.googlegroups.com>...
> Hi everyone,
>
> I am trying to generate random numbers on a unit sphere and further in
> higher dimensions also. Any suggestions to generate these
> efficiently? I know that numerical recipes should have a code for
> this in C, which I cannot find. Any suggestions or pointers or
> mathematical tips for this?
>
> Thank you
> Kanika
>
Try this:
  
  z = 2*rand(1,1000);
  rho = sqrt(1 - z.^2);
  phi = pi*(2*rand(1,1000) - 1);
  x = rho .* cos(phi);
  y = rho .* sin(phi);

The uniform random numbers can be visually verified by:
  plot3(x,y,z, '.');

Have fun

Cheers


   -bernie

Subject: Generating random numbers on unit sphere

From: Bernie Schattka

Date: 31 Jul, 2007 05:28:25

Message: 25 of 31

 kanika <kanika.dhyani@gmail.com> wrote in message <1182872905.154561.195950@o61g2000hsh.googlegroups.com>...
> Hi everyone,
>
> I am trying to generate random numbers on a unit sphere and further in
> higher dimensions also. Any suggestions to generate these
> efficiently? I know that numerical recipes should have a code for
> this in C, which I cannot find. Any suggestions or pointers or
> mathematical tips for this?
>
> Thank you
> Kanika
>
Try this:
  
  z = 2*rand(1,1000);
  rho = sqrt(1 - z.^2);
  phi = pi*(2*rand(1,1000) - 1);
  x = rho .* cos(phi);
  y = rho .* sin(phi);

The uniform random numbers can be visually verified by:
  plot3(x,y,z, '.');

Have fun

Cheers


   -bernie

Subject: Generating random numbers on unit sphere

From: Tom Bryan

Date: 31 Jul, 2007 08:16:33

Message: 26 of 31

kanika wrote:
> Hi everyone,
>
> I am trying to generate random numbers on a unit sphere and further in
> higher dimensions also. Any suggestions to generate these
> efficiently? I know that numerical recipes should have a code for
> this in C, which I cannot find. Any suggestions or pointers or
> mathematical tips for this?
>
> Thank you
> Kanika
>

% From G.S. Watson, Statistics on Spheres, Wiley, New York, 1983:
%
% If the elements of x = [x1 x2 ... xn] are real-valued independent
% Gaussian random variable with mean 0 and variance 1, then the
% random variable defined by z = x/norm(x,2) is distributed
% uniformly and randomly on the unit sphere S_{n-1} in n dimensions.
%
% Here's the MATLAB

% Generate NPOINTS on a sphere in DIMS dimensions:
dims=3;
npoints=1000;
x=randn(dims,npoints);
z=zeros(dims,npoints);
for k=1:npoints
   z(:,k)=x(:,k)/norm(x(:,k),2);
end

% Plot for 2D or 3D
switch(dims)
case 2
   plot(z(1,:),z(2,:),'.')
   axis equal
   figure(gcf)
  case 3
   plot3(z(1,:),z(2,:),z(3,:),'.')
   axis equal
   figure(gcf)
end

Subject: Generating random numbers on unit sphere

From: John D'Errico

Date: 31 Jul, 2007 12:25:26

Message: 27 of 31

"Bernie Schattka" <bernie.schattka@mathworks.com> wrote in message <f8mh9p$ag2$1@fred.mathworks.com>...
> Try this:
>
> z = 2*rand(1,1000);
> rho = sqrt(1 - z.^2);
> phi = pi*(2*rand(1,1000) - 1);
> x = rho .* cos(phi);
> y = rho .* sin(phi);
>
> The uniform random numbers can be visually verified by:
> plot3(x,y,z, '.');

While this will do fine in 3 dimensions,
sampling in higher dimensions will take
more effort to write. Can you perhaps
give an example of the analogous code
to sample on the surface of a 5 or 12
dimensional sphere?

This is why the randn solution is a
better choice.

John

Subject: Generating random numbers on unit sphere

From: roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson)

Date: 31 Jul, 2007 17:22:33

Message: 28 of 31

In article <f8n9nl$3vc$1@fred.mathworks.com>,
John D'Errico <woodchips@rochester.rr.com> wrote:
>"Bernie Schattka" <bernie.schattka@mathworks.com> wrote in message <f8mh9p$ag2$1@fred.mathworks.com>...
>> Try this:

>> z = 2*rand(1,1000);
>> rho = sqrt(1 - z.^2);
>> phi = pi*(2*rand(1,1000) - 1);
>> x = rho .* cos(phi);
>> y = rho .* sin(phi);

>While this will do fine in 3 dimensions,
>sampling in higher dimensions will take
>more effort to write.


[Note: Bernie Schattka is one of my co-workers, a hall over from me;
we discussed this earlier today.]

First off, Bernie had a typo; the first line should be

z = 2*rand(1,1000) - 1;

Secondly, I do not believe that this will produce a uniform random
distribution over the sphere: I believe that this will produce
points that are denser towards z = +/- 1.

Suppose that the random number generator generates two z with
the same value. The rho for those will be the same, and will
be the radius of the circle which is the sphere projected on to
the xy plane at that z. The random phi is the angle around
that circle, and the x and y calculations are standard polar to
cartesian transformations of polar(phi,rho). Hence, for
those two duplicate z, the maximum distance between the generated
points would be 2 * rho (i.e., 180 degrees apart, the diameter
of the circle), and the average distance between the two generated
points would be rho. But when these duplicate z values are near
z = +/- 1 then rho is relatively small, and when the duplicate
z values are near z = 0, then rho is relatively large. Therefore
the average distance between generated points with duplicated z
values is small near the z = +/- 1 pole, and bigger near z = 0,
so the distribution of points will not be uniform random.

Now, exact duplicate random numbers are possible but not
especially high probability, but the argument can be extended
for z values that are nearly the same, using a 3D distance
function rather than the 2D function I used above: z values
that are nearly the same and are close to +/- 1 will result
in generated points closer together than z values that are nearly
the same and are close to 0. And since the order of generation
is irrelevant, you could sort(z) and observe the distance variation
fairly directly.


function [x,y,z] = rsphere(N)
z = 2 * rand(1,N) - 1;
rho = sqrt(1-z.^2);
phi = pi * (2 * rand(1,N) - 1);
x = rho .* cos(phi);
y = rho .* sin(phi);


Then

[x,y,z] = rsphere(100000);
T = sortrows([x;y;z]',3);
d = sqrt(sum((T(2:end,:,:) - T(1:end-1,:,:)).^2,2));
plot(d);

If you use only 1000 points, as in Bernie's posting, then
you get a fairly high variance but you can see the general shape
of the distance distribution; by the time you get to 100000
points, the variance is pretty much invisible.
--
Programming is what happens while you're busy making other plans.

Subject: Generating random numbers on unit sphere

From: ellieandrogerxyzzy@mindspring.com.invalid (Roger Stafford)

Date: 31 Jul, 2007 17:34:37

Message: 29 of 31

In article <f8mh9p$ag2$1@fred.mathworks.com>, "Bernie Schattka"
<bernie.schattka@mathworks.com> wrote:

> Try this:
>
> z = 2*rand(1,1000);
> rho = sqrt(1 - z.^2);
> phi = pi*(2*rand(1,1000) - 1);
> x = rho .* cos(phi);
> y = rho .* sin(phi);
>
> The uniform random numbers can be visually verified by:
> plot3(x,y,z, '.');
>
> Have fun
> Cheers
> -bernie
------------------
  Bernie, you had better change your first line to read:

   z = 2*rand(1,1000)-1;

Roger Stafford

Subject: Generating random numbers on unit sphere

From: ellieandrogerxyzzy@mindspring.com.invalid (Roger Stafford)

Date: 31 Jul, 2007 18:34:45

Message: 30 of 31

In article <f8nr4p$gcf$1@canopus.cc.umanitoba.ca>,
roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote:

> In article <f8n9nl$3vc$1@fred.mathworks.com>,
> John D'Errico <woodchips@rochester.rr.com> wrote:
> >"Bernie Schattka" <bernie.schattka@mathworks.com> wrote in message
<f8mh9p$ag2$1@fred.mathworks.com>...
> >> Try this:
>
> >> z = 2*rand(1,1000);
> >> rho = sqrt(1 - z.^2);
> >> phi = pi*(2*rand(1,1000) - 1);
> >> x = rho .* cos(phi);
> >> y = rho .* sin(phi);
>
> >While this will do fine in 3 dimensions,
> >sampling in higher dimensions will take
> >more effort to write.
>
>
> [Note: Bernie Schattka is one of my co-workers, a hall over from me;
> we discussed this earlier today.]
>
> First off, Bernie had a typo; the first line should be
>
> z = 2*rand(1,1000) - 1;
>
> Secondly, I do not believe that this will produce a uniform random
> distribution over the sphere: I believe that this will produce
> points that are denser towards z = +/- 1.
>
> Suppose that the random number generator generates two z with
> the same value. The rho for those will be the same, and will
> be the radius of the circle which is the sphere projected on to
> the xy plane at that z. The random phi is the angle around
> that circle, and the x and y calculations are standard polar to
> cartesian transformations of polar(phi,rho). Hence, for
> those two duplicate z, the maximum distance between the generated
> points would be 2 * rho (i.e., 180 degrees apart, the diameter
> of the circle), and the average distance between the two generated
> points would be rho. But when these duplicate z values are near
> z = +/- 1 then rho is relatively small, and when the duplicate
> z values are near z = 0, then rho is relatively large. Therefore
> the average distance between generated points with duplicated z
> values is small near the z = +/- 1 pole, and bigger near z = 0,
> so the distribution of points will not be uniform random.
>
> Now, exact duplicate random numbers are possible but not
> especially high probability, but the argument can be extended
> for z values that are nearly the same, using a 3D distance
> function rather than the 2D function I used above: z values
> that are nearly the same and are close to +/- 1 will result
> in generated points closer together than z values that are nearly
> the same and are close to 0. And since the order of generation
> is irrelevant, you could sort(z) and observe the distance variation
> fairly directly.
>
>
> function [x,y,z] = rsphere(N)
> z = 2 * rand(1,N) - 1;
> rho = sqrt(1-z.^2);
> phi = pi * (2 * rand(1,N) - 1);
> x = rho .* cos(phi);
> y = rho .* sin(phi);
>
>
> Then
>
> [x,y,z] = rsphere(100000);
> T = sortrows([x;y;z]',3);
> d = sqrt(sum((T(2:end,:,:) - T(1:end-1,:,:)).^2,2));
> plot(d);
>
> If you use only 1000 points, as in Bernie's posting, then
> you get a fairly high variance but you can see the general shape
> of the distance distribution; by the time you get to 100000
> points, the variance is pretty much invisible.
-------------------------
  Walter, I think you are mistaken. Provided Bernie makes the indicated
correction on his first line, his code will produce a uniform distribution
over the surface of the unit sphere. It will NOT be more dense near the
poles.

  Your argument about average distances between points with the same z
value is faulty. The criterion for uniform distribution of points within
a two-dimensional surface is that the expected number of points falling
within any region of the surface, no matter what size or shape, be
proportional to that region's area. With Bernie's cylindrical coordinates
we have for an infinitesimal area dA corresponding to infinitesimal
changes dz and dphi, regarding rho as a function of z,

 dA = 1*dphi*rho*sqrt(dz^2+drho^2) =
      dphi*rho*sqrt(1+(drho/dz)^2)*dz =
      dphi*rho*sqrt(1 + z^2/rho^2)*dz =
      dphi*rho*sqrt(1/rho^2)*dz = dphi*dz.

In other words, the Jacobian with respect to z and phi is a constant unit
value. Since Bernie has provided each of z and phi with independent,
uniformly distributed 'rand' values, then it follows that the surface will
be uniformly distributed, areawise.

Roger Stafford

Subject: Generating random numbers on unit sphere

From: Richard Willey

Date: 31 Jul, 2007 15:17:04

Message: 31 of 31

For what its worth, Cleve Moler's textbook on numerical computing has an
entire chapter on random number generators.
Several of the topics that he discusses touch (tangentially) on the question
at hand.
He also has a number of M files that address similar topics.

The text is available online at http://www.mathworks.com/moler/random.pdf

Richard Willey
rwilley@mathworks.com

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us