How can I create randomly scattered points between two circles?

How can I create randomly scattered points between two circles?

I would like to have 3 inputs (assuming that the two circles are concentric, with center at (0,0)):

1. the number of scattered points
2. the radius of the inner circle
3. the radius of the outer circle

% create two concentric circles with radii r=1 and r=2
figure('Color',[1 1 1]);
for r = 1:2
th = 0:pi/50:2*pi;
xc = r * cos(th);
yc = r * sin(th);
hold on
plot(xc,yc,'color','r','LineWidth',2)
end
axis equal
set(gca,'Visible','off')

[Sim]

I am not sure if this is correct:

% create the concentric circles
figure('Color',[1 1 1]);
for r = 1:2
th = 0:pi/50:2*pi;
xc = r * cos(th);
yc = r * sin(th);
hold on
plot(xc,yc,'color','r','LineWidth',1)
end

% generate N random numbers in the interval (a,b) with the formula r = a + (b-a).*rand(N,1)
N=500;
t = 2 * pi * rand(N,1);
g = 1 + (2-1).*rand(N,1);
x = g.*cos(t);
y = g.*sin(t);
plot(x,y,'o','MarkerFaceColor','b','MarkerEdgeColor','b')

axis equal
set(gca,'Visible','off')

Answer 1

[Adam Danz]

Your solution is fine.

a = 4;
b = 10;
N = 1000;

r = a + (b-a) * rand(N,1);
th = 2*pi * rand(N,1); 
x = r.*cos(th);
y = r.*sin(th);

where a<=b.

Another way to think about it is in polar coordinates

[x, y] = pol2cart(th, r);

Update

As @John D'Errico points out below, the solution above will not produce a uniform distribution of coordinates in 2D space. This is because points will be compressed as the radii decreases within the disc. To compensate for this, you could use inverse transform sampling to generate random r values from a nonuniform distribution that accounts for the increasing circumference size.

r2 = a + (b-a) * rand(N,1).^.5; 
x2 = r2.*cos(th);
y2 = r2.*sin(th);

Compare these results

The first method clusters points near the inner radius. The second method is closer to uniform in 2D space but has a negative bias arond the inner radius. See comments below for more on this topic.

figure()
tiledlayout(2,2)
plot(nexttile(),x,y,'b.')
title('Uniform th and r sampling')
axis equal tight
grid on
rectangle('Position',[-a -a 2*a 2*a],'Curvature',1,'EdgeColor','r')
rectangle('Position',[-b -b 2*b 2*b],'Curvature',1,'EdgeColor','r')

plot(nexttile(),x2,y2,'b.')
title('Inverse transform sampling')
axis equal tight
grid on
rectangle('Position',[-a -a 2*a 2*a],'Curvature',1,'EdgeColor','r')
rectangle('Position',[-b -b 2*b 2*b],'Curvature',1,'EdgeColor','r')

N = 100000; 
r = a + (b-a) * rand(N,1);
th = 2*pi * rand(N,1); 
x = r.*cos(th);
y = r.*sin(th);

r2 = a + (b-a) * rand(N,1).^.5; 
x2 = r2.*cos(th);
y2 = r2.*sin(th);

histogram2(nexttile(),x,y,20,'DisplayStyle','tile','EdgeColor','none')
colorbar
axis equal 
histogram2(nexttile(),x2,y2,20,'DisplayStyle','tile','EdgeColor','none')
colorbar
axis equal

[Sim]

Thanks a lot!! :-)

