Simulating a Poisson point process on a disk

Sometimes one needs to simulate a Poisson point process on a disk. I know I often do. A disk or disc, depending on your spelling preference, is isotropic or rotationally-invariant, so a lot of my simulations of a Poisson point process happen in a circular simulation window when I am considering such a setting. For example, maybe you want to consider a single wireless receiver in a Poisson network of wireless transmitters, which only cares about the distance to a transmitter. Alternatively, maybe you want to randomly sprinkle a virtual cake. What to do? A Poisson point process on a disk is the answer.

In this post I will simulate a Poisson point process with intensity $$\lambda>0$$ on a disk with radius $$r>0$$.

Steps

The simulation steps are very similar to those in the previous post where I simulated a  homogeneous Poisson point process on a rectangle, and I suggest going back to that post if you are not familiar with the material. The main difference between simulation on a rectangle and a disk is the positioning of the points, but first we’ll look at the number of points.

Number of points

The number of points of a Poisson point process falling within a circle of radius $$r>0$$ is a Poisson random variable with mean  $$\lambda A$$, where $$A=\pi r^2$$ is the area of the disk. As in the rectangular case, this is the most complicated part of the simulation procedure. But as long as your preferred programming language can produce (pseudo-)random numbers according to a Poisson distribution, you can simulate a homogeneous Poisson point process on a disk.

To do this in MATLAB,  use the poissrnd function with the argument $$\lambda A$$. In R, it is done similarly with the standard  function rpois . In Python, we can use either the scipy.stats.poisson or numpy.random.poisson function from the SciPy or NumPy libraries.

Locations of points

We need to position all the points randomly and uniformly on a disk. In the case of the rectangle, we worked in Cartesian coordinates. It is then natural that we now work in polar coordinates.  I’ll denote the angular and radial coordinate respectively by $$\theta$$ and $$\rho$$. To generate the random angular (or $$\theta$$) values, we simply produce uniform random variables between zero and one, which is what all standard (pseudo-)random number generators do in programming languages. But we then multiply all these numbers by $$2\pi$$, meaning that all the numbers now fall between $$0$$ and $$2\pi$$.

To generate the random radial (or $$\rho$$) values, 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.

Before multiplying uniform random variables by the radius, we must take the square root of all the random numbers. We then multiply them by the radius, generating random variables between $$0$$ and $$r$$. (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.) 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$$.

In summary, if $$U$$ and $$V$$ are two independent uniform random variables on $$(0,1)$$, then random point located uniformly on a disk of radius $$r$$ has the polar coordinates $$(r\sqrt(U), 2\pi V)$$.

From polar to Cartesian coordinates

That’s it. We have generated polar coordinates for points randomly and uniformly located in a disk. But to plot the points, often we have to convert coordinates back to Cartesian form. This is easily done in MATLAB with the pol2cart function. In languages without such a function, trigonometry comes to the rescue: $$x=\rho\cos(\theta)$$ and $$y=\rho\sin(\theta)$$.

Equal x and y axes

Sometimes the plotted points more resemble points on an ellipse than a disk due to the different scaling of the x and y axes. To fix this in MATLAB, run the command: axis square. In Python, set axis(‘equal’) in your plot; see this page for a demonstration.

Code

I’ll now give some code in MATLAB, R and Python, which, as you can see, are all very similar

MATLAB
%Simulation window parameters
xx0=0; yy0=0; %centre of disk

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

%Point process parameters
lambda=100; %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

%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;
R

Note: it is a bit tricky to write “<-” in the R code (as it automatically changes to the html equivalent in the HTML editor I am using), so I have usually used “=” instead of the usual “<-”.

#Simulation window parameters
xx0=0; yy0=0; #centre of disk

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

#Point process parameters
lambda=100; #intensity (ie mean density) of the Poisson process

#Simulate Poisson point process
numbPoints=rpois(1,areaTotal*lambda);#Poisson number of points
theta=2*pi*runif(numbPoints);#angular  of Poisson points

#Convert from polar to Cartesian coordinates
xx=rho*cos(theta);
yy=rho*sin(theta);

#Shift centre of disk to (xx0,yy0)
xx=xx+xx0;
yy=yy+yy0;

#Plotting
par(pty="s")
plot(xx,yy,'p',xlab='x',ylab='y',col='blue');

Of course, with the amazing spatial statistics library spatstat, simulating a Poisson point process in R is even easier.

#create Poisson "point pattern" object
plot(ppPoisson); #Plot point pattern object
#retrieve x/y values from point pattern object
xx=ppPoisson$x; yy=ppPoisson$y;

