I need this seldom enough to forget how it's done - but then it's annoying to have to think/google for the solution again when I do need it... So I'll document here.
The task is to generate uniformly distributed numbers within a circle of radius R in the (x,y) plane. At first polar coordinates seems like a great idea, and the naive solution is to pick a radius r uniformly distributed in [0, R], and then an angle theta uniformly distributed in [0, 2pi]. BUT, you end up with an exess of points near the origin (0, 0)! This is wrong because if we look at a certain angle interval, say [theta, theta+dtheta], there needs to be more points generated further out (at large r), than close to zero. The radius must not be picked from a uniform distribution, but one that goes as
pdf_r = (2/R^2)*r
That's easy enough to do by calculating the inverse of the cumulative distribution, and we get for r:
r = R*sqrt( rand() )
where rand() is a uniform random number in [0, 1]. Here is a picture:
some matlab code here.
The thinking for generating random points on the surface of a sphere in 3D is very similar. If I get inspired I will do a post on that later, meanwhile you can go read these lecture notes.
2 Comments
Why is pdf_r= (2/R^2)*r ?
In polar coordinates an infinitesimal area-element is: dr * r * dtheta
if we want points uniformly distributed in this area the pdf for r should "clearly" be proportional to r.
we also want the pdf to be normalized so that the integral from 0 to R of the pdf is 1. That leads to the pre-factor (2/R^2)
Thanks. I too made the mistake but realized later that it would be more concentrated closer to the origin.
Your code could be a lot more efficient by vectorizing instead of looping:
% prepare
Nmax = 1e3;
R = 5;
% Wrong
r1 = R*rand(Nmax,1);
theta1 = 2*pi*rand(Nmax,1);
x1 = r1.*cos(theta1);
y1 = r1.*sin(theta1);
% Right
r2 = R*sqrt(rand(Nmax,1));
theta2 = 2*pi*rand(Nmax,1);
x2 = r2.*cos(theta2);
y2 = r2.*sin(theta2);
% plot it
subplot(1,2,1);
plot(x1, y1, 'r.');
axis square
title('Left')
subplot(1,2,2);
plot(x2, y2, 'g.');
axis square
title('Right')
Post a Comment