How to choose your next holidays destination – Uniform distribution on a sphere

November 14, 2012
By

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

If you want to choose randomly your next holidays destination, you are likely to process in a way which is certainly biased. Especially if you choose randomly the latitude and the longitude. A bit like they do in this lovely advertising. (For those of you who do not speak French, this is about a couple who won the national gamble prize and have to decide their next travel. The husband randomly pick Australia and the wife is complaining : "Not again!" she says.) So let me help you to chose your next travel!


If we were able to generate uniformly distributed variables on [0 ; 1], we could easily generate variables on other space such as [0 ; 1] x [0 ; 1] which is a square of length 1. It can be simply done by using two uniform on [0 ; 1] and independent variables  X1 and X2 and to consider X = (X1 , X2). Then X is a random variable evenly distributed on the square [0 ; 1] x [0 ; 1].


 However, this generation may be a bit more complicated if we work with more complicated shapes. Such as, for example, a sphere. Indeed, for a sphere of a given radius, say R = 1, it is possible to project the variables of a square on the sphere.

The first method I thought of, is wrong. I was tempted to generate two variables Y1 and Y2 where Y1 follows a uniform on  [0 ; Pi] and Y2 follows a uniform on [-Pi ; Pi]. In other words :
- Y1 = Pi * X1
- Y2 = 2Pi * (X2 - 0.5).
 Then I wanted to consider (Y1, Y2) as the spherical coordinated of the Sphere.

However, this method does NOT generate a uniform on the sphere. Indeed, it has a tendency to over generate north and south poles. The reason is simple, this method generate, in average the same amount of variable in each latitude. North and South poles are of smaller area than Equator. Therefore, the closer we are from the Equator latitude the less variables there are.

On the following graph, we observe a high density of point in the "middle" of the sphere, this is the North Pole. It shows that this method does not offer a uniform distribution.

A wrong method to generate evenly distributed variables on the sphere. The North Pole  (in the middle of the graph) as well as the South Pole are over represented.


Many methods have been proposed to generate evenly distributed random on a sphere. We propose one of them here. We consider the couple z = (u,v) defined as :

- u = 2 * Pi * X1
- v = arc-cos(2 * (X2 - 0.5))

A correct simulation of a uniform distribution on a sphere. There is no over represented area.
In this case, a theorem shows that z = (u , v) generate evenly distributed variables. It can be observed on the next graph. There is no irregularity in the distribution of the random variables.

Other methods exist. I like this one since it is really simple and use a uniform distribution at the beginning. The idea of this post is to show that generated a uniform distribution can be adapted to many shapes and cases. However, to do so, a previous analytical study has to be done to find the correct transformation. 

So this method would help you to avoid being too many times in the chilly places such as North Pole and South Pole since it does not over represent the extreme latitudes.

The program (R) :

# import package to plot in 3D
install.packages("rgl", dependencies = TRUE)
library(rgl)

################################################################
# Uniform distribution in a square
################################################################

size = 10000
x1 = runif(size)
x2 = runif(size)

# the option pch = '.' change the symbol for the graph into dot.
# cex = 2 doubles the size of the dots
plot(x1,x2, col = 'blue', pch = '.', cex = 2)

################################################################
# Wrong solution for the sphere
################################################################

y1 = pi * x1
y2 = 2* pi * (x2-0.5)
y = matrix(0, nrow = 2, ncol = size)
y[1,] = y1
y[2,] = y2
plot(y1, y2)

# This function transform the spherical coordinates into cartesian coordinates
sphereToCartesian = function(matrice){
  x= matrix(0,nrow = 3, ncol = length(matrice[1,]))
  x[1,] = sin(matrice[2,]) * cos(matrice[1,])
  x[2,] = sin(matrice[2,]) * sin(matrice[1,])
  x[3,] = cos(matrice[2,])
  return(x)
}

a = sphereToCartesian(y)
plot3d(x = a[1,], y = a[2,], z = a[3,]))

################################################################
# Correct solution for the sphere
################################################################

uniformSphere = function(length){
  x1 = runif(length, 0,1)
  x2 = runif(length, 0,1)
  u = 2*pi*x1
  v = acos(2*x2- 1)
  z =matrix(0, ncol = length, nrow = 2)
  z[1,] = u
  z[2,] = v
  return(z)
}

z = uniformSphere(size)
b = sphereToCartesian(z)
plot3d(x = b[1,], y = b[2,], z = b[3,])


To leave a comment for the author, please follow the link and comment on his blog: ProbaPerception.

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,ecdf, 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...

Tags: , , , , ,

Comments are closed.