3
$\begingroup$

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.

Lune

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 plotenter image description here

1 Answers 1

1

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:

results

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]
  • 0
    Sorry I cant locate it, could you edit my question?2011-06-19
  • 0
    @yayu: Rather than editing it into your question, I put the edited code into my answer.2011-06-19
  • 0
    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?2011-06-19
  • 0
    My mistake was to equate the random variables to the CDF while I equated the PDF's.. that's why I am missing the 12011-06-19
  • 0
    @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).2011-06-19
  • 0
    @Isaac thanks, here's my [answer](http://math.stackexchange.com/questions/44689/how-to-find-a-random-axis-or-unit-vector-in-3d/46234#46234) I've attributed you and made it CW2011-06-19
  • 0
    @yayu: I actually hadn't thought much about this in terms of the original question. I suspect, from the looks of the 3d plot, that what you're generating is skewed toward the poles (or some poles), as the original question suggests can happen, rather than being uniform.2011-06-19
  • 0
    If the method itself is wrong. hmm.. I would have to sleep on this one.2011-06-19
  • 0
    @yayu: I'm not exactly an expert on this kind of probability—once probability gets to PDFs and CDFs, my knowledge and understanding fades quickly; my answer here was based entirely on debugging the Mathematica code and my suggestion about the uniformity is purely based on how the plot looks.2011-06-19