Uniformly Sampling a Point Inside a 3D Unit Sphere
Published:
three ways with their corresponding codes
Introduction
When sampling a random point uniformly inside a 3D unit sphere, several methods can be used:
- Rejection Sampling:
- Generate a point in the cube: \[ (x, y, z) \sim \operatorname{Unif}[-1, 1]^3 \]
- Test if the point is inside the sphere: \[ \text{Accept if } x^2 + y^2 + z^2 \leq 1 \]
- Inverse CDF Sampling (Polar)
- Sample radius \( r \) from: \[ r = U_1^{\frac{1}{3}}, \quad U_1 \sim \operatorname{Unif}[0, 1] \]
- Sample angles: \[ \theta \sim \operatorname{Unif}[0, 2\pi] \] \[ \phi = \cos^{-1}(U_2), \quad U_2 \sim \operatorname{Unif}[-1, 1] \] Note that the formula of \(\phi\) is a little bit non-intuitive because the area of spherical cap using \(r\) and \(\phi\) is \[ A = 2 \pi r^2 (1-\cos\phi) \]
- Convert to Cartesian coordinates: \[ x = r \cos\theta \sin\phi, \quad y = r \sin\theta \sin\phi, \quad z = r \cos\phi \]
- Inverse CDF Sampling (Cartesian):
- This method can be challenging for beginners.
- Computing the conditional CDF is complex.
Inverse CDF Sampling in a Sphere
-
Step 1: We denote the unit ball by
\[B = \{(x, y, z) \in \mathbb{R}^3 : x^2 + y^2 + z^2 \le 1\}\]From geometry, we know
\[\text{Volume}(B) = \frac{4}{3}\pi\]Now consider the CDF \(F_X : [-1, 1] \to [0, 1]\) defined by
\[F_X(x) = \frac{\int_{-1}^x \pi(1 - t^2) \ dt}{\text{Volume}(B)} = \frac{\pi (x+1) - \frac{\pi}{3} (x^3 + 1)}{\frac{4}{3}\pi} = - \frac{1}{4} x^3 + \frac{3}{4} x + \frac{1}{2}\]Since \(F_X\) is continuous and strictly increasing, it has an inverse
\[F_X^{-1} : [0, 1] \to [-1, 1]\]Then for \(U_1 \sim \operatorname{Unif}[0, 1]\), we have \(F_X^{-1}(U_1) \sim F_X\), and we can get a sample of \(X = F_X^{-1}(U_1)\).
-
Step 2: When \(X\) is fixed, we can sample \(Y\) and \(Z\) uniformly from the disc in the plane orthogonal to the unit vector \((1, 0, 0)\).
We denote the disc by
\[D(x) = \{(y, z) \in \mathbb{R}^2 : y^2 + z^2 \le 1 - x^2\}\]From geometry, we know
\[\quad \text{Area}(D(x)) = \pi(1 - x^2)\]Now consider the CDF \(F_{Y \mid X} : [-\sqrt{1 - x^2}, \sqrt{1 - x^2}] \to [0, 1]\) defined by
\[F_{Y \mid X}(y \mid x) = \frac{\int_{-\sqrt{1 - x^2}}^y 2\sqrt{1 - x^2 - t^2} \ dt}{\text{Area}(D(x))} = \frac{y \sqrt{1 - x^2 - y^2}}{\pi(1 - x^2)} + \frac{\arcsin\left(\frac{y}{\sqrt{1 - x^2}}\right)}{\pi} + \frac{1}{2}\]Since \(F_{Y \mid X}\) is continuous and strictly increasing, it has an inverse
\[F_{Y \mid X}^{-1} : [0, 1] \to [-\sqrt{1 - x^2}, \sqrt{1 - x^2}]\]Then for \(U_2 \sim \operatorname{Unif}[0, 1]\), we have \(F_{Y \mid X}^{-1}(U_2) \sim F_{Y \mid X}\), and we can get a sample of \(Y = F_{Y \mid X}^{-1}(U_2)\).
-
Step 3: When \(X\) and \(Y\) are fixed, we can sample \(Z\) uniformly from the line segment in the plane orthogonal to the X-Y plane.
We denote the line segment by
\[L(x, y) = \{z \in \mathbb{R} : -\sqrt{1 - x^2 - y^2} \le z \le \sqrt{1 - x^2 - y^2}\}\]From geometry, we know
\[\text{Length}(L(x, y)) = 2\sqrt{1 - x^2 - y^2}\]Now consider the CDF \(F_{Z \mid X,Y} : [-\sqrt{1 - x^2 - y^2}, \sqrt{1 - x^2 - y^2}] \to [0, 1]\) defined by
\[F_{Z \mid X,Y}(z \mid x,y) = \frac{\int_{-\sqrt{1 - x^2 - y^2}}^z 1 \ dt}{\text{Length}(L(x, y))} = \frac{z + \sqrt{1 - x^2 - y^2}}{2\sqrt{1 - x^2 - y^2}}\]Since \(F_{Z \mid X,Y}\) is continuous and strictly increasing, it has an inverse
\[F_{Z \mid X,Y}^{-1} : [0, 1] \to [-\sqrt{1 - x^2 - y^2}, \sqrt{1 - x^2 - y^2}]\]Then for \(U_3 \sim \operatorname{Unif}[0, 1]\), we have \(F_{Z \mid X,Y}^{-1}(U_3) \sim F_{Z \mid X,Y}\), and we can get a sample of \(Z = F_{Z \mid X,Y}^{-1}(U_3)\).
Python Implementation
-
Rejection Sampling:
import numpy as np RNG = np.random.default_rng() unif = RNG.random def ball_reject(n: int) -> np.ndarray: i = 0 samples = [] while i < n: sample = -1 + 2 * unif(3) if np.linalg.norm(sample) < 1: samples.append(sample) i += 1 return np.array(samples) -
Inverse CDF Sampling (Polar):
import numpy as np RNG = np.random.default_rng() unif = RNG.random def ball_cdf_polar(n: int) -> np.ndarray: r = unif(n) ** (1/3) theta = 2 * np.pi * unif(n) phi = np.arccos(-1 + 2 * unif(n)) x = r * np.cos(theta) * np.sin(phi) y = r * np.sin(theta) * np.sin(phi) z = r * np.cos(phi) return np.moveaxis([x, y, z], 0, -1) -
Inverse CDF Sampling (Cartesian):
import numpy as np from typing import Callable RNG = np.random.default_rng() unif = RNG.random def sample_cdf_bounded(cdf: Callable[[float], float], a: float, b: float, epsilon: float = 0.01, size: int = 1) -> np.ndarray: result = np.zeros(size) count = 0 while count < size: target = unif() left = a right = b # make sure the error is within epsilon while right - left > epsilon: mid = (left + right) / 2 if cdf(mid) < target: left = mid elif cdf(mid) > target: right = mid else: break result[count] = mid count += 1 return result def F_X(x: float) -> float: return -0.25 * x**3 + 0.75 * x + 0.5 def F_Y_given_X(y: float, x: float) -> float: a = np.sqrt(1 - x**2) if a == 0: return 1.0 if y >= 0 else 0.0 term1 = (y * np.sqrt(a**2 - y**2)) / (np.pi * a**2) term2 = np.arcsin(y / a) / np.pi return term1 + term2 + 0.5 def F_Z_given_XY(z: float, x: float, y: float) -> float: b = np.sqrt(1 - x**2 - y**2) if b == 0: return 1.0 if z >= 0 else 0.0 return (z + b) / (2 * b) def ball_cdf_cartesian(n: int) -> np.ndarray: res = np.zeros((n, 3)) for i in range(n): # Step 1: Sample X from F_X using inverse CDF x = sample_cdf_bounded(F_X, -1.0, 1.0, epsilon=1e-6, size=1)[0] # Step 2: Sample Y from F_{Y|X} using inverse CDF a = np.sqrt(1 - x**2) if a == 0: y = 0.0 else: def F_Y_given_X_fixed(y_val): return F_Y_given_X(y_val, x) y = sample_cdf_bounded(F_Y_given_X_fixed, -a, a, epsilon=1e-6, size=1)[0] # Step 3: Sample Z from F_{Z|XY} using inverse CDF b = np.sqrt(1 - x**2 - y**2) if b == 0: z = 0.0 else: def F_Z_given_XY_fixed(z_val): return F_Z_given_XY(z_val, x, y) z = sample_cdf_bounded(F_Z_given_XY_fixed, -b, b, epsilon=1e-6, size=1)[0] res[i] = [x, y, z] return res