I want to spatially distribute 1000 mobile devices in a network according to Poisson Point Process.

I have 1000 mobile devices or Users. I want to spatially distribute them in a Network(Area detail given in next line) according to Poisson Point Process where the density rate or Lambda is 100. Network is a circular Area of Radius 1000m. poissrnd function only gives me a random variable. What can I do with it to spatially distribute them?. Please inform regarding it. I would be obliged.

[KalMandy]

Hi, have you solved this?

Answer 1

[H. Paul Keeler]

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]

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]

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]

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.

[Bidyarani]

r = 1; is it correct for radius 1000 meter? or it should be r = 1 only?

Answer 2

[Astha]

Can anyone tell how to generate voronoi diagram with centers generated according to poisson point process and then users are substituted in each voronoi cell again according to poisson point process. If matlab code could be provided then it would be of great help. Thanks in advance.

[H. Paul Keeler]

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.

[Astha]

Thank you so much. Just trying it. Will get back once understood.