Random points on the Earth

December 7, 2013
By

(This article was first published on Freakonometrics » R-english, and kindly contributed to R-bloggers)

The problem with puzzles is that you keep it in your head for days, until you find an answer. Or at least some ideas about a possible answer. This is what happened to me a few weeks ago, when a colleague of mine asked me the following question : Consider http://latex.codecogs.com/gif.latex?n points uniformly distributed on a sphere. What is the probability that the http://latex.codecogs.com/gif.latex?n points lie on a same hemisphere, for some hemisphere (there is no south or north here) ?

Analogously, what is the probability to see the http://latex.codecogs.com/gif.latex?n points on the Earth, at the same time, from somewhere in the galaxy ? (even extremely far away, so we can see a complete hemisphere) I wanted to use Monte Carlo simulations to estimate that probability, for some http://latex.codecogs.com/gif.latex?n. But it was difficult. I mean, given http://latex.codecogs.com/gif.latex?n points on the sphere, in  can you easily determine if they lie on a common hemisphere, or not ? I did try with distance, or angle, but I could not find a simple answer. So I tried a technique I did learn a few years ago : if you cannot do something in dimension 3, try first in dimension 2.

Again, I could not find a simple answer. But there is a simple technique.

  1. Draw http://latex.codecogs.com/gif.latex?n points on the unit sphere, which simply means generate http://latex.codecogs.com/gif.latex?n random variables http://latex.codecogs.com/gif.latex?%20U_1,\ldots%20,U_n\sim\mathcal{U}([0,2\pi])
  2. Try to find http://latex.codecogs.com/gif.latex?\theta\in[0,2\pi] such that, after a rotation (with angle http://latex.codecogs.com/gif.latex?\theta) all the points lie in the upper part (the North hemisphere, for instance)

So the question here is simply : is

http://latex.codecogs.com/gif.latex?\max_{\theta\in[0,2\pi]}\left\{\sum_{i=1}^n%20\boldsymbol{1}(\sin(U_i+\theta)\geq%200)\right\}

equal to http://latex.codecogs.com/gif.latex?n ? Of course, from a computational point of view, it is slightly more complex, since this function is not differentiable.

n=5
Theta=runif(n)*2*pi
top=Vectorize(function(theta) sum(sin(Theta+theta)>=0))

So a simple strategy can be to compute those values on a finite grid, and to check if, for some http://latex.codecogs.com/gif.latex?\theta, all the points lie in the upper part

max(top(seq(0,2*pi,length=6001)))==n

Hence, with the following sample, all the points cannot be on a common hemisphere

set.seed(2)
Theta=runif(5)*2*pi

while, for another sample of points, it was possible

set.seed(7)
Theta=runif(5)*2*pi

With this simple code, we get get the probability, but only in dimension 2 (so far)

SIM=Vectorize(function(n) simul(n,1000))
plot(3:10,SIM(3:10))

In dimension 3, it is still possible to use also a polar representation. Things are easier to generate, but also it is simple to consider rotations. And again, a simple algorithm can be derived,

But there is a simple technique.

  1. Draw http://latex.codecogs.com/gif.latex?n points on the unit sphere, which simply means generate http://latex.codecogs.com/gif.latex?2n random variables http://latex.codecogs.com/gif.latex?%20U_1,\ldots%20,U_n\sim\mathcal{U}([0,2\pi]) and http://latex.codecogs.com/gif.latex?V_1,\ldots,V_n\sim\mathcal{U}([-\pi/2,\pi/2])
  2. Try to find http://latex.codecogs.com/gif.latex?\theta and http://latex.codecogs.com/gif.latex?\varphi such that, after a rotation (with angles http://latex.codecogs.com/gif.latex?\theta and http://latex.codecogs.com/gif.latex?\varphi) all the points lie in some given part (say  http://latex.codecogs.com/gif.latex?\{x\geq0\})

As mentioned by Dominque, using this technique, points are not uniformly distributed on the sphere. Instead we can use

  1. Draw http://latex.codecogs.com/gif.latex?n points on points from a trivaraite Gaussian distribution, and normalize it
  2. Get the polar coordinates, and try to find http://latex.codecogs.com/gif.latex?\theta and http://latex.codecogs.com/gif.latex?\varphi such that, after a “rotation” (with angles http://latex.codecogs.com/gif.latex?\theta and http://latex.codecogs.com/gif.latex?\varphi) all the points lie in some given part (say  http://latex.codecogs.com/gif.latex?\{x\geq0\})

So the question here is simply : is

http://latex.codecogs.com/gif.latex?\max_{(\theta,\varphi)\in[0,2\pi]\times[-\pi/2,\pi/2]}\left\{\sum_{i=1}^n%20\boldsymbol{1}(\sin(U_i+\theta)\cos(V_i+\varphi)\geq%200)%20\right\}

(I am not sure about the set of angles for the rotations, so I tried a larger one, just in case). Again, it would be complex, or more complex than before, because we need here a joint grid. For instance

MZ=matrix(rnorm(n*3),n,3)
d=apply(MZ,1,function(z) sqrt(sum(z^2)))
X=MZ[,1]/d; Y=MZ[,2]/d; Z=MZ[,3]/d; 
Theta=acos(Z)
Phi=acos(X/sqrt(X^2+Y^2))*(Y>=0)+(2*pi-acos(X/sqrt(X^2+Y^2)))*(Y<0)
top=function(theta,phi) sum(sin(Theta+theta)*cos(Phi+phi)>=0)
TOP=mapply(top,rep(seq(0,2*pi,length=1001),1001),rep(seq(-pi,pi,length=1001),each=1001))
max(TOP)==n

As we can see with the red curve, there might be some problems here. Because, (as mention in kmath327), this problem was solved in any dimension in Wendel (1962) with the following simple expression : in dimension http://latex.codecogs.com/gif.latex?d, the probability that the http://latex.codecogs.com/gif.latex?n points (uniformly distributed on a sphere) lie on a same hemisphere is exactly

http://latex.codecogs.com/gif.latex?p(d,n)=%20\frac{1}{2^{n-1}}\sum_{i=0}^{d-1}%20\binom{n-1}{i}

p=function(d,n) .5^(n-1) * sum(choose(n-1,0:(d-1)))

Note that Leonard Savage proved (a few years before) that

http://latex.codecogs.com/gif.latex?p(d,d+1)=1%20-\frac{1}{2^d}

(which can be obtained easily actually). I do not see what’s wrong with my Monte Carlo simulations… and if anyone has a nice Monte Carlo strategy to get that probability, I’d be glad to hear it !

To leave a comment for the author, please follow the link and comment on his blog: Freakonometrics » R-english.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.