Simulating a homogeneous Poisson point process on a rectangle

This is the first of a series of posts about simulating Poisson point processes. I’ll start with arguably the simplest Poisson point process on two-dimensional space, which is the homogeneous one defined on a rectangle. Let’s say that we we want to simulate a Poisson point process with intensity \(\lambda>0\) on a (bounded) rectangular region, for example, the rectangle \([0,w]\times[0,h]\) with dimensions \(w>0\) and \(h>0\) and area \(A=wh\). We assume for now that the bottom left corner of the rectangle is at the origin.

Steps

Number of points

The number of points in the rectangle  \([0,w]\times[0,h]\) is a Poisson random variable with mean \(\lambda A\). In other words, this random variable is distributed according to the Poisson distribution with parameter \(\lambda A\), and not just \(\lambda\), because the number of points depends on the size of the simulation region.

This is the most complicated part of the simulation procedure. 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. There’s a couple of different ways used to simulate Poisson random variables, but we will skip the details. In MATLAB, it is done by using 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.

Location of points

The points now need to be positioned randomly, which is done by using Cartesian coordinates. For a homogeneous Poisson point process, the \(x\) and \(y\) coordinates of each point are independent uniform points, which is also the case for the binomial point process, covered in a previous post. For the rectangle \([0,w]\times[0,h]\), the \(x\) coordinates are uniformly sampled on the interval \([0,w]\), and similarly for the \(y\) coordinates. If the bottom left corner of rectangle is located at the point \((x_0,y_0)\), then we just have to shift the random \(x\) and \(y\) coordinates by respectively adding \(x_0\) and \(y_0\).

Every scientific programming language has a random uniform number generator because it is the default random number generator. In MATLAB, R and SciPy, it is respectively rand, runif and scipy.stats.uniform.

Code

Here is some code that I wrote for simulating a homogeneous Poisson point process on a rectangle. You will notice that in all the code samples the part that simulates the Poisson point process requires only three lines of code: one line for the number of points and two lines lines for the \(x\) and \(y\) coordinates  of the points.

MATLAB
%Simulation window parameters
xMin=0;xMax=1;
yMin=0;yMax=1;
xDelta=xMax-xMin;yDelta=yMax-yMin; %rectangle dimensions
areaTotal=xDelta*yDelta;

%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
xx=xDelta*(rand(numbPoints,1))+xMin;%x coordinates of Poisson points
yy=yDelta*(rand(numbPoints,1))+yMin;%y coordinates of Poisson points

%Plotting
scatter(xx,yy);
xlabel('x');ylabel('y');
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
xMin=0;xMax=1;
yMin=0;yMax=1;
xDelta=xMax-xMin;yDelta=yMax-yMin; #rectangle dimensions
areaTotal=xDelta*yDelta;

#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
xx=xDelta*runif(numbPoints)+xMin;#x coordinates of Poisson points
yy=yDelta*runif(numbPoints)+yMin;#y coordinates of Poisson points

#Plotting
plot(xx,yy,'p',xlab='x',ylab='y',col='blue');
Python

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

import numpy as np
import scipy.stats
import matplotlib.pyplot as plt

#Simulation window parameters
xMin=0;xMax=1;
yMin=0;yMax=1;
xDelta=xMax-xMin;yDelta=yMax-yMin; #rectangle dimensions
areaTotal=xDelta*yDelta;

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

#Simulate Poisson point process
numbPoints = scipy.stats.poisson( lambda0*areaTotal ).rvs()#Poisson number of points
xx = xDelta*scipy.stats.uniform.rvs(0,1,((numbPoints,1)))+xMin#x coordinates of Poisson points
yy = yDelta*scipy.stats.uniform.rvs(0,1,((numbPoints,1)))+yMin#y coordinates of Poisson points
#Plotting
plt.scatter(xx,yy, edgecolor='b', facecolor='none', alpha=0.5 )
plt.xlabel("x"); plt.ylabel("y")
Julia

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

Results

MATLAB

PoissonMatlab

R

PoissonR

Python

PoissonPython

Higher dimensions

If you want to simulate a Poisson point process in a three-dimensional box (typically called a cuboid or rectangular prism), you just need two modifications.

For a box \([0,w]\times[0,h]\times[0,\ell]\), the number of points now a Poisson random variable with mean \(\lambda V\), where \(V= wh\ell\) is the volume of the box. (For higher dimensions, you need to use $n$-dimensional volume.)

To position the points in the box, you just need an additional uniform variable for the extra coordinate. In other words, the \(x\), \(y\) and \(z\) coordinates are uniformly and independently sampled on the respective intervals \([0,w]\), \([0,h]\), \([0,\ell]\). The more dimensions, the more uniform random variables.

And that’s it.

Further reading

For simulation of point processes, see, for example, the books Statistical Inference and Simulation for Spatial Point Processes by Møller and Waagepetersen, or Stochastic Geometry and its Applications by Chiu, Stoyan, Kendall and Mecke.

There are books written by spatial statistics experts such as Stochastic Simulation by Ripley and Spatial Point Patterns: Methodology and Applications with R by Baddeley, Rubak and Turner, where the second book covers the spatial statistics R-package spatstat. Kroese and Botev also have a good introduction in the edited collection Stochastic Geometry, Spatial Statistics and Random Fields : Models and Algorithms by Schmidt, where the relevant chapter (number 12) is also freely available online.

