试用软件

How can I create randomly scattered points between two circles?

24 次查看（过去 30 天）

显示更早的评论

Sim 2023-12-28

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2064826-how-can-i-create-randomly-scattered-points-between-two-circles

编辑： David Goodmanson 2024-1-17

在 MATLAB Online 中打开

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)):

the number of scattered points
the radius of the inner circle
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')

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Sim 2023-12-28

编辑：Sim 2023-12-28

在 MATLAB Online 中打开

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')

请先登录，再进行评论。

请先登录，再回答此问题。

回答（3 个）

Adam Danz 2023-12-28

3
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2064826-how-can-i-create-randomly-scattered-points-between-two-circles#answer_1379611

编辑：Adam Danz 2023-12-29

在 MATLAB Online 中打开

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

10 个评论
显示 8更早的评论隐藏 8更早的评论

Adam Danz 2023-12-29

在 MATLAB Online 中打开

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.

Paul 2023-12-29

编辑：Paul 2023-12-29

在 MATLAB Online 中打开

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.

Paul 2023-12-30

在 MATLAB Online 中打开

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)')

Adam Danz 2023-12-30

编辑：Adam Danz 2023-12-30

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 2023-12-30

编辑：Sim 2023-12-30

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...

请先登录，再进行评论。

David Goodmanson 2023-12-30

3
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2064826-how-can-i-create-randomly-scattered-points-between-two-circles#answer_1380646

编辑：David Goodmanson 2024-1-17

在 MATLAB Online 中打开

(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.

2 个评论
显示无隐藏无

John D'Errico 2023-12-30

编辑：John D'Errico 2023-12-30

在 MATLAB Online 中打开

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 2024-1-16

编辑：Sim 2024-1-16

@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! :-) :-) :-)

请先登录，再进行评论。

John D'Errico 2023-12-30

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2064826-how-can-i-create-randomly-scattered-points-between-two-circles#answer_1380456

编辑：John D'Errico 2023-12-30

在 MATLAB Online 中打开

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,

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.
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.
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.
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),

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Paul 2023-12-30

在 MATLAB Online 中打开

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 =

请先登录，再进行评论。

请先登录，再回答此问题。

类别

MATLAB Graphics Formatting and Annotation Axes Appearance Combine Multiple Plots Subplots

在 Help Center 和 File Exchange 中查找有关 Subplots 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Translated by

试用软件