Uniformly Sampling a Point Inside a 3D Unit Sphere

three ways with their corresponding codes

Introduction

When sampling a random point uniformly inside a 3D unit sphere, several methods can be used:

1. 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, \ \text{else repeat generating} \]
2. Inverse CDF Sampling (Polar)
   • Sample radius \( r \) from: \[ r = U_1^{\frac{1}{3}}, \ \text{where } \ U_1 \sim \operatorname{Unif}[0, 1] \]
   • Sample angles: \[ \theta \sim \operatorname{Unif}[0, 2\pi] \] \[ \phi = \cos^{-1}(U_2), \ \text{where } \ U_2 \sim \operatorname{Unif}[-1, 1] \] Note that the formula of \(\phi\) is not that 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 \]
3. 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

  \[\begin{aligned} 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} \\ \end{aligned}\]

  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

  \[\begin{aligned} 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} \\ \end{aligned}\]

  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

  \[\begin{aligned} 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}} \end{aligned}\]

  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

1. 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)

2. 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)

3. 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