# Random points on the Earth

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 $n$ points uniformly distributed on a sphere. What is the probability that the $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 $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 $n$. But it was difficult. I mean, given $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 $n$ points on the unit sphere, which simply means generate $n$ random variables $U_1,\ldots ,U_n\sim\mathcal{U}([0,2\pi])$
2. Try to find $\theta\in[0,2\pi]$ such that, after a rotation (with angle $\theta$) all the points lie in the upper part (the North hemisphere, for instance)

So the question here is simply : is

$\max_{\theta\in[0,2\pi]}\left\{\sum_{i=1}^n \boldsymbol{1}(\sin(U_i+\theta)\geq 0)\right\}$

equal to $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 $\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 $n$ points on the unit sphere, which simply means generate $2n$ random variables $U_1,\ldots ,U_n\sim\mathcal{U}([0,2\pi])$ and $V_1,\ldots,V_n\sim\mathcal{U}([-\pi/2,\pi/2])$
2. Try to find $\theta$ and $\varphi$ such that, after a rotation (with angles $\theta$ and $\varphi$) all the points lie in some given part (say  $\{x\geq0\}$)

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

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

So the question here is simply : is

$\max_{(\theta,\varphi)\in[0,2\pi]\times[-\pi/2,\pi/2]}\left\{\sum_{i=1}^n \boldsymbol{1}(\sin(U_i+\theta)\cos(V_i+\varphi)\geq 0) \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 $d$, the probability that the $n$ points (uniformly distributed on a sphere) lie on a same hemisphere is exactly

$p(d,n)= \frac{1}{2^{n-1}}\sum_{i=0}^{d-1} \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

$p(d,d+1)=1 -\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 !

# Circular or spherical data, and density estimation

I few years ago, while I was working on kernel based density estimation on compact support distribution (like copulas) I went through a series of papers on circular distributions. By that time, I thought it was something for mathematicians working on weird spaces…. but during the past weeks, I saw several potential applications of those estimators.

• circular data density estimation

Consider the density of an angle say, i.e. a function such that

with a circular relationship, i.e. . It can be seen as an invariance by rotation.
von Mises proposed a parametric model in 1918 (see here or there), assuming that

where is Bessel modified function of order 1,

(which is simply a normalization parameter). There are two parameters here, (some concentration parameter) and mu a direction.
From a series of observed angles, the maximum likelihood estimator for kappa is solution of

where

and

and where , where those functions are modified Bessel functions. Well, that estimator is biased, but it is possible to improve it (see here or there). This can be done easily in R (actually Jeff Gill – here – used that package in several applications). But I am not a big fan of that technique….

• density estimation for hours on simulated data

A nice application can be on the estimation of the daily density of a temporal events (e.g. phone calls as we’ll see later on, or email arrival time). Let is the time (in hours) for the th observation (the th phone call received). Then set

The time is now seen as an angle. It is possible to consider the equivalent of an histogram,

set.seed(1)
library(circular)
X=rbeta(100,shape1=2,shape2=4)*24
Omega=2*pi*X/24
Omegat=2*pi*trunc(X)/24
plot(Ht, stack=FALSE, shrink=1.3, cex=1.03,
points(Ht, rotation = "clock", zero =c(rad(90)),
col = "1", cex=1.03, stack=TRUE )

rose.diag(Ht-pi/2,bins=24,shrink=0.33,xlim=c(-2,2),ylim=c(-2,2),
axes=FALSE,prop=1.5)

or a kernel based estimation of the density (the gray line on the right).

circ.dens = density(Ht+3*pi/2,bw=20)
plot(Ht, stack=TRUE, shrink=.35, cex=0, sep=0.0,
axes=FALSE,tol=.8,zero=c(0),bins=24,
xlim=c(-2,2),ylim=c(-2,2), ticks=TRUE, tcl=.075)
lines(circ.dens, col="darkgrey", lwd=3)
text(0,0.8,"24", cex=2); text(0,-0.8,"12",cex=2);
text(0.8,0,"6",cex=2); text(-0.8,0,"18",cex=2)

The code looks rather simple. But I am not very comfortable using codes that I do not completely understand. So I did my own. The first step was to get a graph similar to the one we have on the right, except that I prefer my own kernel based estimator. The idea is that instead of estimating the density on , we estimate it on the sample . Then we multiply by 3 to get the density only on . For the bandwidth, I took the same as the one that we would have taken on

The code is simply the following

U=seq(0,1,by=1/250)
O=U*2*pi
U12=seq(0,1,by=1/24)
O12=U12*2*pi
X=rbeta(100,shape1=2,shape2=4)*24
OM=2*pi*X/24
XL=c(X-24,X,X+24)
d=density(X)
d=density(XL,bw=d$bw,n=1500) I=which((d$x>=6)&(d$x<=30)) Od=d$x[I]/24*2*pi-pi/2
Dd=d$y[I]/max(d$y)+1

plot(cos(O),-sin(O),xlim=c(-2,2),ylim=c(-2,2), type="l",axes=FALSE,xlab="",ylab="") for(i in pi/12*(0:12)){ abline(a=0,b=tan(i),lty=1,col="light yellow")} segments(.9*cos(O12),.9*sin(O12),1.1*cos(O12),1.1*sin(O12)) lines(Dd*cos(Od),-Dd*sin(Od),col="red",lwd=1.5) text(.7,0,"6"); text(-.7,0,"18") text(0,-.7,"12"); text(0,.7,"24") R=1/24/max(d$y)/3+1 lines(R*cos(O),R*sin(O),lty=2) Note that it is possible to stress more (visually) on hours having few phone calls, or a lot (compared with an homogeneous Poisson process), e.g. plot(cos(O),-sin(O),xlim=c(-2,2),ylim=c(-2,2), type="l",axes=FALSE,xlab="",ylab="") for(i in pi/12*(0:12)){ abline(a=0,b=tan(i),lty=1,col="light yellow")} segments(2*cos(O12),2*sin(O12),1.1*cos(O12),1.1*sin(O12), col="light grey") segments(.9*cos(O12),.9*sin(O12),1.1*cos(O12),1.1*sin(O12)) text(.7,0,"6") text(-.7,0,"18") text(0,-.7,"12") text(0,.7,"24") R=1/24/max(d$y)/3+1
lines(R*cos(O),R*sin(O),lty=2)
AX=R*cos(Od);AY=-R*sin(Od)
BX=Dd*cos(Od);BY=-Dd*sin(Od)
COUL=rep("blue",length(AX))
COUL[R<Dd]="red"
CM=cm.colors(200)
a=trunc(100*Dd/R)
COUL=CM[a]
segments(AX,AY,BX,BY,col=COUL,lwd=2)
lines(Dd*cos(Od),-Dd*sin(Od),lwd=2)

We get here those two graphs,

To be honest, I do not really like that representation – even if it looks nice. If we compare that circular representation to a more classical one (from 0:00 till 23:59 one the graph on the left, below), I do have a problem to interpret the areas in blue and pink.

density of wind direction

On the left, we compare two densities, so the area in pink is the same as the area in blue. But here, it is no longer the case: the area in pink is always larger to the one in blue. So it might help so see when we have a difference, but there is a scaling issue that we cannot discuss further… But less us see if we can use that estimation technique to several problems.

A standard application when studying angles is wind direction. For instance, in Montréal, it is possible to find hourly observations, starting in 1974 (we just need a R robot to pick up the information, but I’ll tell more about that in another post, someday). Here, we have directly an angle. So we can use a code rather similar to the one used above to estimate the distribution of wind direction in Montréal.

density of 911 phone calls

Note that our estimate is consistent with several graphs that can be found on meteorological websites (e.g. the one above on the right, that was found here).

In a recent post (here) I wanted to check about the “midnight crime” myth, using hours of 911 phone calls in Montréal.

That was for all phone calls. But if we look more specifically, for burglaries, we have the distribution on the left, and for conflicts the one on the right

We do clearly observe that gun shots occur a bit before midnight. See also here for another study, but this time in NYC (thanks @PAC for the link).while for gun shots, we have the distribution on the left, and for “troubles” (basically people making too much noisy in parties) or “noise” the one on the right

• density of earth temperatures, or earthquakes

Of course it is also possible to work in higher dimension. Before, we went from densities on to densities on the unit circle . But similarly, it is possible to go from to the unit sphere . A nice application being global climate studies,

The idea being that point on the left above are extremely close to the one on the right. An application can be e.g. on earthquakes occurrence. Data can be found here.

library(ks)
X=cbind(EQ$Longitude,EQ$Latitude)
Hpi1 = Hpi(x = X)
DX=kde(x = X, H = Hpi1)
library(maps)
map("world")
cbind(X[,1]+360,X[,2]+180),cbind(X[,1]-360,X[,2]+180), cbind(X[,1],X[,2]-180),cbind(X[,1]+360, X[,2]-180),cbind(X[,1]-360,X[,2]-180)) DY=kde(x = Y, H = Hpi1) library(maps) plot (DY,add=TRUE,col="purple")