Generating random points uniformly on a disk [duplicate]

Question

I have to randomly generate 1000 points over a unit disk such that are uniformly distributed on this disk. Now, for that, I select a radius $r$ and angular orientation $\alpha$ such that the radius $r$ is a uniformly distributed variate from $r \in [0,1]$ while $\alpha$ is a uniformly distributed variate from $\alpha \in [0, 2\pi]$ using the following code

r <- runif(1000, min=0, max=1) 
alpha <- runif(1000, min=0, max=2*pi)
x <- r*cos(alpha)
y <- r*sin(alpha)
plot(x,y, pch=19, col=rgb(0,0,0,0.05), asp=1)

Then I look at my sample space and it looks like this:

This obviously doesn't look like a sample with uniform distribution over the disk. Hence, I guessed that the problem might be occurring as a result of a lack of independence between the variables $r$ and $\alpha$ in contingency to how they've been linked computationally.

To take care of that I wrote a new code.

rm(list=ls())
r <- runif(32, min=0, max=1)
df_res <- data.frame(matrix(c(-Inf, Inf), byrow = T, nrow = 1))
for (i in 1:32) {
  for (j in 1:32) {
    alpha <- runif(32, min=0, max=2*pi)
    r <- runif(32, min=0, max=1)
    df <- data.frame(matrix(c(r[i],alpha[j]), byrow = T, nrow = 1))
    df_res <- rbind(df_res,df)
  }
}
df_res <- subset(df_res, df_res$X1 != -Inf)
x<- df_res$X1 *cos(df_res$X2)
y <- df_res$X1 *sin(df_res$X2)
plot(x,y, pch=19, col=rgb(0,0,0,0.05), asp=1)

And, yet again the sample looks non-uniformly distributed over the disk

I'm starting to suspect that there is a deeper mathematical problem going on in the vicinity. Could someone help me write code that would create a sample space uniformly distributed over the disk, or explain the mathematical fallacy if any in my reasoning?

Answer 1

The problem is due to the fact that the radius is not uniformly distributed. Namely, if $(X,Y)$ is uniformly distributed over $$\left\{ (x,y);\ x^2+y^2\le 1\right\}$$ then the (polar coordinates) change of variables $$R=(X^2+Y^2)^{1/2}\qquad A=\text{sign}(Y)\arccos(X/R)$$ has the density $$\frac{1}{\pi} \mathbb{I}_{(0,1)}(r)\left|\frac{\text{d}(X,Y)}{\text{d}(R,A)}(r,\alpha)\right|\mathbb{I}_{(0,2\pi)}(\alpha)$$ Using $x = r \cos \alpha$ and $y = r \sin \alpha$ leads to $$\left|\frac{\text{d}(X,Y)}{\text{d}(R,A)}(r,\alpha)\right|=r(\sin^2\alpha+\cos^2\alpha)=r$$ Therefore, the angle $A$ is distributed uniformly over $(0,2\pi)$ but the radius $R$ has density $f(r)=2r\mathbb{I}_{(0,1)}(r)$ and cdf $F(r)=r^2$ over $(0,1)$. As one can check by running

r <- sqrt(runif(1000, min=0, max=1) )
alpha <- runif(1000, min=0, max=2*pi)
x <- r*cos(alpha)
y <- r*sin(alpha)
plot(x,y, pch=19, col=rgb(0,0,0,0.05), asp=1)

where the radius is simulated by the inverse cdf representation, which makes it the square root of a Uniform variate, the random repartition of the 10³ simulated points is compatible with a uniform:

Answer 2

The simplest and least error-prone approach would be rejection sampling: generate uniformly distributed points in the square around your circle, and only keep those that are in the circle.

nn <- 1e4
radius <- 1
set.seed(1) # for reproducibility
foo <- cbind(runif(nn,-radius,radius),runif(nn,-radius,radius))
plot(foo[rowSums(foo^2)<radius^2,],pch=19,cex=0.6,xlab="x",ylab="y")

Of course, you will only keep a fraction of your generated data points, around $\frac{\pi}{4}$ (which is the ratio of the areas of the circumscribed square to the disk). So you can either start with $\frac{4n}{\pi}$ points, or generate points until you keep your target number $n$ of them.

Answer 3

You can find the mathematics of this situation in a related question here. The method is set out in Xi'an's excellent answer, and it can be summarised by the following requirements:

$$\begin{matrix} R^2 \sim \text{U}(0,1) \quad \ \ & & & X = R \cos (\theta), \\[6pt] \theta \sim \text{U}(0, 2\pi) & & & Y = R \sin(\theta). \\[6pt] \end{matrix}$$

Following on from that other answer, when you come up with these solutions, it is often useful to try to generalise them into functions that can generate random values for a particular class of problems. A natural generalisation in this case is to look at randomly generated points on a circle with an arbitrary centre and radius. Using the same basic method as in the existing answer, here is a general function to produce random points uniformly over a circle with arbitrary centre and radius.

runifcircle <- function(n, centre = c(0, 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
  OUT      <- matrix(0, nrow = 2, ncol = n)
  rownames(OUT) <- c('x', 'y')
  
  #Generate uniform values on circle
  r2       <- runif(n, min = 0, max = radius^2)
  theta    <- runif(n, min = 0, max = 2*pi)
  OUT[1, ] <- center[1] + sqrt(r2)*cos(theta)
  OUT[2, ] <- center[2] + sqrt(r2)*sin(theta)
  
  OUT }

Creating this function allows you to easily generate any number of points over an arbitrary circle. (If you want an interesting exercise that extends this problem, try modifying the above function to create a new function runifball that generates uniform random values on a hypersphere with an arbitrary centre and radius.) We can easily verify that this function works correctly by plotting the results for a large number of sample values.

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

#Plot the points
plot(UNIF[1, ], UNIF[2, ], 
     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)