# JUNTO Practice: Probability, Points on the Unit Sphere

Discussed on August 12, 2020.

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

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
``````