# Trying to pick a random point on sphere end up picking from a lune

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

-

I'm fairly certain that you're just missing a -1 inside the last ] of the first line. When I ran the code, I saw what you saw, but with a whole mess of errors about complex numbers where there shouldn't be any. Looking over your code, it seemed they must be coming from the ArcSin[], since it looks like you're taking the arcsine of numbers in the interval $[0,2]$ instead of $[-1,1]$. Hence, my thought about a -1, which for me seems to generate a more-covered sphere:

edit: Here's the code with the modification:

theta := ArcSin[2*RandomReal[UniformDistribution[{0, 1}]]-1 (* Here's the -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]

-
Sorry I cant locate it, could you edit my question? –  kuch nahi Jun 19 '11 at 5:23
@yayu: Rather than editing it into your question, I put the edited code into my answer. –  Isaac Jun 19 '11 at 5:25
I just realized what I did wrong. It is in the derivation... if I am able to correct it should I leave it as an edit in the question or post it as an answer? –  kuch nahi Jun 19 '11 at 5:26
My mistake was to equate the random variables to the CDF while I equated the PDF's.. that's why I am missing the 1 –  kuch nahi Jun 19 '11 at 5:27
@yayu: I'd say post it as an answer, since it's a resolution of the issue you were asking about (and you can mark it as accepted to indicate that your question is answered). –  Isaac Jun 19 '11 at 5:28