In the previous post I outlined how to simulate or sample a homogeneous Poisson point process on the surface of a sphere. Now I will consider a homogeneous Poisson point process on the \((n-1)-\) sphere, which is the surface of the Euclidean ball in \(n\) dimensions.
This is a short post because it immediately builds off the previous post. For positioning the points uniformly, I will use Method 2 from that post, which uses normal random variables, as it immediately gives a fast method in \(n\) dimensions.
I wrote this post and the code more for curiosity than any immediate application. But simulating a Poisson point process in this setting requires placing points uniformly on a sphere. And there are applications in that, such as Monte Carlo integration methods, as mentioned in this post, which nicely details different sampling methods.
Steps
As is the case for other shapes, simulating a Poisson point process requires two steps.
Number of points
The number of points of a Poisson point process on the surface of a sphere of radius \(r>0\) is a Poisson random variable. The mean of this random variable is \(\lambda S_{n-1}\), where \(S_{n-1}\) is the surface area of the sphere. For a ball embedded in \(n\) dimension, the area of the corresponding sphere is given by
$$S_{n-1} = \frac{2 \pi ^{n/2} }{\Gamma(n/2)} r^{n-1},$$
where \(\Gamma\) is the gamma function, which is a natural generalization of the factorial. In MATLAB, we can simply use the function gamma. In Python, we need to use the SciPy function scipy.special. gamma.
Locations of points
For each point on the sphere, we generate \(n\) standard normal or Gaussian random variables, say, \(W_1, \dots, W_n\), which are independent of each other. These random variables are the Cartesian components of the random point. We rescale the components by the Euclidean norm, then multiply by the radius \(r\).
For \(i=1,\dots, n\), we obtain
$$X_i=\frac{rW_i}{(W_1^2+\cdots+W_n^2)^{1/2}}.$$
These are the Cartesian coordinates of a point uniformly scattered on a sphere with radius \(r\) and a centre at the origin.
How does it work?
In the post on the circle setting, I gave a more detailed outline of the proof, where I said the method is like the Box-Muller transform in reverse. The joint density of the normal random variables is from a multivariate normal distribution with zero correlation. This joint density a function of the Cartesian equation for a sphere. This means the density is constant on the sphere, implying that the angle of the point \((W_1,\dots, W_n)\) will be uniformly distributed.
The vector formed from the normal variables \((W_1,\dots,W_n)\) is a random variable with a chi distribution. But the final vector, which stretches from the origin to the point \((X_1,\dots,X_n)\), has length one, because we rescaled it with the Euclidean norm.
Code
The code for all my posts is located online here. For this post, the code in MATLAB and Python is here.
Further reading
I recommend this blog post, which discusses different methods for randomly placing points on spheres and inside spheres (or, rather, balls) in a uniform manner. (Embedded in two dimensions, a sphere is a circle and a ball is a disk.)
Our Method 2 for positioning points uniformly, which uses normal variables, comes from the paper:
- 1959, Muller, A note on a method for generating points uniformly on n-dimensional spheres.
Two recent works on this approach are the papers:
- 2010, Harman and Lacko, On decompositional algorithms for uniform sampling from -spheres and -balls;
- 2017, Voelker, Gosman, Stewart, Efficiently sampling vectors and coordinates.