I need to generate a uniformly random point within a circle of radius R.
I realize that by just picking a uniformly random angle in the interval [0 ... 2π), and uniformly random radius in the interval (0 ... R) I would end up with more points towards the center, since for two given radii, the points in the smaller radius will be closer to each other than for the points in the larger radius.
I found a blog entry on this over here but I don't understand his reasoning. I suppose it is correct, but I would really like to understand from where he gets (2/R2)×r and how he derives the final solution.
Update: 7 years after posting this question I still hadn't received a satisfactory answer on the actual question regarding the math behind the square root algorithm. So I spent a day writing an answer myself. Link to my answer.
Here is my Python code to generate
num
random points from a circle of radiusrad
:Let's approach this like Archimedes would have.
How can we generate a point uniformly in a triangle ABC, where |AB|=|BC|? Let's make this easier by extending to a parallelogram ABCD. It's easy to generate points uniformly in ABCD. We uniformly pick a random point X on AB and Y on BC and choose Z such that XBYZ is a parallelogram. To get a uniformly chosen point in the original triangle we just fold any points that appear in ADC back down to ABC along AC.
Now consider a circle. In the limit we can think of it as infinitely many isoceles triangles ABC with B at the origin and A and C on the circumference vanishingly close to each other. We can pick one of these triangles simply by picking an angle theta. So we now need to generate a distance from the center by picking a point in the sliver ABC. Again, extend to ABCD, where D is now twice the radius from the circle center.
Picking a random point in ABCD is easy using the above method. Pick a random point on AB. Uniformly pick a random point on BC. Ie. pick a pair of random numbers x and y uniformly on [0,R] giving distances from the center. Our triangle is a thin sliver so AB and BC are essentially parallel. So the point Z is simply a distance x+y from the origin. If x+y>R we fold back down.
Here's the complete algorithm for R=1. I hope you agree it's pretty simple. It uses trig, but you can give a guarantee on how long it'll take, and how many
random()
calls it needs, unlike rejection sampling.Here it is in Mathematica.
The area element in a circle is dA=rdr*dphi. That extra factor r destroyed your idea to randomly choose a r and phi. While phi is distributed flat, r is not, but flat in 1/r (i.e. you are more likely to hit the boundary than "the bull's eye").
So to generate points evenly distributed over the circle pick phi from a flat distribution and r from a 1/r distribution.
Alternatively use the Monte Carlo method proposed by Mehrdad.
EDIT
To pick a random r flat in 1/r you could pick a random x from the interval [1/R, infinity] and calculate r=1/x. r is then distributed flat in 1/r.
To calculate a random phi pick a random x from the interval [0, 1] and calculate phi=2*pi*x.
It really depends on what you mean by 'uniformly random'. This is a subtle point and you can read more about it on the wiki page here: http://en.wikipedia.org/wiki/Bertrand_paradox_%28probability%29, where the same problem, giving different interpretations to 'uniformly random' gives different answers!
Depending on how you choose the points, the distribution could vary, even though they are uniformly random in some sense.
It seems like the blog entry is trying to make it uniformly random in the following sense: If you take a sub-circle of the circle, with the same center, then the probability that the point falls in that region is proportional to the area of the region. That, I believe, is attempting to follow the now standard interpretation of 'uniformly random' for 2D regions with areas defined on them: probability of a point falling in any region (with area well defined) is proportional to the area of that region.
There is a linear relationship between the radius and the number of points "near" that radius, so he needs to use a radius distribution that is also makes the number of data points near a radius
r
proportional tor
.Think about it this way. If you have a rectangle where one axis is radius and one is angle, and you take the points inside this rectangle that are near radius 0. These will all fall very close to the origin (that is close together on the circle.) However, the points near radius R, these will all fall near the edge of the circle (that is, far apart from each other.)
This might give you some idea of why you are getting this behavior.
The factor that's derived on that link tells you how much corresponding area in the rectangle needs to be adjusted to not depend on the radius once it's mapped to the circle.
Edit: So what he writes in the link you share is, "That’s easy enough to do by calculating the inverse of the cumulative distribution, and we get for r:".
The basic premise is here that you can create a variable with a desired distribution from a uniform by mapping the uniform by the inverse function of the cumulative distribution function of the desired probability density function. Why? Just take it for granted for now, but this is a fact.
Here's my somehwat intuitive explanation of the math. The density function f(r) with respect to r has to be proportional to r itself. Understanding this fact is part of any basic calculus books. See sections on polar area elements. Some other posters have mentioned this.
So we'll call it f(r) = C*r;
This turns out to be most of the work. Now, since f(r) should be a probability density, you can easily see that by integrating f(r) over the interval (0,R) you get that C = 2/R^2 (this is an exercise for the reader.)
Thus, f(r) = 2*r/R^2
OK, so that's how you get the formula in the link.
Then, the final part is going from the uniform random variable u in (0,1) you must map by the inverse function of the cumulative distribution function from this desired density f(r). To understand why this is the case you need to find an advanced probability text like Papoulis probably (or derive it yourself.)
Integrating f(r) you get F(r) = r^2/R^2
To find the inverse function of this you set u = r^2/R^2 and then solve for r, which gives you r = R * sqrt(u)
This totally makes sense intuitively too, u = 0 should map to r = 0. Also, u = 1 shoudl map to r = R. Also, it goes by the square root function, which makes sense and matches the link.