Return to JUNTO

JUNTO Practice: Probability, Points on the Unit Sphere

Discussed on August 12, 2020.

From "Twenty problems in probability":

Four points are chosen on the unit sphere. What is the probability that the origin lies inside the tetrahedron determined by the four points?

Please provide an answer and a justification.

For more information, see 3Blue1Brown's video "The hardest problem on the hardest test".

Solutions

Click to see:

Oscar Martinez

I believe that the probability that the origin lies inside the tetrahedron determined by the four points is P ≈ 1/8.

import numpy as np
from numba import jit


@jit
def get_vectors(n_vectors=4, n_dim=3):
    vectors = np.random.randn(n_vectors, n_dim)
    vectors_norm = vectors.transpose() / np.linalg.norm(
        vectors.transpose(), axis=0
    )
    vectors_norm = vectors_norm.transpose()
    return vectors_norm


@jit
def new_coordinate_system(vectors, p=np.array([0, 0, 0])):
    mat = vectors[1:] - vectors[0]
    mat = mat.T
    inv_mat = np.linalg.inv(mat)

    p_new_coord = inv_mat.dot(p - vectors[0])

    return (
        np.all(p_new_coord >= 0)
        and np.all(p_new_coord <= 1)
        and np.sum(p_new_coord) <= 1
    )


@jit
def simulate(n=10000):
    sims = [
        new_coordinate_system(get_vectors())
        for _ in range(n)
    ]
    return np.mean(sims)


simulate(n=1_000_000)
0.124961

John Lekberg

I believe that the probability that the origin lies inside the tetrahedron determined by the four points is P ≈ 0.125.


I wrote a Python function, random_origin_in_tetrahedron, that simulates creating a tetrahedron from points on the sphere and then checks if the origin lies inside:

import numpy as np
from numpy import sin, cos, pi


def random_origin_in_tetrahedron(N):
    """Sample this boolean random variable:
    
    > The origin lies inside the tetrahedron
    > determined by 4 randomly selected points
    > on the unit sphere.
    
    N -- int. the number of samples to take.
        E.g. 537.
    
    Here's how the code works:
    
    1. Generate spherical coordinates (points on
       the sphere).
    2. Convert spherical coordinates to cartesian
       coordinates.
    3. Use the cartesian coordinates to define a
       tetrahedron.
    4. Convert the origin (0, 0, 0) into
       barycentric coordinates using the
       tetrahedron.
    5. The origin lies inside the tetrahedron iff
       the barycentric coordinates are all positive.
    """
    π = pi

    # Spherical coordinates.
    r = 1
    θ = np.random.uniform(0, π, [N, 4])
    φ = np.random.uniform(0, 2 * π, [N, 4])

    # Cartesian coordinates.
    x = r * sin(θ) * cos(φ)
    y = r * sin(θ) * sin(φ)
    z = r * cos(θ)

    # Absolute barycentric coordinates.
    R = np.stack([x, y, z, np.ones([N, 4])], axis=1)
    origin = np.array([[0, 0, 0, 1]])
    λ = np.linalg.solve(R, origin)

    # Does the tetrahedron contain the origin?
    contained = (λ > 0).all(axis=1)

    return contained


random_origin_in_tetrahedron(1_000_000).mean()
0.125061

Here's a version of the function that only uses ASCII identifiers:

def random_origin_in_tetrahedron(N):
    # Spherical coordinates.
    r = 1
    theta = np.random.uniform(0, pi, [N, 4])
    phi = np.random.uniform(0, 2 * pi, [N, 4])

    # Cartesian coordinates.
    x = r * sin(theta) * cos(phi)
    y = r * sin(theta) * sin(phi)
    z = r * cos(theta)

    # Absolute barycentric coordinates.
    R = np.stack([x, y, z, np.ones([N, 4])], axis=1)
    origin = np.array([[0, 0, 0, 1]])
    L = np.linalg.solve(R, origin)

    # Does the tetrahedron contain the origin?
    contained = (L > 0).all(axis=1)

    return contained

Daniel Bassett

My answer is 1/8, resulting from the simulation below. Using linear algebra, if the determinants of the random points and center point have the same sign, we know that the center is within the tetrahedron.

import numpy as np

expected = []
np.random.seed(seed=2)


def is_valid(x, y, z, center):
    matrix = np.column_stack((x, y, z, [1, 1, 1, 1]))
    determinant = np.linalg.det(matrix)
    for point in range(4):
        sim_matrix = matrix.copy()
        sim_matrix[point] = np.hstack((center, 1))
        if (
            np.sign(np.linalg.det(sim_matrix) * determinant)
            < 0
        ):
            result = False
            break
        else:
            result = True
    return result


smin = 1
smax = 4
sval = np.logspace(smin, smax, 100, dtype=int)
for sim in sval:
    valids = 0
    for points in range(sim):
        phi = np.random.uniform(
            low=0.0, high=2 * np.pi, size=4
        )
        theta = np.random.uniform(
            low=0.0, high=np.pi, size=4
        )
        p1 = np.sin(theta) * np.cos(phi)
        p2 = np.sin(theta) * np.sin(phi)
        p3 = np.cos(theta)
        if is_valid(p1, p2, p3, [0, 0, 0]):
            valids += 1

    expected.append(valids / sim)

print(sum(expected) / len(expected))
0.1264443725491392