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
    @yay$u$: I'm not exactl$y$ an expert on this kind of probabilit$y$—once probabilit$y$ 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