Uniform distribution on sphere

So, we want to generate uniformly distributed random numbers on a unit sphere. This came up today in writing a code for molecular simulations. Spherical coordinates give us a nice way to ensure that a point $(x, y, z)$

x = r sin (ϕ) cos (θ)

y = r sin (ϕ) sin (θ)

z = r cos (ϕ) .

In spherical coordinates, $r$

A tempting way to generate uniformly distributed numbers in a sphere is to generate a uniform distribution of $θ$

#include<random>
#include<cmath>
#include<chrono>

int main(int argc, char *argv[]) {
    // Set up random number generators
    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
    std::mt19937 generator (seed);
    std::uniform_real_distribution<double> uniform01(0.0, 1.0);

    // generate N random numbers
    int N = 1000;

    // the incorrect way
    FILE * incorrect;
    incorrect = fopen("incorrect.csv", "w");
    fprintf(incorrect, "Theta,Phi,x,y,z\n");
    for (int i = 0; i < N; i++) {
        // incorrect way
        double theta = 2 * M_PI * uniform01(generator);
        double phi = M_PI * uniform01(generator);
        double x = sin(phi) * cos(theta);
        double y = sin(phi) * sin(theta);
        double z = cos(phi);
        fprintf(incorrect, "%f,%f,%f,%f,%f\n", theta, phi, x, y, z);
    }
    fclose(incorrect);
}

The distribution in the $θ$

After mapping these points in the $θ$

The reason for this is that the area $d A$

So, close to the poles of the sphere ( $ϕ = 0$

Our goal is to find and then draw samples from the probability distribution $f (θ, ϕ)$

Let $v$

f (v) d A = 1 4 π d A = f ( θ , ϕ ) d θ d ϕ ,

since $f (v) d A$

Marginalizing the joint distribution to get the p.d.f of $θ$

f (θ) = \int π 0 f (θ, ϕ) d ϕ = 1 2 π

f (ϕ) = \int 2 π 0 f (θ, ϕ) d θ = sin ( ϕ ) 2 .

We see that $θ$

Now, how can we sample numbers $ϕ$

F (ϕ) = \int ϕ 0 f (ϕ^) d ϕ^= 1 2 ( 1 - cos ( ϕ ) ) .

Keep in mind that $F (ϕ)$

Let $U$

P r (U \leq F (ϕ)) = F (ϕ) .

This is a property of the uniform random variable $U [0, 1]$

P r (F - 1 (U) \leq ϕ) = F (ϕ) .

Aha! This shows that $F (ϕ)$

Generate a uniform random number $u$
Compute $ϕ$
Take this $ϕ$

In our case, $F^{- 1} (u) = \arccos (1 - 2 u)$

The algorithm below in C++ shows how to generate uniformly distributed numbers on the sphere using this method:

#include<random>
#include<cmath>
#include<chrono>

int main(int argc, char *argv[]) {
    // Set up random number generators
    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
    std::mt19937 generator (seed);
    std::uniform_real_distribution<double> uniform01(0.0, 1.0);

    // generate N random numbers
    int N = 1000;

    // the correct way
    FILE * correct;
    correct = fopen("correct.csv", "w");
    fprintf(correct, "Theta,Phi,x,y,z\n");
    for (int i = 0; i < N; i++) {
        // incorrect way
        double theta = 2 * M_PI * uniform01(generator);
        double phi = acos(1 - 2 * uniform01(generator));
        double x = sin(phi) * cos(theta);
        double y = sin(phi) * sin(theta);
        double z = cos(phi);
        fprintf(correct, "%f,%f,%f,%f,%f\n", theta, phi, x, y, z);
    }
    fclose(correct);
}

We then get the following distribution of points in the $(θ, ϕ)$

And, finally, a uniform distribution of points on the sphere.

Alternative method 1

An alternative method to generate uniformly disributed points on a unit sphere is to generate three standard normally distributed numbers $X$

The following Julia code implements this. We have to be careful in the case that the vector has a norm close to zero, in which we must worry about floating point precision by dividing by a very small number. This is the reason for the while loop.

n = 100

f_normal = open("normal.csv", "w")
write(f_normal, "x,y,z\n")

for i = 1:n
    v = [0, 0, 0]  # initialize so we go into the while loop

    while norm(v) < .0001
        x = randn()  # random standard normal
        y = randn()
        z = randn()
        v = [x, y, z]
    end
    
    v = v / norm(v)  # normalize to unit norm

    @printf(f_normal, "%f,%f,%f\n", v[1], v[2], v[3])
end

To prove this, note that the standard normal distribution is:

f (x) = 1 2 π--\sqrt e - 1 2 x 2 .

As $X$

f (x, y, z) = f (v) = (1 2 π--\sqrt e - 1 2 x 2 ) ( 1 2 π--\sqrt

With some algebra:

f (x, y, z) = f (v) = 1 ( 2 π ) 3 2 e - 1 2 ( x 2 + y 2 + z

This shows that the probability distribution of $v$

Alternative method 2

Credit to FX Coudert for pointing this out in a comment, another method is to generate uniformly distributed numbers in the cube $[- 1, 1]^{3}$

Alternative method 1 can be used to efficiently generate uniformly distributed numbers on a hypersphere– i.e. in higher dimensions. On the other hand, as the number of dimensions grows, the ratio of the volume of the edges of a cube to the volume of the ball inside of it grows (see Wikipedia article). Hence, a larger and larger fraction of the uniformly generated numbers in the cube will be rejected using Alternative method 2, and so this algorithm will be inefficient.

Mathemathinking

Math, Data, Machine learning, Chemical engineering

Generating uniformly distributed numbers on a sphere

Alternative method 1

Alternative method 2

导航

公告