random points within circle with fixed density per m^2 (matlab)

Question

I have to generate in matlab, random data points uniformly distributed within a circle centered at [0,0] and radius r=30m I want these points to have density 0.3 points/m^2.

The code I wrote is this:

n = 848;
Rc = 30;
Xc = 0;
Yc = 0;
theta = rand(1,n)*(2*pi);
r = Rc*sqrt(rand(1,n));  
x = Xc + r.*cos(theta);
y = Yc + r.*sin(theta);
plot(x,y,'.'); axis square;

I used a fixed number of datapoints (848 points), to ensure that the density INSIDE the circle area remains 0.3.

However, after plotting the resulting coordinates I see that they are not equally distributed per m^2 and I wonder what am I doing wrong??

The problem is that I can not fix the density per m^2... can anybody help please??

Could you clarify a bit more what kind of density you want the points to have? Just the fixed number of points inside the circle, or something else? If the first, by visual inspection of the plot, the result of your code seems fine to me. — mikkola, Jan 03 '16 at 21:57
thanks for the reply. I want every posible squared meter of the circle area to have density 0.3 and not only the whole circle (which in my case is ok). this means that if I take one random squared meter I want to find 0.3 points inside (not possible in one :P but just generalize this to 10m^2 or more) Is that the clarification you wanted?? — K.S., Jan 03 '16 at 23:00
You can (i) generate a 2-D grid of points of a particular density, all the points being equidistant, (ii) add to the (x,y) coordinates of each point a particular amount of randomness, and then (iii) exclude the points outside the circle. This will ensure that, on average, the density per meter square is what you have defined initially. — marsei, Jan 03 '16 at 23:45

brainkz · Answer 1 · 2016-01-03T23:55:49.723

Imagine an extensible ring made of rubber. Uniformly distribute several points on this ring. When it has small radius, the density of points is high. But when you try to increase the radius or ring by stretching it the density of points on the ring would fall. That is the main problem in your algorithm.

Instead, try the following approach: 1. Draw a biggest possible circle on a square piece of paper 2. Put a point randomly on the paper 3. If it's outside a circle reject point and go to step 2 4. Repeat steps 2-3 until you get the necessary density

The above algorithm is implemented using complex numbers with loops:

n = 848;
Rc = 30;
Xc = 0;
Yc = 0;
x = 2*(rand(1,n) + 1i*rand(1,n)) - 1 - 1i;
for i = 1:length(x)
    while abs(x(i)-Xc-1i*Yc) > 1
        x(i) = 2*(rand + 1i*rand) - 1 - 1i;
    end
end
x = Rc * x + Xc + 1i*Yc;
plot(x,'.');
t = linspace(0,2*pi);
hold on
plot(Rc*cos(t)+Xc,Rc*sin(t)+Yc)
axis equal;
hold off;

The 5th line generates 848 random points with x and y coordinates between -1 and 1. Then, the loop checks each random point for being inside unit circle and if not, replaces it with another random point. Finally, we multiply your random points by 30 and shift them along X and Y axes (zero shift in your case). This you get as a result:

UPD to OP's comment: it's not exactly 0.3 per m^2 for any m^2, but it's close to this number. If you want this to be checked using squares, then you would define it as some range along x and some range along y. Then you can use:

xmin = 5;
xmax = 20;
ymin = 5;
ymax = 20;

if all(abs([xmin + ymin*1i, xmin + ymax*1i, xmax + ymin*1i, xmax + ymax*1i]) < Rc)
    density = nnz(real(x) > xmin & real(x) < xmax & imag(x) > ymin & imag(x) < ymax)/(xmax-xmin)/(ymax-ymin)
else
    error('Square is outside circle')
end

density =

    0.2800

So, you manually define your x and y limits that form your test square. Then, the square are checked to be inside circle and if it is, the density is checked.

using `i` as a loop variable is [generally frowned upon in MATLAB](http://stackoverflow.com/questions/14790740/using-i-and-j-as-variables-in-matlab), even when not using complex numbers, but now it's just begging for confusing and errors. — Adriaan, Jan 03 '16 at 21:18
thank you for your time but Isn't that the same thing that I did before? How can you prove that the density is 0,3 per m^2 in every possible m^2 within the circle? — K.S., Jan 03 '16 at 23:04
@Adriaan thanks for remark, but you focus on the grammar, instead of focusing on content. I read the linked question. In my code I did what's recommended and what I usually do: used `1i` instead of `i`. This action alone is enough to prevent most associated errors. — brainkz, Jan 03 '16 at 23:07
I'm tempted to downvote your answer because of `i` (intentionally perpetuating a harmful practice), but I can't honestly do that since your answer is otherwise "useful". But please do consider changing from `i` to `k` in your loop. Accidentally writing `2*(rand + i*rand) - 1 - 1i` in there would give you a very different result, and it would be a pain to track down such a subtle bug. Of course, you can say "I never make any errors", but then you should just call every variable `a, b, c` etc. for brevity (just to use a bit of a [slippery slope](https://yourlogicalfallacyis.com/slippery-slope)) — Andras Deak -- Слава Україні, Jan 04 '16 at 11:56

KKS · Accepted Answer · 2016-01-04T05:40:38.853

In order to see density shape, you need sufficient samples.

One 848 sample set would not be uniform, but if you increase your samples, you can see uniformly distributed shape.

I will show you how to evaluate your code with below code.

    n = 848;
    Rc = 30;
    Xc = 0;
    Yc = 0;

    minX=-30;
    maxX=30;
    minY=-30;
    maxY=30;
    delta=5; %%% area 5m^2


    sum_PDF=0;

    iterN=100; % 100 times average to get density

    for a=1:iterN;


%%% generate uniform distribution in circle%%%

        theta = rand(1,n)*(2*pi);
        r = Rc*sqrt(rand(1,n));  
        x = Xc + r.*cos(theta);
        y = Yc + r.*sin(theta);

        d=[x;y];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%get density function per each iteration%%%

        axisX=minX:delta:maxX;
        axisY=minY:delta:maxY;

        nPDF_X=length(axisX);
        nPDF_Y=length(axisY);

        PDF=zeros(nPDF_X,nPDF_Y);

        temp=0;
        count_i=1;
        count_j=1;

        for i=axisX;

             lowlimit_x=i-delta/2;
             upperlimit_x =i+delta/2;

                 for j=axisY;

                    lowlimit_y=j-delta/2;
                    upperlimit_y =j+delta/2;
                    temp=0;

                    for k=1:length(d(1,:)); 

                        if lowlimit_x<=d(1,k) & d(1,k)<upperlimit_x
                             if lowlimit_y<=d(2,k) & d(2,k)<upperlimit_y
                                 temp=temp+1;
                             else 
                             end
                        else
                        end
                    end

                    PDF(count_i,count_j)=temp;
                    count_j=count_j+1;
                 end
            count_i=count_i+1;
            count_j=1;
        end

        randVar_X=minX:delta:maxX;
        randVar_Y=minY:delta:maxY;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%Summing all density functions%%%%%%%%%

        sum_PDF=sum_PDF+PDF;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    end

%%%%%%%%%%%%%%%number per 1m^2%%%%%%%%%%%%%%%%

    meanPDF=(1/delta^2)*(1/iterN)*sum_PDF; 

    surf(randVar_X,randVar_Y,meanPDF')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I calculate number in 5m^2 over all area, 100 times. and, I averaged all density function and estimated number in 1m^2. As a result, your code well generates uniform distribution in circle with 0.3 in 1m^2.

The final density function is as follows.

random points within circle with fixed density per m^2 (matlab)

2 Answers2