Actually, you can even do it all in two lines: one for loading the spatstat library and one for creating and plotting the point pattern object.

Python

Note: “lambda” is a reserved word in Python (and other languages), so I have used “lambda0” instead.

import numpy as np; #NumPy package for arrays, random number generation, etc
import matplotlib.pyplot as plt #for plotting

#Simulation window parameters
xx0=0; yy0=0; #centre of disk
areaTotal=np.pi*r**2; #area of disk

#Point process parameters
lambda0=100; #intensity (ie mean density) of the Poisson process

#Simulate Poisson point process
numbPoints = np.random.poisson(lambda0*areaTotal);#Poisson number of points
theta=2*np.pi*np.random.uniform(0,1,numbPoints); #angular coordinates

#Convert from polar to Cartesian coordinates
xx = rho * np.cos(theta);
yy = rho * np.sin(theta);

#Shift centre of disk to (xx0,yy0)
xx=xx+xx0; yy=yy+yy0;

#Plotting
plt.scatter(xx,yy, edgecolor='b', facecolor='none', alpha=0.5 );
plt.xlabel("x"); plt.ylabel("y");
plt.axis('equal');
Julia

After writing this post, I later wrote the code in Julia. The code is here and my thoughts about Julia are here.

Results

Python

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 Stoyan, 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. There is the book Spatial Point Patterns: Methodology and Applications with R written by spatial statistics experts  Baddeley, Rubak and Turner, which covers the spatial statistics (and point process simulation) R-package spatstat.

(Visited 3,717 times, 1 visits today)

21 thoughts on “Simulating a Poisson point process on a disk”

1. Ash says:

Can you please tell how to generate voronoi diagram in Matlab 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.

1. MATLAB already has an inbuilt function for this:

https://au.mathworks.com/help/matlab/ref/voronoi.html

I recommend you look at their examples first. I have code that implements MATLAB’s voronoi function and does something similar to you what you want, but it needs tidying up. Perhaps I’ll post it sometime in the near future.

For the second part, I believe a Poisson point process over the entire simulation window/region corresponds to an (individual) 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 intensity different from the first Poisson point processed used to create the Voronoi tessellation/diagram).

2. yassine says:

Hi Paul,
Please consider the event “only the nearest point to the center of the disc, is displaced in a random direction, independently from the other points”. Can you please tell if the new point process remains a PPP.

1. Paul Keeler says:

Hi Yassine,

Doubtful. You are changing the homogeneous PPP in the region around the center of the disk, where there is now *less* or *more* chance (depending on the random displacement) of finding a point. Although perhaps the change is very small.

The only way I see it still being a PPP is by positioning it with the same random distribution that it was originally positioned with: (in polar coordinates) the angular variable $\Theta$ is uniform on $(0,2\pi)$ and the radial component $R$ has a Rayleigh distribution $P(R\leq r) = 1- exp(- \lambda\pi r^2 )$ (which is the distribution of the nearest-neighbour distance ). But then you haven’t changed the Poisson point process (in distribution) as it statistically behaves the same.

1. Paul Keeler says:

This code supports my claim:

numbSim=10^4; %number of simulations

%Simulation window parameters
xx0=0; yy0=0; %centre of disk
areaTotal=pi*r^2; %area of disk

%Point process parameters
lambda=5; %intensity (ie mean density) of the Poisson process

countPointsMid=zeros(numbSim,1);
rMid=r/2; %arbitrary point along the radius
%loop through for many simulations
for ss=1:numbSim
%Simulate Poisson point process
numbPoints=poissrnd(areaTotal*lambda);%Poisson number of points
theta=2*pi*(rand(numbPoints,1)); %angular coordinates

[rhoMin, indexMin]=min(rho);

%change rho and theta -- comment next 4 commands for regular Poisson
%generate
thetaNew=2*pi*rand(1);
rhoNew=r*rand(1); %uniform rho
%rhoNew=sqrt(-log(rand(1))/lambda/pi); %Rayleigh rho
%replace new variables
theta(indexMin)=thetaNew;
rho(indexMin)=rhoNew;