[John D'Errico]

Um, be careful!

This does not generate UNIFORMLY random points between two circles! Yes, it does generate random points. But they will be more densely populated closer to the center than near the outer radius. For example...

a = 1;
b = 10;
N = 5000;
r = a + (b-a) * rand(N,1);
th = 2*pi * rand(N,1); 
x = r.*cos(th);
y = r.*sin(th);

plot(x,y,'.')
axis equal

You should see that the points near the inner circle are far more dense. Instead, you need to more carefully sample in radius. The uniform sampling in r chosen by @Adam Danz is not correct. Instead, you want an increasing fraction for large values of the radius if your goal is a uniform sample.

[Adam Danz]

The cartesian coordinates are not uniform but the polar coordinates are uniform. This is, of course, due to a compression of space around the circumference of a circle as the radius shrinks.

a = 1;
b = 10;
N = 5000;
r = a + (b-a) * rand(N,1);
th = 2*pi * rand(N,1); 
x = r.*cos(th);
y = r.*sin(th);

tiledlayout(2,2)
histogram(nexttile(),r,50)
title('distribution of r')

histogram(nexttile(),th,50)
title('distribution of th')

histogram(nexttile(),x,50)
title('distribution of x')

histogram(nexttile(),y,50)
title('distribution of y')

John's probably right that the intent is to create a uniform distribution in cartesian coordinates which I hadn't thought of. One direction you could take is to produce random r values from a distribution that considers the linearly increasing circumference size. This is known as inverse transform sampling.

a = 1;
b = 10;
N = 5000;
r = a + (b-a) * rand(N,1).^.5; % <-- updated
th = 2*pi * rand(N,1); 
x = r.*cos(th);
y = r.*sin(th);

figure()
plot(x,y,'.')
axis equal
grid on

Now the x and y distributions are closer to uniform. Note the falloff at the edges of the x and y distributions. This is expected since area of the disc within each 1D bin decreases as the bins get closer to the edges of the disc. Also note the dips in the x and y distributions around 0. This is, of course, due to the hold in the middle of the disc which doesn't contain points.

figure()
tiledlayout(2,2)
histogram(nexttile(),r,50)
title('distribution of r')

histogram(nexttile(),th,50)
title('distribution of th')

histogram(nexttile(),x,50)
title('distribution of x')

histogram(nexttile(),y,50)
title('distribution of y')

I'll update my answer to include this method.

[Sim]

Thanks a lot @John D'Errico for your post and I will bear this in mind! In case you have a short code that might show an uniform random scatter of points in the annulus, I might change the accepted answer (if everyone is OK with that).

[Paul]

Hi Sim,

The term "uniform random scatter of points" is ambiguous.

At a minimum, the problem statement needs to define the random variables. Then, based on those random variables, the problem statement needs to define the probability density function, either before or after conditioning.

For example, we could say

a) Let R and Theta be independent random variables with a uniform joint density function f(r,theta) that is constant over a < r < b , 0 < theta < 2*pi and zero otherwise. Generate samples of R and Theta from this density. This definition is assumed in the first solution of @Adam Danz.

b = 10;
a = 5;
N = 1e7;
rng(100)

r = a + (b-a) * rand(N,1);
th = 2*pi * rand(N,1); 
x = r.*cos(th);
y = r.*sin(th);

figure
histogram2(x,y,40,'Normalization','pdf','DisplayStyle','tile','EdgeColor','none')
xlabel('x'),ylabel('y')
title('f(x,y | a^2 < x^2 + y^2 < b^2)')
colorbar

As has already been shown, the density in x-y space is not constant over the annulus, even though the joint density of R and Theta is uniform (and R and Theta are uniform)

or

b) Let X and Y be independent random variables with uniform joint density function f(x,y) that is constant over -b < x < b, -b < y < b and zero otherwise. With these definitions, X and Y both have marginal uniform densities, i.e., viewed in isolation X is uniform and so is Y. Let f(x,y | a^ < x^2 + y^2 < b^2) be the joint density function of x and y conditioned on (x,y) lying in the annulus. Generate samples from this conditional density.

clearvars -except a b N

x = -b + 2*b*rand(N,1);
y = -b + 2*b*rand(N,1);
r = sqrt(x.^2 + y.^2);
keep = r <= b & r >= a;
x = x(keep);
y = y(keep);

Over the annulus, the joint conditional density should be constant with a value of

1/(pi*(b^2-a^2))

ans = 0.0042

figure
histogram2(x,y,40,'Normalization','pdf','DisplayStyle','tile','EdgeColor','none')
colorbar
xlabel('x'),ylabel('y')
title('f(x,y | a^2 < x^2 + y^2 < b^2)')

Here, the conditional joint density is uniform in x-y space over the annulus (I think the small edges are vagaries of histogram2?) However, the marginal densities of X and Y after condiitoning are NOT uniform (and look like something from Lord of the Rings!)

