114 views (last 30 days)

H. Paul Keeler
on 1 Oct 2018

Edited: H. Paul Keeler
on 26 Aug 2019

If the simulation windown/region has 1000 (or non random n) points, then it is not a Poisson process. It is a binomial point process, which happens to be also a homogeneous Poisson point process conditioned on having n points.

You treat this problem in polar coordinates. For each point, you uniformly choose a random angular coordinate (or component) on the interval . For the radial coordinate (or component), you first choose a random number on the unit interval , and then you --- and this is very important -- square root that random number and multiply it by the radius of the disk. You have to take the square root to assure uniform points (because the area of a disk/sector is proportional to the radius squared, not the radius).

I recently answered this question in detail on my blog, where I discuss how to simulate a (homogeneous) Poisson point process on a disk. The code is here:

%Simulation window parameters

r=1; %radius of disk

xx0=0; yy0=0; %centre of disk

areaTotal=pi*r^2; %area of disk

%Point process parameters

lambda=1000; %intensity (ie mean density) of the Poisson process

%Simulate Poisson point process

numbPoints=poissrnd(areaTotal*lambda);%Poisson number of points

theta=2*pi*(rand(numbPoints,1)); %angular coordinates

rho=r*sqrt(rand(numbPoints,1)); %radial coordinates

%Convert from polar to Cartesian coordinates

[xx,yy]=pol2cart(theta,rho); %x/y coordinates of Poisson points

%Shift centre of disk to (xx0,yy0)

xx=xx+xx0;

yy=yy+yy0;

%Plotting

scatter(xx,yy);

xlabel('x');ylabel('y');

axis square;

Hai Ngo Thanh
on 4 Sep 2019

I think that rho ~ U(0,1) so

rho=r*sqrt(rand(numbPoints,1)); instead of rho=r*sqrt(rand(numbPoints,1));

I don't still understand x,y ~ U(0,1) or rho,theta ~ Uniform distribution?

I would appreciate if you give me an explanation.

H. Paul Keeler
on 5 Sep 2019

I think you meant to write something else, as it looks like you wrote the same command twice.

Of course, this all depends if you want to simulate on a rectangle (use Cartesian coordinates) or disk (use polar coordinates).

Yes, in Cartesian both random variables are uniform, but they are NOT both uniform in polar coordinates.

In polar coordinates, rho is not uniformly distributed. Only theta is uniform. A reasonable guess would be to do the same as before and generate uniform random variables between zero and one, and then multiply them by the disk radius r. But that would be wrong.

We must take the square root because the area element of a sector or disk is proportional to the radius squared, and not the radius (recall that the area element is dA=r dr dtheta). These random numbers do not have a uniform distribution, due to the square root, but in fact their distribution is an example of the triangular distribution, which is defined with three real-valued parameters a, b and c, and for our case, set a=0 and b=c=r.

The third edition of the classic book Stochastic Geometry and its Applications by Chiu, Stoyan, Kendall and Mecke details on page 54 how to uniformly place points on a disk, which they call the radial way. The same simulation section appears in the previous edition by Chiu, Kendall and Mecke (Chiu didn’t appear as an author until the current edition), though these books in general have little material on simulation methods.

As I mentioned before, I answered this question in detail on my blog, where I discuss how to simulate a (homogeneous) Poisson point process on a disk.

Noman Aftab
on 20 Jan 2020

r=1; %radius of disk

xx0=0; yy0=0; %centre of disk

areaTotal=pi*r^2; %area of disk

%Point process parameters

lambda=10/areaTotal; %intensity (ie mean density) of the Poisson process

%Simulate Poisson point process

numbPoints=poissrnd(areaTotal*lambda)%Poisson number of points

theta=2*pi*(rand(numbPoints,1)); %angular coordinates

rho=r*sqrt(rand(numbPoints,1)) %radial coordinates

%Convert from polar to Cartesian coordinates

[xx,yy]=pol2cart(theta,rho); %x/y coordinates of Poisson points

%Shift centre of disk to (xx0,yy0)

xx=xx+xx0;

yy=yy+yy0;

%Plotting

scatter(xx,yy);

xlabel('x');ylabel('y');

axis square;

we have to divide lambda by area of disc.

Sign in to comment.

Astha
on 8 Jun 2019

H. Paul Keeler
on 27 Aug 2019

The first part is not too difficult. For example, for a homogeneous Poisson point on a rectangle:

%Simulation window parameters

xMin=0;xMax=1;

yMin=0;yMax=1;

xDelta=xMax-xMin;yDelta=yMax-yMin; %rectangle dimensions

areaTotal=xDelta*yDelta; %area of simulation window

%Point process parameters

lambda=10; %intensity (ie mean density) of the Poisson process

%Simulate Poisson point process

numbPoints=poissrnd(areaTotal*lambda);%Poisson number of points

xx=xDelta*(rand(numbPoints,1))+xMin;%x coordinates of Poisson points

yy=xDelta*(rand(numbPoints,1))+yMin;%y coordinates of Poisson points

xxyy=[xx(:) yy(:)]; %combine x and y coordinates

%Perform Voronoi tesseslation using built-in function

[vertexAll,cellAll]=voronoin(xxyy);

Variables: vertexAll is an array with the Cartesian ccordinates of all the vertices of the Voronoi tessellation. cellAll structure array, where each entry is an array of indices of the vertices describing the Voronoi tesselation; see MATLAB function voronoin.

If you want to plot a Voronoi diagram, you can just run the command:

voronoi(xx,yy); %create voronoi diagram on the PPP

For the second part, I would think a Poisson point process over the entire simulation window/region corresponds to an (invididual) independent Poisson point process in each Voronoi window. This should be true because of a defining property of the Poisson point process, which is sometimes called the scattering property (ie the number of points of a Poisson point process in non-overlapping/disjoint regions are independent variables.) This will NOT be true if it's not a Poisson point process.

In short, just simulate another Poisson point process over the entire region for the users (it can have a different intensity).

It becomes more tricky if you want to place a single point in each cell -- that will NOT be a Poisson point process over the entire simulation window. It's still possible to do this. Here are two methods.

1) The simplest and crudest method is to bound a Voronoi cell with a rectangle or disk, then place/sample the point uniformly on the rectangle or disk, as is done when sampling Poisson points on these geometric objects. If the single point doesn't randomly land inside the rectangle or disk, then start again. Crude, slightly inefficient, but simple.

2) The more elegant way is to first partition (ie divide) the Voronoi cells into triangles. (Each side of a Voronoi cell corresponds to one side of a triangle (ie two points), where the third point is the original Poisson point corresponding to the Voronoi cell.) Given the n triangles that form the Voronoi cell, randomly choose a triangle based on the areas of the triangle (ie the probabilities are the areas renormalized by the total area of the cell). In MATLAB, you would do this:

cdfTri=cumsum(areaTri)/areaCell; %create triangle CDF

indexTri=find(rand(1)<=cdfTri,1); %use CDF to randomly choose a triangle

Then for that single triangle, uniformly place a point on it. Repeat for all (bounded) Voronoi cells.

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 1 Comment

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/300022-i-want-to-spatially-distribute-1000-mobile-devices-in-a-network-according-to-poisson-point-process#comment_421922

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/300022-i-want-to-spatially-distribute-1000-mobile-devices-in-a-network-according-to-poisson-point-process#comment_421922

Sign in to comment.