Uniform and Stratified Sampling from Manifolds in Euclidean Space

Set up

This chunk attaches the installed version of the package, saves the default graphics settings and initializes the random number generator:

library(tdaunif)
oldpar <- par(no.readonly = TRUE)
set.seed(0)

Volume-preserving parameterizations

In principle, any parameterization \(f : S \to X\) can be viewed as the composition of two steps:

• a warping map \(\omega : S \to \tilde{S}\) from the parameter space to a “warped” parameter space, and
• a volume-preserving map \(\phi : \tilde{S} \to X\) from the “warped” parameter space to the coordinate space.

In the best case, the warping map can be inverted and the volume-preserving map, itself a parameterization, extracted as \(\psi = \phi \circ \omega^{-1}\). Then points can be uniformly sampled from the warped parameter space and mapped, without distortion, to the manifold.

The upper hemisphere of \(\mathbb{S}^2\) is such a case. The function sample_2hemisphere() follows Arvo’s procedure (p. 50–51), which begins with a conventional parameterization to spherical coordinates and extracts this area-preserving parameterization from it:

\[\begin{align*} h : [0,1]\times[0,1] & \to \mathbb{R}^3 \\ \left[\begin{array}{c} u \\ v \end{array}\right] & \mapsto \left[\begin{array}{c} x \\ y \\ z \end{array}\right] = \left[\begin{array}{c} \sqrt{u (2 - u)}\cos{2\pi v} \\ \sqrt{u (2 - u)}\sin{2\pi v} \\ 1 - u \end{array}\right] \end{align*}\]

hemi <- sample_2hemisphere(180L)
pairs(hemi, asp = 1, pch = 19, cex = .5)

Notice the classical pattern of random clustering: Samples from uniform distributions are not themselves uniformly distributed (in the way that lattices are), so they don’t match our intuitive expectations. If we want greater uniformity, volume-preserving maps offer a solution.

sq1 <- sample_strat_square(100L, bins = 1L)
sq4 <- sample_strat_square(100L, bins = 4L)
sq10 <- sample_strat_square(100L, bins = 10L)
par(mfrow = c(1, 3))
plot(sq1, asp = 1, pch = 19, cex = .5)
plot(sq4, asp = 1, pch = 19, cex = .5)
plot(sq10, asp = 1, pch = 19, cex = .5)

hemi10 <- sample_2hemisphere(225L, bins = 10L)
pairs(hemi10, asp = 1, pch = 19, cex = .5)

Rejection sampling

In lieu of an analytic solution to the warping problem, rejection sampling provides a fallback that shifts the burden to additional sampling. The idea, explained in detail by Diaconis, Holmes, and Shahshahani and demonstrated in more detail in this blog post, is to oversample points throughout the parameter space and discard them with a probability that compensates for local distortions. Since the rejection is probabilistic, we can’t know in advance how large a sample will be necessary, so we repeatedly sample and discard until the total number of points meets or exceeds the desired number (then discard the excess points).

One application in this package is to sample from a famous class of embeddings of the circle \(\mathbb{S}^1\) in \(\mathbb{R}^3\): the trefoil, one of many such classes called knots that cannot be continuously deformed in \(\mathbb{R}^3\) to a geometric circle (the “unknot”).

\[\begin{align*} g : \mathbb{S}^1 & \to \mathbb{R}^3 \\ \left[\begin{array}{c} \theta \end{array}\right] & \mapsto \left[\begin{array}{c} x \\ y \\ z \end{array}\right] = \left[\begin{array}{c} \sin\theta + 2\sin 2\theta \\ \cos\theta - 2\cos 2\theta \\ -(\sin 3\theta) \end{array}\right] \end{align*}\]

We generate a sample below using the function sample_trefoil() and plot it in the three coordinate planes of \(\mathbb{R}^3\). The parameterization above suggests a sort of symmetry within the \(xy\)-plane, and the plot shows that the embedding is “knotted” by way of undulations along the \(z\)-axis:

tref <- sample_trefoil(120L)
pairs(tref, asp = 1, pch = 19, cex = .5)

disk_parameterization <- function(s, t) cbind(
  x = s,
  y = t * sqrt(1 - s^2)
)
disk_jacobian <- function(s, t) sqrt(1 - s^2)
disk_sampler <- make_rejection_sampler(
  parameterization = disk_parameterization,
  jacobian = disk_jacobian,
  min_params = c(-1, -1), max_params = c(1, 1),
  max_jacobian = 1
)
disk1 <- disk_sampler(360L)

disk2 <- sample_disk(360L)
par(mfrow = c(1L, 2L))
plot(disk1, asp = 1, pch = 19, cex = 1/3)
plot(disk2, asp = 1, pch = 19, cex = 1/3)

par(oldpar)