Uniform Distribution on a Sphere: Part 3b (The Solution by Calculus)

(This post continues where we left off in Part 2 and gives a more direct solution than Part 3a. However this post uses multivariable/vector calculus, and I do not provide many explanations for readers unfamiliar with that material.)

As we saw in Part 2, the problem with picking our spherical angles \theta and \phi uniformly was that the area on our two-dimensional grid changed when we transformed that surface into the sphere in three dimensions. In particular if we wanted to calculate the probability mass in a region R of our spherical angles, we have the formula from Part 0:

\int_R\!f_{\theta,\phi}(\theta,\phi)\,\mathrm{d}\theta\,\mathrm{d}\phi

where f_{\theta,\phi} is our probability density function. However when we move to the sphere our area elements are no longer little rectangles with area \mathrm{d}\theta\cdot\mathrm{d}\phi; they are little patches on the surface of a sphere. What area do these corresponding patches have?

Spherical Coordinates
Figure 1: The spherical coordinate system

In Part 1 we defined spherical coordinates graphically using the scheme shown in Figure 1. With some basic trigonometry we can describe the transformation from spherical to Euclidean coordinates by:

\begin{align*}x&=r\sin\phi\cos\theta \\ y&=r\sin\phi\sin\theta \\ z&=r\cos\phi.\end{align*}

And so we can parameterize the surface of the unit sphere as a vector valued function

 \mathbf{\vec{x}}(\theta,\phi) = \left( \begin{array}{c} \sin\phi\cos\theta \\ \sin\phi\sin\theta \\ \cos\phi\end{array} \right).

Now we consider our two-dimensional space of spherical angles \theta and \phi that we can think of as a rectangle where the coordinates of each point specify the longitude and co-latitude. Take any region R in this space. The formula for surface integration gives us that

 \mathrm{Area}_{\mathbf{\vec{x}}(R)}=\int_{\mathbf{\vec{x}}(R)} \mathrm{d}\mathbf{\vec{x}(R)} = \iint_R \left\| \frac{\partial\mathbf{\vec{x}}}{\partial\theta}\times\frac{\partial\mathbf{\vec{x}}}{\partial\phi} \right\|\,\mathrm{d}\theta\,\mathrm{d}\phi,

where \mathrm{Area}_{\mathbf{\vec{x}}(R)} is the area of the region R after it has been transformed by the function \mathbf{\vec{x}}(\theta,\phi) onto the surface of the sphere in three-dimensional space, \mathrm{d}\mathbf{\vec{x}}(R) represents the area of the area of a little patch on the surface of the sphere and \times denotes the vector cross product.[1]

We won’t perform all the calculations here, but you can differentiate and use the rule \sin^2\alpha+\cos^2\alpha=1 to get that

 \left\| \frac{\partial\mathbf{\vec{x}}}{\partial\theta}\times\frac{\partial\mathbf{\vec{x}}}{\partial\phi} \right\| = |\sin\phi| = \sin\phi,

since \phi\in[0,\pi].

Therefore we get that the small patches on the surface of the sphere at a given longitude and co-latitude (\theta,\phi) have area \sin\phi\cdot\mathrm{d}\theta\cdot\mathrm{d}\phi. This fits with what we saw in Part 2. At the equator (\phi=\pi/2) we have that the area element is \sin\frac{\pi}{2}\cdot\mathrm{d}\theta\cdot\mathrm{d}\phi=\mathrm{d}\theta\cdot\mathrm{d}\phi. That means that at the equator, the area element on the surface of the sphere corresponds exactly with the area element on our equirectangular map. As we move away from the equator, we see that \sin\phi decreases and so the equally sized rectangles on our map will be transformed into increasingly smaller patches on the sphere. Once we reach the poles the area element vanishes.

Now that we have figured out exactly how area on our rectangular grid corresponds to area on the surface of the sphere, how do we solve the original problem? Recall from Part 0 that a probability density function defined on a sample space A must:

  1. be non-negative at every point in A
  2. integrate to 1 over A.

Furthermore the uniform distribution on A is such that the probability distribution is constant on A. Therefore for a uniform distribution on the unit sphere S, we must have the probability density function f_S=c on the sphere where c=\left(\int_S\mathrm{d}S\right)^{-1}=1/(4\pi), a constant. The formula for integrating a function over a surface gives us that

 \int_{\mathbf{\vec{x}}(R)} \!c\,\mathrm{d}\mathbf{\vec{x}(R)} = \iint_R \!c\sin\phi\,\mathrm{d}\theta\,\mathrm{d}\phi.

We can see that the integral on the left gives the mass contained by the uniform distribution on the sphere over the region \mathbf{\vec{x}}(R). The quantity on the right then must be the mass in the region R on our map of \theta and \phi under the density f(\theta,\phi)=c\sin\phi.[2] Therefore if we could draw (\theta,\phi) from this density function, we would be able to generate uniformly random points on the sphere.

Generating random samples from a given density is not a trivial problem. We will use a trick that will allow us to do this as long as we have a method to generate uniformly random samples on an interval (we do). If we go back to our surface integration formula, we can perform a change of variables to get

\int_{\mathbf{\vec{x}}(R)} \!c\,\mathrm{d}\mathbf{\vec{x}(R)} = \iint_R \!c\sin\phi\,\mathrm{d}\theta\,\mathrm{d}\phi = \iint_{R^\prime} \!c\,\mathrm{d}\theta\,\mathrm{d}(-\cos\phi).

In other words, we see that drawing from our desired density is equivalent to drawing samples uniformly in \theta and -\cos\phi. Since -\cos\phi takes all values in the interval [-1,1] for \phi\in[0,\pi], we merely need to draw samples uniformly in the rectangle [0,2\pi]\times[-1,1]. This is exactly as we did in Part 3a where that rectangle was the Lambert cylindrical map.

While we have solved the original problem, this solution is cumbersome and extension to higher dimensions requires extra work. Continue on to Part 4 to see a different approach to the problem that extends well to any dimension and is computationally appealing.

  1. [1] Note that \left\| \frac{\partial\mathbf{\vec{x}}}{\partial\theta}\times\frac{\partial\mathbf{\vec{x}}}{\partial\phi} \right\| is equivalent to the determinant of the Jacobian of the spherical coordinates transformation evaluated at r=1. I chose to show the solution to this problem directly using surface integrals, but the Jacobian is the standard way to handle changes of coordinates in integrals.
  2. [2] How do we know this is actually a density function? It is non-negative as \sin\phi>0 for \phi\in[0,\pi]. It integrates to 1 over the entire region of spherical angles by the surface integral formula.