countPointsMid(ss)=sum(rho

3. Hawat says:

how can we simulate a poisson point rocess in R^3

1. Paul Keeler says:

Hi Hawat,

You would need to do it in a cube or rectangular prism. It’s just a variation of doing it on a rectangle, but you use the volume (instead of the area) of the cuboid (or rectangular prism). And then for each point you generate a third uniform random variable Z for the three Cartesian coordinate. Look at this post, where the Poisson point process is simulated on a rectangle:

https://hpaulkeeler.com/poisson-point-process-simulation/

I have included a section at the end detailing how to do it in a bounded three-dimensional box (that is, cuboid or rectangular prism).

Alternatively, you could also do it in a 3-D sphere. This post does it for a sphere (the *surface* of a ball), but I detail at the end of the post how you can do position points uniformly *inside* a ball:

https://hpaulkeeler.com/simulating-a-poisson-point-process-on-a-sphere/

4. Mots'oane Nkhabu says:

Hi Paul!

I have a two-tier network comprising devices, gateways and the base station. Gateways are distributed according to PPP while the devices themselves are distributed according MCP. How do I achieve this in MATLAB assuming a circular area of coverage for each gateway generated, within which the respective devices are distributed wrt MCP?

1. Shuchi Tripathi says:

Sir, thank you for sharing this article. I have used the MATLAB code to generate a base station PPP inside B(0,R) where each BS has a mobile user at a fixed distance ‘r’ from its associated BS. To define mobile users, I have used

rhoM = rho + r;
thetaM = theta + unifrnd(0, 2*pi);

After that I have used pol2cart operator to get the coordinates of the mobile user but I am not getting the required result. May you suggest anything regarding this problem?

I am also thinking to generate a high density mobile user point process separatly and then used thing process to get the target mobile user the given fixed distance. Can I simulate this problem in this way?

Thanks and Regards

1. Paul Keeler says:

I don’t understand what you’re doing. In polar coordinates, you cannot simply add the radial (rho) and angular (theta) components. You do that in Cartesian coordinates. If you look at my MATLAB code, that’s what I do:

%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;

https://github.com/hpaulkeeler/posts/blob/master/PoissonDisk/PoissonDisk.m

5. Shuchi Tripathi says:

Sir, thank you for sharing this article. I have used the MATLAB code to generate a base station PPP inside B(0,R) where each BS has a mobile user at a fixed distance ‘r’ from its associated BS. To define mobile users, I have used

rhoM = rho + r;
thetaM = theta + unifrnd(0, 2*pi);

After that I have used pol2cart operator to get the coordinates of the mobile user but I am not getting the required result. May you suggest anything regarding this problem?

I am also thinking to generate a high density mobile user point process separatly and then used thining to get the target mobile user the at the fixed distance ‘r’. Can I simulate this problem in this way?

Thanks and Regards

6. Omar THABET says:

Hello Prof. many thanks for sharing these codes
I assumed I have BS at the center of the the circle and users should be deployed at least 35m from the BS. as well as I assumed this circle is compromised of two circles, which means there are inner cercle used for deplying BS-users and out circle used to deply another users. So, can I impement something like that assumtion ?

7. Wray says:

I also have a question about the simulation. In PPP, the number of points is a variable, maybe we should not use the mean of the variable to replace the points number. In each cycle of executing the codes, the number of points should be different.

1. Paul Keeler says:

Thank you for your comment. In the code, the number of points is a random variable. It changes every simulation. It is not a constant.

But in the code you pass a constant (the mean) to a function that generates Poisson random variables (or more correctly, variates). In other words, the input is constant, but the output is random. To do this in MATLAB, use the poissrnd function with the argument λA. In R, it is done similarly with the standard function rpois. In Python, you can use either the scipy.stats.poisson or numpy.random.poisson function from the SciPy or NumPy libraries.

1. Wray says:

Thanks for your detailed interpretation with patience, Sir. It’s really helpful to my work.

1. Paul Keeler says:

No problem.

8. Amira says:

Hi ,Sir, thanks for publishing this helpful article, but I have a problem, I made a simulation for uplink system in square and I want to make the same on disk but the results are not the same, is the shape changes the result.

1. Paul Keeler says:

The results will be different due to “edge effects” (you are approximating the infinite plane with a square or disk). But the results should be the same as the disk and square grow in size.

Hi Paul,

I have a two tiers wireless network , drones and ground base stations, both are distributed according to BPP in a circular area with radius R. With a number of users distributed randomly on the ground. How to simulate in MATLAB please? Your help is really appreciated.

1. Paul Keeler says:

Thank you for your comment. For the drones (or base stations), you just need to simulate the binomial point process on a disk. You take the the simulation method for a Poisson point process on a disk and replace the random (Poisson) variable of points, say $N$, with a constant number of points, say $n$.

See the posts:

https://hpaulkeeler.com/simulating-a-poisson-point-process-on-a-disk/

https://hpaulkeeler.com/binomial-point-process/

Here is the corresponding code for the two above posts:

https://github.com/hpaulkeeler/posts/tree/master/PoissonDisk

https://github.com/hpaulkeeler/posts/tree/master/BinomialRectangle