Sunday, March 30, 2008

Uniform deviates on the surface of a sphere

Here is a simple example of how to generate points randomly on the surface of the sphere that are uniformly distributed. The reason I'm posting it here is that it shows how convenient mayavi's mlab is for this sort of thing.

Anyone who has studied pseudo random numbers should know that the following will not generate points uniformly on the sphere. Here is a demonstration of the fact:

$ ipython -wthread
In [1]: from numpy import *
In [2]: from enthought.mayavi import mlab
In [3]: p = random.random(10000)*2*pi
In [4]: t = random.random(10000)*pi
In [5]: x, y, z = sin(t)*cos(p), sin(t)*sin(p), cos(t)
In [6]: g = mlab.points3d(x, y, z, z, mode='point')
In [7]: g.actor.property.point_size = 2

The point size is made larger so you can see the points more clearly. Here is the picture it produces:


Notice the clustering at the poles. OTOH, the following works nicely,

In [3]: p = random.random(10000)*2*pi
In [10]: mlab.clf()
In [12]: t = arccos(random.random(10000)*2.0 - 1.0)
In [13]: x, y, z = sin(t)*cos(p), sin(t)*sin(p), cos(t)
In [14]: g1 = mlab.points3d(x, y, z, z, mode='point')
In [15]: g1.actor.property.point_size = 2


This is the picture this one produces:



As you can see the transformation has worked and generates what appears to be a uniform distribution on the surface of a sphere. Visual testing isn't good enough uniformness but that isn't the point here.

2 comments:

Nicolas said...

Hi,
another way is to sample Gaussian vectors of dimension d, with mean 0_d, and covariance matrix I_d, and to normalise them:
x = randn(n)
u = x/sqrt(sum(x))

Code is more compact, and looks faster, but I have tried to compare both approaches.
Nicolas CHOPIN, ENSAE

Nicolas said...

oops, sorry, code is not correct:
x = randn(d,n)
u = x/sqrt(sum(x**2,0))
if d==2:
plot(u[0,:],u[1,:])