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