figure
subplot(121)
histogram(x,'Normalization','pdf','EdgeColor','none');
xlabel('x'),ylabel('f_X(x | a^2 < x^2 + y^2 < b^2)')
subplot(122)
histogram(y,'Normalization','pdf','EdgeColor','none');
xlabel('y'),ylabel('f_Y(y | a^2 < x^2 + y^2 < b^2)')

c) For completeness and comparison, here is the method in the Update of this answer.

clearvars -except a b N

th = 2*pi * rand(N,1); 
r2 = a + (b-a) * rand(N,1).^.5; 
x2 = r2.*cos(th);
y2 = r2.*sin(th);

figure
histogram2(x2,y2,40,'Normalization','pdf','DisplayStyle','tile','EdgeColor','none')
colorbar
xlabel('x'),ylabel('y')
title('f(x,y | a^2 < x^2 + y^2 < b^2)')

figure
subplot(121)
histogram(x2,'Normalization','pdf','EdgeColor','none');
xlabel('x'),ylabel('f_X(x | a^2 < x^2 + y^2 < b^2)')
subplot(122)
histogram(y2,'Normalization','pdf','EdgeColor','none');
xlabel('y'),ylabel('f_Y(y | a^2 < x^2 + y^2 < b^2)')

I'm not quite sure in what sense this solution is 'uniform.' Perhaps @Adam Danz can comment.

[Torsten]

https://stackoverflow.com/questions/13064912/generate-a-uniformly-random-point-within-an-annulus-ring

[Paul]

Trying the methods from @Torsten's link.

d) use the inverse sampling method. Maybe this is what the Updated method in the answer was intended to be? Seems to match option (b) above.

b = 10;
a = 5;
N = 1e7;
rng(100)

th = 2*pi*rand(N,1);
r  = sqrt(rand(N,1)*(b^2-a^2)+a^2);
x = r.*cos(th);
y = r.*sin(th);

figure
histogram2(x,y,40,'Normalization','pdf','DisplayStyle','tile','EdgeColor','none')
colorbar
xlabel('x'),ylabel('y')
title('f(x,y | a^2 < x^2 + y^2 < b^2)')

figure
subplot(121)
histogram(x,'Normalization','pdf','EdgeColor','none');
xlabel('x'),ylabel('f_X(x | a^2 < x^2 + y^2 < b^2)')
subplot(122)
histogram(y,'Normalization','pdf','EdgeColor','none');
xlabel('y'),ylabel('f_Y(y | a^2 < x^2 + y^2 < b^2)')

e) I'm not sure what to call this method. Seems to match (b) and (d).

clearvars -except a b N

u = rand(N,1)*(b-a) + a;
v = rand(N,1)*(b+a);
r = nan*u; tf = v < u;
r(tf)  = u(tf);
r(~tf) = a + b - u(~tf);
th = rand(N,1)*2*pi;

x = r.*cos(th);
y = r.*sin(th);

figure
histogram2(x,y,40,'Normalization','pdf','DisplayStyle','tile','EdgeColor','none')
colorbar
xlabel('x'),ylabel('y')
title('f(x,y | a^2 < x^2 + y^2 < b^2)')

figure
subplot(121)
histogram(x,'Normalization','pdf','EdgeColor','none');
xlabel('x'),ylabel('f_X(x | a^2 < x^2 + y^2 < b^2)')
subplot(122)
histogram(y,'Normalization','pdf','EdgeColor','none');
xlabel('y'),ylabel('f_Y(y | a^2 < x^2 + y^2 < b^2)')

[David Goodmanson]

( this comment moved to an answer)

[Adam Danz]

Great discusion here. @Paul and @David Goodmanson, please consider moving your solutions to a new answer to increase visibility. I can't un-accept my answer at the moment and there are much bettern approaches in this thread.

[Sim]

Many many thanks to everyone!

I did not expect that such a question could generate so many interesting comments/discussions/answers!

In these days I am not able to read carefully all your answers, but I will do it as soon as possible!

Meanwhile: have a great start into the New Year!

