16
$\begingroup$

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:

Sample_space_1

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

Sample_space_2

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?

$\endgroup$
  • 3
    $\begingroup$ I suggest that the title is edited to "in a disk" rather than "in a circle", which can be misinterpreted. $\endgroup$ – pglpm Aug 5 at 8:24
  • 1
    $\begingroup$ @pglpm: Good point - edited. $\endgroup$ – Ben Aug 5 at 12:52
  • 1
    $\begingroup$ By taking an angle you're splitting the circle into slices. The problem is that these slices are narrower at the centre, and this keeps being a problem no matter how thin the slices are. Then you select a radius, which will, with equal probability, put points into the thin part at the centre and the thicker parts further from the centre. So then the area in the centre will have more points compared to some area further from the centre. This applies to the first approach at least, I don't know R (or whatever) well enough to follow the second approach. $\endgroup$ – Bernhard Barker Aug 5 at 16:12
  • 2
    $\begingroup$ This is an extremely common problem, and/ has been asked/answered literally dozens of times on stackoverflow: google.com/… $\endgroup$ – BlueRaja - Danny Pflughoeft Aug 5 at 17:50
  • 2
    $\begingroup$ @jreese Thank you for the links. You might enjoy searching our site, too. Some of the hits are stats.stackexchange.com/questions/462820, stats.stackexchange.com/questions/310346, stats.stackexchange.com/questions/7977, and stats.stackexchange.com/questions/258303. $\endgroup$ – whuber Aug 5 at 18:36
21
$\begingroup$

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:

enter image description here

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

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.

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.

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

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)

enter image description here

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

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