I was inspired by this question to play around a little bit (its a weekend). I was pretty confident of my derivation and thought it might be nice to supplement it with a pretty picture. However, instead of picking random points from the entire sphere the result I get are a bunch of random points from two opposite lunes. I've been scratching my head for an hour, and the only suspicion is that the distributions I have derived must be wrong (i.e not the code, I think). I am posting here what I would have hoped to be an answer.
If $\Omega$ is the solid angle, then the probability that a point lies in an infinitesimal cone is $p(\Omega)d\Omega$ On normalizing, we get $\displaystyle \int_{\mbox{unit sphere}} p(\Omega)d\Omega = 1$ $p(\Omega) = \frac{1}{4\pi}$
We can also write $p(\Omega)d\Omega= p(\theta,\phi) d\theta d\phi$ Where $\theta$ is the polar angle and $\phi$ is the azimuthal angle. Using $\displaystyle d\Omega = \sin \theta d\theta d\phi$ we get $p(\theta,\phi) = \frac{\sin\theta}{4\pi}$
So $p(\theta) = \int_0^{2\pi}p(\theta,\phi)d\phi = \frac{\sin \theta}{2}$ $p(\phi) = \int_0^\pi p(\theta,\phi)\sin\theta d\theta = \frac{1}{2\pi}$.
To implement this we choose two random variables $u,v$ which are uniform on the interval $(0,1)$ $\theta = \sin^{-1}(2u)$ $\phi = 2\pi v$
For programming ease I map these points into the cartesian system with the standard transformations.
Mathematica Code:
theta := ArcSin[2*RandomReal[UniformDistribution[{0, 1}]]]; phi := 2 Pi RandomReal[UniformDistribution[{0, 1}]]; plotPts = Table[{phi, theta}, {5000}]; cartesian[{theta_, phi_}] := {Sin[theta]*Cos[phi], Sin[theta]*Sin[phi], Cos[theta]}; cPts = cartesian /@ plotPts; points = Graphics3D[Point /@ cPts]; Show[points]
Mistake: Random variables should be equated to the CDF (probability) not the PDF. On integrating the PDFs to get the CDFs, the correct random variables are $\theta = \cos^{-1}(2u-1)$ $\phi = 2\pi v$
And the resulting plot