P.S. I thought to temporarily unaccept the @Adam Danz's answer (and please apologise Adam!! Your answer is really great!) just to encourage other replies and give more visibility to everyone...

Answer 2

[David Goodmanson]

(At Adam's suggestion I moved my previous comment to this answer)

What tells the tale is the differential area element in the (r, theta) plane, which is r dr dtheta.

Picking theta from a uniform random variable on 0<theta<2pi and letting x = r*cos(theta), y = r*sin(theta) takes care of the angle part in the usual way. That leaves picking r.

To fill the plane uniformly in r dr with draws from a uniformly distributed variable rho, then those two quantities must be proportonal, so

r dr = c drho                                            ***

for some constant c. Everything else follows from that expression.

Now

r dr = d(r^2/2) = c drho    and consequently
r^2/2 = c rho + d                                        (4)

for a couple of constants c and d. For an annulus of radii a and b,

r = a  @  rho = 0  and
r = b  @  rho = 1

which fixes c and d, and

r^2 = a^2 + (b^2-a^2) rho       so
r = sqrt(a^2 + (b^2-a^2) rho)                             (3)

Then for the choosing of N random points,

r = sqrt(a^2 + (b^2-a^2)*rand(1,N))                       ***

which is the correct result and is one of the previously suggested solutions. What does not work are two other previous suggestions,

r = a + (b-a) sqrt(rho)        % Incorrect                (2)
r^2/2 = (1/2) ( a^2 + 2a(b-a)sqrt(rho) + (b-a)^2 rho )

r = a + (b-a) rho               % Incorrect               (1)
r^2/2 = (1/2) (a^2 + 2a(b-a) rho + (b-a)^2 rho^2 )

neither of which is of the required form (4).

It's interesting to take a look at the point densities that are produced by these three cases. The density is, for N points (assuming angular symmetry),

sigma = dn/dA = N drho / (2pi r dr) = N / ( 2pi r dr/drho )  

Ignoring the N/2pi factor for comparison purposes,then

rho = 0:.001:1;
a = 1;
b = 2;
% case 1
r1 = a+(b-a)*rho;
drdrho1 = (b-a);
sigma1 = 1./(r1.*drdrho1);
% case 2
r2 =  a+(b-a)*sqrt(rho);
drdrho2 = (1/2)*(b-a)./sqrt(rho);
sigma2 = 1./(r2.*drdrho2);
% case 3
r3 = sqrt(a^2+(b^2-a^2)*rho);
drdrho3 = (1/2)./sqrt(a^2+(b^2-a^2)*rho) * (b^2-a^2);  
sigma3 = 1./(r3.*drdrho3);
plot(r1,sigma1,r2,sigma2,r3,sigma3)
xlabel('r')
ylabel('density' )
legend('[1]  a+(b-a)*rho','[2]  a+(b-a)*sqrt(rho)', ...
    '[3]  sqrt(a^2+(b^2-a^2)*rho)','location','south')

As has been previously mentioned, case 1 throws too many points to the inside, and case2 throws too many points to the outside. Case 3 is just right.

oo As a side note, N points go into a specified area. Replacing the draws from the two random uniform variables with integrals over rho and theta, then for case (3), and using 2pi for the dtheta integral due to angular symmetry

r = sqrt(a^2 + (b^2-a^2) rho)
dr = (1/2) (b^2-a^2) 1/sqrt(a^2 + (b^2-a^2) rho) drho
rdr = (1/2) (b^2-a^2) drho
Integral (rdr dtheta) = (1/2) (b^2-a^2) (2pi) = pi (b^2-a^2)                                 % annulus

which is the correct area. That's not a ringing endorsement of (3), though, since *any* well behaved function r = f(rho) such that f(0) = a and f(1) = b (including (1) and (2) ), will produce the same result, albeit with an incorrect density of points.

[John D'Errico]

I think the incorrect form stems from an intuitive leap, because when a==0, so the circle will be completely filled, then this

r = a + (b-a)*sqrt(rho)

reduces to

r = b*sqrt(rho)

And that form should fill the circle in a uniform way.

As well, the generation of random points across a triangle uses a similar form with the sqrt. So this is a not uncommon idea. It just does not work for an annulus with a>0.

[Sim]

@Adam Danz, @John D'Errico, @Paul, @Torsten many many thanks to all of you for all your answers and comments, that I found extremely interesting and insightful! I have really appreciated them!

.....and now, which one to accept ? This is quite embarassing, since all of your answers are top-notch solutions, and they all deserve to be accepted and used by the community! :-) :-) :-)

