Generate random points inside a circle

7 次查看（过去 30 天）
Mahesh Ramaraj
Mahesh Ramaraj 2011-1-23
编辑： Torsten 2022-5-8
I am trying to generate random points inside a circle given the radius and the center. If anyone has a sample code or can help me with this, thanks for the help.

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

采纳的回答

Paulo Silva
Paulo Silva 2011-1-23
function [x y]=cirrdnPJ(x1,y1,rc)
%the function, must be on a folder in matlab path
a=2*pi*rand
r=sqrt(rand)
x=(rc*r)*cos(a)+x1
y=(rc*r)*sin(a)+y1
end
%to test the function
clf
axis equal
hold on
x1=1
y1=1
rc=1
[x,y,z] = cylinder(rc,200);
plot(x(1,:)+x1,y(1,:)+y1,'r')
for t=1:1000 %loop until doing 1000 points inside the circle
[x y]=cirrdnPJ(x1,y1,rc)
plot(x,y,'x')
%pause(0.01) %if you want to see the point being drawn
end
  7 个评论
显示 5更早的评论
Alan Weiss
Alan Weiss 2012-2-27
You might be interested that essentially this same code is in an obscure example in MATLAB documentation:
http://www.mathworks.com/help/toolbox/bioinfo/ug/bs3tbev-1.html#bs3tbev-15

请先登录，再进行评论。

更多回答（6 个）

