Wendel’s theorem on the sphere

Solution. By rotational symmetry, a uniform random point $P$ on the sphere and its antipode $-P$ are identically distributed. For four independent uniform points $P_1, P_2, P_3, P_4$, the $2^4 = 16$ configurations $$(\varepsilon_1 P_1, \varepsilon_2 P_2, \varepsilon_3 P_3, \varepsilon_4 P_4), \qquad \varepsilon_i \in \{\pm 1\},$$ obtained by independently reflecting each point through the centre, all have the same joint distribution.

Claim: among these $16$ configurations, exactly $2$ produce a tetrahedron containing the centre. Given the claim, the probability of interest is $2/16 = 1/8$.

Proof of the claim. The tetrahedron with vertices $Q_1, Q_2, Q_3, Q_4$ contains the origin $O$ iff no open hemisphere contains all four vertices. Fix $P_1, \ldots, P_4$ in general position, and consider the $16$ sign patterns $\varepsilon \in \{\pm 1\}^4$.

Form the $3\times4$ matrix $P$ whose columns are the $P_i$. General position means every three columns are linearly independent, so $P$ has rank $3$ and its nullspace is one-dimensional. Choose a nonzero vector $a$ with $Pa=0$. Each $a_i$ is nonzero: otherwise the other three columns would be dependent.

For the sign pattern $\varepsilon$, the origin lies in the convex hull iff there are weights $w_i\geq0$, with $\sum_i w_i=1$, such that $\sum_i w_i\varepsilon_iP_i=0$. Thus $(w_i\varepsilon_i)_i$ must be a scalar multiple of $a$. Since every $a_i\ne0$, all $w_i$ must be positive, and this is possible exactly when $\varepsilon_i=\operatorname{sign}(a_i)$ for every $i$, or when all those signs are reversed. Conversely, for either pattern the weights $w_i=|a_i|/\sum_j|a_j|$ work. Exactly two of the sixteen patterns therefore contain the origin.

Hence $$\boxed{\mathbb{P}(O \in \text{tetrahedron}) = \frac{2}{16} = \frac{1}{8}.} \qedhere$$

import numpy as np
# Uniform on sphere via normalised Gaussian
pts = np.random.randn(3, 4, 10**6)
pts /= np.linalg.norm(pts, axis=0)
def svol(a, b, c, d):
    return np.einsum('i...,i...->...', b-a, np.cross((c-a).T, (d-a).T).T)
O = np.zeros_like(pts[:, 0])
s0 = svol(pts[:, 0], pts[:, 1], pts[:, 2], pts[:, 3])
s = [svol(O, pts[:, 1], pts[:, 2], pts[:, 3]),
     svol(pts[:, 0], O, pts[:, 2], pts[:, 3]),
     svol(pts[:, 0], pts[:, 1], O, pts[:, 3]),
     svol(pts[:, 0], pts[:, 1], pts[:, 2], O)]
contains = np.all([np.sign(si) == np.sign(s0) for si in s], axis=0)
print(f"sim: {contains.mean():.5f}   exact: 0.125")
# sim: 0.12548   exact: 0.125