Simulate a uniform distribution on a disc I was attempting to simulate injection of random points within a circle, such that any part of the circle has the same probability of having a defect.  I expected the count per area of the resulting distribution to follow a Poisson distribution if I break up the circle into equal area rectangles.  
Since it requires only placing points within a circular area, I injected two uniform random distributions in polar coordinates: $R$ (radius) and $\theta$ (polar angle).
But after doing this injection, I clearly get more points in the center of the circle compared to the edge. 

What would be the correct way to perform this injection across the circle such that the points are randomly distributed across the cirlce?
 A: Rejection Sampling can be used. This means we can sample from 2D uniform distribution, and select samples that satisfy the disc condition.
Here is an example. 
x=runif(1e4,-1,1)
y=runif(1e4,-1,1)

d=data.frame(x=x,y=y)
disc_sample=d[d$x^2+d$y^2<1,]
plot(disc_sample)


A: I'll give you a general n-dimensional answer that works for two-dimensional case too, of course. In three dimensions an analog of a disk is a volume of a solid ball (sphere).
There are two approaches I'm going to discuss. One of them I would call "precise", and you'll get a complete solution with it in R. The second one I call heuristic, and it's only the idea, no complete solution is provided.
"Precise" solution
My solution is based on Marsaglia and Muller's works. Basically, it happens so that the Gaussian vector normalized to its norm would give you the uniformly distributed points on a d-dimensional hypersphere:

It's the same as uniformely distributed points on a circle in two dimensions. To extend this to the entire surface of a disk, you need to further scale them by radius. The square of a radius is from uniform distribution in two dimensions, or raised to power $d$ in d-dimensions. So, you raise to power $1/d$ a uniform random number to get the properly distributed radius. Here's a complete code  in R for two dimensions, which you can easily extend to any number of dimensions:
n <- 1e4
rho <- sqrt(runif(n))
# d - # of dimensions of hyperdisk
d = 2
r = matrix(rnorm(n*d),nrow=n,ncol=d)
x = r/rep(sqrt(rowSums(r^2))/rho,1)
plot(x[,1], x[,2], pch=19, cex=0.6, col="#00000020")


Here's a code snippet for 3d case, i.e. a solid ball:
library(scatterplot3d)
n <- 1e3
# d - # of dimensions of hyperdisk

d=3
rho <- (runif(n))^(1/d)
r = matrix(rnorm(n*d),nrow=n,ncol=d)
x = r/rep(sqrt(rowSums(r^2))/rho,1)

scatterplot3d(x[,1], x[,2], x[,3])


Heuristic approach
This approach is based on a not so obvious fact that the ration of the volume of the unit hypersphere over the volume of a unit hypercube that encloses it shrinks to zero when the number of dimensions increases to infinity. This can easily be seen from the expression for a volume of a hypersphere:
$$V_n(R) = \frac{\pi^\frac{n}{2}}{\Gamma\left(\frac{n}{2} + 1\right)}R^n$$
Here, you can see how the coefficient in front of $R^n$ quickly decreases to zero. This is another manifestation of the phenomenon that is linked to what's known as a dimensionality curse in machine learning.
Why is this relevant to our problem at hand? Suppose, you want to generate d random uniform numbers, these would the random points inside d-dimensional hypercube. Next, you apply rejection sampling to pick the points inside the hypersphere (aka n-ball): $\sum_{i=1}^d x_i^2<R^2$. The problem is that for high number of dimensions d, almost all points will be outside the sphere! You'll end up throwing out vast majority of your samples.
The solution I propose is to use the rejection sampling with oversampling the points near the center. It turns out that if you were watching one of the cartesian coordinates of the random uniform sample from inside the ball, its distribution would be converging to a Gaussian with variance $\frac 1 {\sqrt{d+2}}$. So, instead of picking points uniformly from the cube, we be sampling the cartesian coordinate using the Gaussian, then apply rejection sampling on them. This way we would not be wasting as many generated random variates. This would be a form of importance sampling technique.
A: You want the proportion of points to be uniformly proportional to area rather than distance to the origin.  Since area is proportional to the squared distance, generate uniform random areas and take their square roots; scale the results as desired. Combine that with a uniform polar angle.
This is quick and simple to code, efficient in execution (especially on a parallel platform), and generates exactly the prescribed number of points.
Example
This is working R code to illustrate the algorithm.
n <- 1e4
rho <- sqrt(runif(n))
theta <- runif(n, 0, 2*pi)
x <- rho * cos(theta)
y <- rho * sin(theta)
plot(x, y, pch=19, cex=0.6, col="#00000020")


A: Here is an alternative solution in R:
n <- 1e4
## r <- seq(0, 1, by=1/1000)
r <- runif(n)
rho <- sample(r, size=n, replace=T, prob=r)
theta <- runif(n, 0, 2*pi)
x <- rho * cos(theta)
y <- rho * sin(theta)
plot(x, y, pch=19, cex=0.6, col="#00000020")