Bruno Luong
Bruno Luong 2011-1-23
Careful, mathematically there is always a tiny probability the rejection method (Matt's proposal) returns number of points less than anticipated and it could lead to crash, or run forever if safety looping to fill would be implemented. In practice that might never happens, but if I was a designer of the code for flying an airplane or controlling a robot for brain surgery, then I would never use such code for a peace of mind (mine not the patient).
Here is a direct method
% Data
n = 10000;
radius = rand;
xc = randn;
yc = randn;
% Engine
theta = rand(1,n)*(2*pi);
r = sqrt(rand(1,n))*radius;
x = xc + r.*cos(theta);
y = yc + r.*sin(theta);
% Check
plot(x,y,'.')
% Bruno
Matt Fig
Matt Fig 2011-1-23
Here is a function to do this without skew:
function [X,Y] = rand_circ(N,x,y,r)
% Generates N random points in a circle.
% RAND_CIRC(N) generates N random points in the unit circle at (0,0).
% RAND_CIRC(N,x,y,r) generates N random points in a circle with radius r
% and center at (x,y).
if nargin<2
x = 0;
y = 0;
r = 1;
end
Ns = round(1.28*N + 2.5*sqrt(N) + 100); % 4/pi = 1.2732
X = rand(Ns,1)*(2*r) - r;
Y = rand(Ns,1)*(2*r) - r;
I = find(sqrt(X.^2 + Y.^2)<=r);
X = X(I(1:N)) + x;
Y = Y(I(1:N)) + y;
Now from the command line:
[x,y] = rand_circ(10000,4,5,3);
plot(x,y,'*r')
axis equal
  5 个评论
显示 3更早的评论
Amit Goel
Amit Goel 2022-5-7
why Ns = round(1.28*N + 2.5*sqrt(N) + 100); % 4/pi = 1.2732??
please share the logic behid
Torsten
Torsten 2022-5-7
编辑：Torsten 2022-5-8
If you choose 4/pi * N points between -r and r for X and Y, then approximately N points are inside the circle X^2+Y^2 <= r^2 and N*(4/pi-1) points are outside since the area of circle and square are pi*r^2 resp. 4*r^2. The 2.5*sqrt(N) + 100 - term is not clear to me (maybe a confidence level for the uniform distribution).

请先登录，再进行评论。

Bruno Luong
Bruno Luong 2011-1-23
Here is another non-rejection method based on Gaussian RANDN. To goal is to avoid calling sine/cosine. It seems only slightly faster according my test (2011A PreRe)
% Data
n = 1e6;
radius = rand;
xc = randn;
yc = randn;
% Engine 1, rand/cos/sin
tic
theta = rand(1,n)*(2*pi);
r = sqrt(rand(1,n))*radius;
x = xc + r.*cos(theta);
y = yc + r.*sin(theta);
t1=toc;
% Engine 2, randn
tic
x = randn(1,n);
y = randn(1,n);
r2 = x.^2+y.^2;
r = sqrt(rand(1,n)./r2)*radius;
x = xc + r.*x;
y = yc + r.*y;
t2=toc;
fprintf('rand/sin/cos method: %g [s]\n', t1);
fprintf('randn method : %g [s]\n', t2);
% Check
plot(x,y,'.')
% Bruno
Bruno
  1 个评论
Jeff Mason
Jeff Mason 2011-5-17
How do we know that gaussian x & y give uniform coverage of the circle?

请先登录，再进行评论。

Antti
Antti 2011-5-18
Here is a multidimensional version of the Brunos code:
%%Modified randn method for hypercircle
% Data
n = 1e6;
n_dim = 3;
pc = randn(1, n_dim);
% Engine 3, multidimensional randn
tic
p = randn(n, n_dim);
r2 = sum(p.^2, 2);
r = sqrt(rand(n, 1) ./ r2)*radius;
p = repmat(pc, n, 1) + repmat(r, 1, n_dim) .* p;
t3 = toc;
fprintf('randn, 3-dim, method : %g [s]\n', t3);
% Check
plot3(p(:,1), p(:,2), p(:,3), '.')
% Antti
  2 个评论
John D'Errico
John D'Errico 2019-1-24
I just saw this answer. It is incorrect, because it creates r incorrectly.
r = sqrt(rand(n, 1) ./ r2)*radius;
Use of sqrt there is a bad idea, since it will not produce uniform sampling in other numbers of dimensions than 2.
Do NOT use the above code in higher dimensions, as it is NOT a uniform sampling.
Bruno Luong
Bruno Luong 2019-2-16
编辑：Bruno Luong 2019-2-16
To generate uniform distribution in spherical shell on higher dimensions
V(a,b) := { X in R^m : a <= |X| <= b }
the correct modification is :
m = 4;
a = 2;
b = 3;
n = 1000; % number of points
s = randn(m,n);
r = (rand(1,n)*(b^m-a^m)+a^m).^(1/m);
c = r./sqrt(sum(s.^2,1));
X = bsxfun(@times, s, c); % (m x n) or simply X = s .* c on MATLAB with auto expansion

请先登录，再进行评论。

James Tursa
James Tursa 2011-1-25
FYI, if you wish to extend methods to higher dimensions, the randn methods are the way to go. Rejection methods in higher dimensions are prone to very small probabilities of being inside the hyper-sphere and thus needing many random numbers generated to get the desired number of samples.
  3 个评论
显示 1更早的评论
Aaron Chacon
Aaron Chacon 2011-5-16
Another problem with rejection methods is that if you reject enough numbers in the pseudo-random number sequence, the next numbers in the sequence can become predictable, and you may not fill the hyper-sphere uniformly.
Walter Roberson
Walter Roberson 2011-5-16
Is this a serious problem with anything more sophisticated than a linear congruential generator? In particular, with Twister being good up to some 600+ dimensions, is it worth thinking about for Twister?

请先登录，再进行评论。

Greg
Greg 2011-5-16
Hi guys, I had a bit of a simple question on the above code:
- has anyone tested whether using Sobol/Halton sequences to generate the random numbers improves (if at all possible) the above results?
I've been doing a bit of MC simulations for a project, and some literature suggests that these sequences are better than the 'rand' function (not sure how to measure whether this is the case or not, but graphic display is the only thing that springs to mind). Also unsure what the benefit of improved results v.s. processing time is.
Greg

类别

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

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by