More general stochastic simulation books that cover relevant material include Uniform Random Variate Generation by Devroye and Stochastic Simulation: Algorithms and Analysis by Asmussen and Glynn.

Point process simulation

Point processes are mathematical objects that can represent a collection of randomly located points on some underlying space. There is a rich range of these random objects, and, depending on the type of points process, there are various steps and methods used to sample, simulate or generate them on computers. (I use these three terms somewhat interchangeably.) This is the first of a series of posts in which I will describe how to simulate some of the more tractable point processes defined on bounded regions of two-dimensional Euclidean space.

Finite point processes

A general family of point processes are the finite point processes, which are, as the name suggests, simply point processes with the property that the total number of points is finite with probability one. Usually it is assumed that both the number of points and the positions of these points are random. These two properties give a natural way to describe and to study point processes mathematically, while also giving an intuitive way to simulate them on computers.

If the number and locations of points can be simulated sequentially, such that the number of points is naturally generated first followed by the placing of the points, then the simulation process is usually easier. Such a point process is a good place to start for an example.

Simple example — Binomial point process

The binomial point process is arguably the simplest point process. It is created by scattering \(n\) points uniformly and independently located in some bounded region, say, \(R\) with area \(|R|\). The number of points located in a region \(B\subset R\) is a binomial random variable, say, \(N\), where its probability parameter \(p=\frac{|B|}{|R|}\) is simply the ratio of the  two areas. This implies that

$$ P(N=k)={n\choose k} p^k(1-p)^{n-k}, $$

The total number of points is non-random number \(n\), so we do not need to generate it as it is given. The independence of point locations means that all the points can be positioned in parallel.

For example, to simulate \(n\) points of a binomial point process on the unit square \([0,1]\times[0,1]\), we just need to independently sample the \(x\) and \(y\) coordinates of the points from a uniform distribution on the unit interval \([0,1]\). In other words, we randomly generate or sample two sets of \(n\) uniform random variables corresponding to the \(x\) and \(y\) of the \(n\) points.

To sample on a general \(w \times h\) rectangle, we just need need to multiple the random \(x\) and \(y\) values by the respective dimensions of the rectangle \(w\) and \(h\). Of course the rectangle can also be shifted up or down by adding or subtracting \(x\) and \(y\) values appropriately.

Essentially every programming language has a function for generating uniform random numbers because it is the default random number generator, which all the other random number generators build off, by employing various techniques like applying transformations. In MATLAB, R and Python, it is respectively rand, runif and scipy.stats.uniform, which all generate uniform points on the open interval \((0,1)\).

Code

Here are some pieces of code that illustrates how to simulate a binomial point process on the unit square. I suggest downloading the code directly here from my repository. (Sometimes my webpage editor mangles simulation code.)

MATLAB

n=10; %number of points

%Simulate binomial point process
xx=rand(n,1);%x coordinates of binomial points
yy=rand(n,1);%y coordinates of binomial points

%Plotting
scatter(xx,yy);
xlabel('x');ylabel('y');

R

n=10; #number of points
#Simulate Binomial point process
xx=runif(n);#x coordinates of binomial points
yy=runif(n);#y coordinates of binomial points

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

Python

&lt;/pre&gt;
import numpy as np
import matplotlib.pyplot as plt

plt.close('all'); # close all figures

numbPoints = 10; # number of points

# Simulate binomial point process
xx = np.random.uniform(0, 1, numbPoints) # x coordinates of binomial points
yy = np.random.uniform(0, 1, numbPoints) # y coordinates of binomial points

# Plotting
plt.scatter(xx, yy, edgecolor='b', facecolor='none', alpha=0.5)
plt.xlabel('x')
plt.ylabel('y')
plt.show()
&lt;pre&gt;

More complicated point processes

The binomial point process is very simple. One way to have a more complicated point process is to have a random number of points. The binomial point process is obtained from conditioning on the number of points of a homogeneous Poisson point process, implying its simulation builds directly off the binomial point process. The homogeneous Poisson point process is arguably the most used point process, and it is then in turned used to construct even more complicated points processes. This suggests that the homogeneous Poisson point process is the next point process we should try to simulate.

Further reading

For simulation of point processes, see, for example, the books Statistical Inference and Simulation for Spatial Point Processes by Møller and Waagepetersen, or Stochastic Geometry and its Applications by Chiu, Stoyan, Kendall and Mecke. There are books written by spatial statistics experts such as Stochastic Simulation by Ripley and Spatial Point Patterns: Methodology and Applications with R by Baddeley, Rubak and Turner, where the second book covers the spatial statistics R-package spatstat. Kroese and Botev also have a good introduction in the edited collection Stochastic Geometry, Spatial Statistics and Random Fields : Models and Algorithms by Schmidt, where the relevant chapter (number 12) is also freely available online. More general stochastic simulation books that cover relevant material include Uniform Random Variate Generation by Devroye and Stochastic Simulation: Algorithms and Analysis by Asmussen and Glynn.

For a mathematical introduction to finite point processes, see the standard reference An Introduction to the Theory of Point Processes by Daley and Vere-Jones (Chapter 5 in either the single-volume edition or the first volume of the second two-volume edition). Finite point processes are also covered in a recent tutorial on Palm calculus and Gibbs point processes.