Random points on the Earth

December 7, 2013
By

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 Dominique, 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 !


Print This Post Print This Post

Tags: , ,

Posted in: Computer, Maths, Geometry

5 Responses to “ Random points on the Earth ”

  1. Eeshan on 11/12/2013 at 4:03 AM

    It’s important to explicitly state what is meant by “Consider n points uniformly distributed on a sphere.” This statement is not unambiguous for a number of reasons.

    1. One way to generate seemingly ‘uniform’ points is to uniformly sample from a the set of azimuthal angles and the polar angles. One (of many possibly) alternative way(s) is to uniformly sample the distance of the equator-parallel plane from the centre of the earth and combine with a uniformly generated longitude.

    These two generate different point densities because one varies linearly from the centre to the poles in a straight line, and the other varies linearly from the equator to the poles *along the surface of the earth*.

    This is not a fault of the analysis, but a missing ingredient in the posed question

    2. This is more of a nitpicky mathematical reason (feel free to ignore this point from a practical standpoint) but saying that you’re going to pick points uniformly from an interval containing an infinity of points does not make sense, since point probability is always zero. You can only assign probabilities to *intervals*. But in a MC simulation you can get around this by effectively claiming that each point representable in the floating point representation is actually pointing to an interval from that value to the next discrete value you can represent, to which the random number generator can assign non-zero probabilities

    • Eeshan on 11/12/2013 at 4:09 AM

      The first problem can perhaps be more easily understood in a 2-D world. Say I try to pick a chord of a circle with uniform probability assigned to each.

      We can pick the length of the chord by uniformly picking either:
      1. the perpendicular distance of the chord from the centre, or
      2. the angle the chord subtends at the centre

      It’s easy to see that the distribution density will be different in both cases (If not, try to calculate the probability of the chord being, say, shorter than the radius), while both methods intuitively feel like a ‘uniform’ distribution over all chords.

      • Arthur Charpentier on 11/12/2013 at 10:48 AM

        thanks Eeshan

        I just found a very interesting discussion on http://stats.stackexchange.com/questions/7977/how-to-generate-uniformly-distributed-points-on-the-surface-of-the-3-d-unit-sph about that !

        I will try to spend some time to get a deeper understanding, cause I find all this discussion interesting !!! I also have to find a old paper I found almost two years ago, by Richard Feynman : he did write a short paper on the uniform distribution on a sphere, in the 60s as far as I remember. I did not understand his problem (I have to stress here that my knowledge in theoretical physics is extremely limited), but now, I start to see why it might be more complicated than I thought…

  2. Dominique on 07/12/2013 at 8:56 AM

    I didn’t deeply investigate how to program à nice Mont Carlo strategy, but it seems that if you parametrize the points in sphere with (is it the case?)
    x=cos(u)*cos(v)
    y=sin(u)*cos(v)
    z=sin(v)
    and if you sample u and v uniformly, you will NOT get an uniform pattern of points. You will see some concentrations in the poles.
    R code :
    set.seed=1
    require(rgl)
    n<-10000
    U<-runif(n,0,2*pi)
    V<-runif(n,-pi/2,pi/2)
    x<-cos(U)*cos(V)
    y<-sin(U)*cos(V)
    z<-sin(V)
    plot3d(x,y,z,col=rgb(0.5,0,0),size=1)

    • Arthur Charpentier on 07/12/2013 at 9:52 AM

      that would be a good explanation, indeed ! this might explain why I overestime the probability ! Let me try to use the sphericallity of the multivariate Gaussian distribution instead

Leave a Reply

Categories

Students

Keep calm…

Keep calm…

Research blogs