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 points uniformly distributed on a sphere. What is the probability that the
points lie on a same hemisphere, for some hemisphere (there is no south or north here) ?
Analogously, what is the probability to see the 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
. But it was difficult. I mean, given
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.
- Draw
points on the unit sphere, which simply means generate
- Try to find
such that, after a rotation (with angle
) all the points lie in the upper part (the North hemisphere, for instance)
So the question here is simply : is
equal to ? 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 , 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.
- Draw
random variables
and
- Try to find
such that, after a rotation (with angles
)
As mentioned by Dominique, using this technique, points are not uniformly distributed on the sphere. Instead we can use
- Draw
- Get the polar coordinates, and try to find
So the question here is simply : is
(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 , the probability that the
points (uniformly distributed on a sphere) lie on a same hemisphere is exactly
p=function(d,n) .5^(n-1) * sum(choose(n-1,0:(d-1)))
Note that Leonard Savage proved (a few years before) that
(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 !
OpenEdition suggests that you cite this post as follows:
Arthur Charpentier (December 7, 2013). Random points on the Earth. Freakonometrics. Retrieved September 19, 2024 from https://doi.org/10.58079/ousw
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
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.
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…
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)
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