5
$\begingroup$

The present question follows on from some other questions on this site asking how to generate uniform points inside a disc (see e.g., here, here and here). The natural extension of that problem is to generate points inside an $m$-dimensional ball with centre $\mathbf{c} \in \mathbb{R}^m$ and radius $r \geqslant 0$. That is, we want to generate IID random variables from the following distribution:

$$\mathbf{X} \sim \text{U}(\mathcal{B}(\mathbf{c},r)) \quad \quad \quad \mathcal{B}(\mathbf{c},r) \equiv \Big\{ \mathbf{x} \in \mathbb{R}^m \Big| ||\mathbf{x} - \mathbf{c}|| \leqslant r \Big\}.$$

How do we generate IID uniform points on this space? Is there a simple way to program this?

$\endgroup$
  • 1
    $\begingroup$ Ray Koopman proposed a nice solution to me here stats.stackexchange.com/q/79919/3277. $\endgroup$ – ttnphns Aug 6 at 0:35
  • $\begingroup$ Although the duplicate question concerns only the case $m=2,$ Aksakal's answer there explicitly addresses the general case. $\endgroup$ – whuber Oct 4 at 12:30
6
$\begingroup$

A simple and efficient method for this problem uses a variation of the well-known Box-Mueller transform, which connects the normal distribution to the uniform distribution on a ball. If we generate a random vector $\mathbf{Z} = (Z_1,...,Z_m)$ composed of IID standard normal random variables and a random variable $U \sim \text{U}(0,1)$ (independent of the first random vector) then we can construct the uniform point of interest as:

$$\mathbf{X} = \mathbf{c} + r \cdot U^{1/m} \cdot \frac{\mathbf{Z}}{||\mathbf{Z}||}.$$

In the code below we create an R function called runifball which implements this method. The function allows the user to generate n random vectors that are points on a ball with arbitrary centre, radius and dimension.

runifball <- function(n, centre = 0, center = centre, radius = 1) {
  
  #Check inputs
  if (!missing(centre) && !missing(center)) {
  if (sum((centre - center)^2) < 1e-15) { 
                 warning("specify 'centre' or 'center' but not both") } else {
                    stop("Error: specify 'centre' or 'center' but not both") } }
  if (radius < 0) { stop("Error: radius must be non-negative") }
  
  #Create output matrix
  m   <- length(center)
  OUT <- matrix(0, nrow = m, ncol = n)
  rownames(OUT) <- sprintf("x[%s]", 1:m)
  
  #Generate uniform values on circle
  UU  <- runif(n, min = 0, max = radius)
  ZZ  <- matrix(rnorm(n*m), nrow = m, ncol = n)
  for (i in 1:n) {
    OUT[, i] <- center + radius*UU[i]^(1/m)*ZZ[, i]/sqrt(sum(ZZ[, i]^2)) }
  
  OUT }

Here is an example using this function to generate random points uniformly over a two-dimensional disk. The plot shows that the points are indeed uniform over the specified ball.

#Generate points uniformly on a disk
set.seed(1)
n      <- 10^5
CENTRE <- c(5, 3)
RADIUS <- 3
UNIF   <- runifball(n, centre = CENTRE, radius = RADIUS)

#Plot the points
plot(UNIF, 
     col = rgb(0, 0, 0, 0.05), pch = 16, asp = 1,
     main = 'Points distributed uniformly over a circle', xlab = 'x', ylab = 'y')
points(x = CENTRE[1], y = CENTRE[2], col = 'red', pch = 16)

Points on a ball

| cite | improve this answer | |
$\endgroup$
4
$\begingroup$

The simplest and least error-prone approach - for low dimensions (see below!) - would still be rejection sampling: pick uniformly distributed points from the $m$-dimensional hypercube circumscribing the sphere, then reject all that fall outside the ball.

runifball <- function(n, centre = 0, center = centre, radius = 1) {
  
  #Check inputs
  if (!missing(centre) && !missing(center)) {
  if (sum((centre - center)^2) < 1e-15) { 
                 warning("specify 'centre' or 'center' but not both") } else {
                    stop("Error: specify 'centre' or 'center' but not both") } }
  if (radius < 0) { stop("Error: radius must be non-negative") }

  n_to_generate <- 2^length(center)*gamma(length(center)/2+1)*n/pi^(length(center)/2) # see below
  
  original_sample_around_origin <- 
      matrix(replicate(length(center),runif(n_to_generate ,-radius,radius)),nrow=n_to_generate )
  index_to_keep <- rowSums(original_sample_around_origin^2)<radius^2
  original_sample_around_origin[index_to_keep,]+
      matrix(center,nrow=sum(index_to_keep),ncol=length(center),byrow=TRUE)
}

Here is an application for the $m=2$-dimensional disk:

#Generate points uniformly on a disk
set.seed(1)
n      <- 10^5
CENTRE <- c(5, 3)
RADIUS <- 3
UNIF   <- runifball(n, centre = CENTRE, radius = RADIUS)

#Plot the points
plot(UNIF, 
     col = rgb(0, 0, 0, 0.05), pch = 16, asp = 1,
     main = 'Points distributed uniformly over a circle', xlab = 'x', ylab = 'y')
points(x = CENTRE[1], y = CENTRE[2], col = 'red', pch = 16)

ball

Once again, we will need to originally generate more points, because we will reject some. Specifically, we expect to keep $\frac{\pi^\frac{m}{2}}{2^m\Gamma(\frac{m}{2}+1)}$, which is the ratio of the volume of the $m$-dimensional ball to the volume of the $m$-dimensional hypercube circumscribing it. So we can either start by generating $\frac{2^m\Gamma(\frac{m}{2}+1)n}{\pi^\frac{m}{2}}$ and expect to end up with $n$ points (this is the approach the code above takes), or just start generating until we have kept $n$.

In either case, the number of points we originally need to draw in the hypercube in order to (expect to) end up with a single point in the ball rises quickly with increasing dimensionality $m$:

number of points to generate originally

(Note the logarithmic vertical axis!)

m <- 2:20
plot(m,2^m*gamma(m/2+1)/pi^(m/2),type="o",pch=19,log="y",
    xlab="Dimension (m)")

This is just a consequence of the fact that for large $m$, most of the volume of the $m$-dimensional hypercube is in the corners, not in the center (where the ball is). So rejection sampling is likely only an option for low dimensions.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ +1 this answer seems to be a dumb thing to do, since the other answer does not need to throw away samples. However, it may perform much faster for small dimensions, and mostly it can be generalized more easily. (this shows the interesting discrepancy between a precise mathematical result and a straight forward pragmatic approach). $\endgroup$ – Sextus Empiricus Aug 6 at 11:30
  • $\begingroup$ @SextusEmpiricus: I also like how it's much harder to make an error in this approach, and much easier to find one you have made. Better to spend five minutes coding something that takes a day to run, rather than spend two days coding and bugfixing something that takes five minutes to run. (Depending on circumstances, of course.) $\endgroup$ – Stephan Kolassa Aug 6 at 13:34
  • $\begingroup$ Personally, in this particular case, I find the probability to make a systematic error not so large. But, I appreciate the idea (of making errors and using simple code/principles) in a wider context. Making errors occurs with me a lot because I tend to make something rigorous (but requiring two days bugfixing, although in this simple case not but point taken) $\endgroup$ – Sextus Empiricus Aug 6 at 14:28

Not the answer you're looking for? Browse other questions tagged or ask your own question.