Answer 3

[John D'Errico]

Let me explain it differently than the rest, with a simple, intuiitive derivation.

Again, the problem is perfectly symmetric in polar angle, so that part is trivial. All we need to understand is the distribution in the radius, r.

Assume you wish to see a uniform sampling on the annulus from a to b in radius. Now consider any fixed radius r, so now just a circle of radius r. We want to see proportionally 2*pi*r*d points around that infinitessimally narrow circular annulus. Essentially, this gives us the radial distribution. The PDF is simply:

syms r d a b
PDF = 2*pi*r*d;

Here d is the packing density, and r will vary from a to b. We need to normalize the pdf so it has unit area on that domain.

PDF = PDF/int(PDF,r,[a,b])

PDF = 

I hope that makes sense, in a very intuitive way. Now that we have the PDF, we can compute the CDF of the radial distribution, as simply the integral over r.

syms R
CDF = int(PDF,r,[a,R])

CDF = 

Again, r will live on the interval [a,b]. And also see that I used the definite integral from a to the radius R. This is important!

For example, with a=1, b=10, we would have:

fplot(subs(PDF,[a,b],[1,10]),[1,10])

fplot(subs(CDF,[a,b],[1,10]),[1,10])

As you can see, the PDF is linear, and the CDF is a simple quadratic polynomial.

How do we sample from a distribution when the CDF is known, and it has a simple inverse? This itself is simple enough. We use the inverse of the CDF. That is, we choose a UNIFORM random sampling on the interval [0,1]. rand does this perfectly. Then we solve for R.

syms u
CDFinv = solve(CDF == u,R)

CDFinv = 

We clearly want the second solution there, since it is the positive one.

CDFinvfun = @(u,a,b) sqrt((b^2-a^2)*u + a^2)

CDFinvfun = function_handle with value:
    @(u,a,b)sqrt((b^2-a^2)*u+a^2)

I just wrote it as a function, since matlabFunction did something screwy to me, and it is too early in the AM right now for me to think.

N = 1000000;
Usample = rand(N,1);
Rsample = CDFinvfun(Usample,1,10);
histogram(Rsample,100,'norm','pdf')

And you see a perfectly linear sample PDF drop out ion the interval [a,b]. Now lets try it out.

N = 10000;
Usample = rand(N,1);
Rsample = CDFinvfun(Usample,1,10);
theta = rand(N,1)*2*pi;
x = Rsample.*cos(theta);
y = Rsample.*sin(theta);

plot(x,y,'.')
axis equal

And that is exactly what we want to see. A nice uniform sampling on the annulus.

If what I did does not make sense, then go through it step by step. I tried to be as complete as possible in my Answer, but I am sure I went too fast at some point.

Addendum: Some additional points of interest I may have glossed over above,

1. The entire solution stems directly from the assumption of a uniform scattering on the annulus. That directly implies the PDF, and then the rest follows.
2. The use of a definite integral to compute the CDF is important, since if you just compute an indefinite integral with no other flags to int, the constant of integration will be lost. The definite integral implicitly builds the (correct) constant of integration into the CDF. I think this is an easily missed point.
3. The above analysis can also trivially be extended to 3 or more dimensions, so with annular spheres. There the transformation will be different, but as easily derived based on the resulting polynomial PDF on the finite annular domain.
4. If you follow the solution I've posed, you can see that whan a==0, and we want to simply sample uniformly over the entire circle of radius b, then the intuitive square root transformation that Adam proposed now appears. But that clearly fails when a>0. When a==0, the inverse CDF reduces to simply b*sqrt(u),

[Paul]

Hi John,

Just came to point out that finverse is another option to get the inverse of the CDF

syms R a b u
CDF = -(R^2-a^2)/(a^2-b^2)

CDF = 

CDFinv = subs(finverse(CDF,R),R,u)

iCDF =