Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

I have posted a previous question, this is related but I think it is better to start another thread. This time, I am wondering how to generate uniformly distributed points inside the 3-d unit sphere and how to check the distribution visually and statistically too? I don't see the strategies posted there directly transferable to this situation.

Please help. Many thanks.

share|improve this question
3  
The techniques in the previous question apply directly once you observe that the number of points within distance $r$ of the origin must be proportional to $r^3$. Thus if you generate an independent uniform variate $u$ in $[0,1]$ along with a point $w$ on the surface of the sphere, scaling $w$ by $u^{1/3}$ does the trick. –  whuber Mar 8 '11 at 16:04
    
@whuber: maybe I just did not get the essence of the previous techniques. Let me try what you described. Additionally, what are the ways to check the uniformity here again? –  Qiang Li Mar 8 '11 at 17:10
1  
@Qiang Ripley's K function and chi-squared tests. You could also check the uniformity of the radial projection of the points on the sphere's surface, the uniformity of the cube of the lengths of the points, and the independence of those two. –  whuber Mar 8 '11 at 17:18
    
For me, it is not so obvious what "uniformly distributed" means... And probably a try to define it will automagically create a generating algorithm (= –  mbq Mar 8 '11 at 19:33
    
@mbq, I think to define the term, we need to have a p.d.f. of $f_{R, \Theta, \Phi}(r,\theta, \phi)=r^2$. –  Qiang Li Mar 8 '11 at 21:12

2 Answers 2

The easiest way is to sample points uniformly in the corresponding hypercube and discard those that do not lie within the sphere. In 3D, this should not happen that often, about 50% of the time. (Volume of the hypercube is 1, volume of the sphere is $\frac{4}{3}\pi r^3 = 0.523...$.)

share|improve this answer
    
+1. This is one of the techniques recommended by the comp.graphics.algorithms FAQ "Uniform random points on sphere". –  David Cary Mar 17 '11 at 2:51

You can also do this in spherical coordinates, in which case there is no rejection. First you generate the radius and the two angles at random, then you use the transition formula to recover $x$, $y$ and $z$ ($x = r \sin \theta \cos \phi$, $y = r \sin \theta \sin \phi$, $z = r \cos \theta$).

You generate $\phi$ unifomly between $0$ and $2\pi$. The radius $r$ and the inclination $\theta$ are not uniform though. The probability that a point is inside the ball of radius $r$ is $r^3$ so the probability density function of $r$ is $3 r^2$. You can easily check that the cubic root of a uniform variable has exactly the same distribution, so this is how you can generate $r$. The probability that a point lies within a spherical cone defined by inclination $\theta$ is $(1-\cos\theta)/2$ or $1 - (1-\cos (-\theta))/2$ if $\theta > \pi/2$. So the density $\theta$ is $sin(\theta)/2$. You can check that minus the arccosine of a uniform variable has the proper density.

Or more simply, we can simulate the cosine of $\theta$ uniformly beteen $-1$ and $1$.

In R this would look as shown below.

n <- 10000 # For example n = 10,000.
phi <- runif(n, max=2*pi)
r <- runif(n)^(1/3)
cos_theta <- runif(n, min=-1, max=1)
x <- r * sqrt(1-cos_theta^2) * cos(phi)
y <- r * sqrt(1-cos_theta^2) * sin(phi)
z <- r * cos_theta

In the course of writing and editing this answer, I realized that the solution is less trivial that I thought.

I think that the easiest and computationally most efficient method is to follow @whuber's method to generate $(x,y,z)$ on the unit sphere as shown on this post and scale them with $r$.

xyz <- matrix(rnorm(3*n), ncol=3)
lambda <- runif(n)^(1/3) / sqrt(rowSums(xyz^2))
xyz <- xyz*lambda
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.