## Simulating a Matérn cluster point process

A Matérn cluster point process is a type of cluster point process, meaning that its randomly located points tend to form random clusters. (I skip the details here, but by using techniques from spatial statistics, it is possible to make the definition of clustering more precise.) This point process is an example of a family of cluster point processes known as Neyman-Scott point processes, which have been used in spatial statistics and telecommunications.

## Overview

Simulating a Matérn cluster point process requires first simulating  a homogeneous Poisson point process with an intensity $$\lambda>0$$ on some simulation window, such as a rectangle, which is the simulation window I will use here. Then for each point of this underlying point process, simulate a Poisson number of points with mean $$\mu>0$$ uniformly on a disk with a constant radius $$r>0$$. The underlying  point process is sometimes called the parent (point) process, and its points are centres of the cluster disks.

The subsequent point process on all the disks is called daughter (point) process and it forms the clusters. I have already written about simulating the homogeneous Poisson point processes on a rectangle and a disk, so  those posts are good starting points, and I will not focus  too much on details for these steps steps.

The main challenge behind sampling this point process, which I originally forgot about in an earlier version of this post, is that it’s possible for daughter points to appear in the simulation window that come from parents points outside the simulation window. In other words, parents points outside the simulation window contribute to points inside the window.

To remove these edge effects, the point processes must be simulated on an extended version of the simulation window. Then only the daughter points within the simulation window are kept and the rest are removed.  Consequently, the points are simulated on an extended window, but we only see the points inside the simulation window.

To create the extended simulation window, we can add a strip of width $$r$$ all around the simulation window. Why? Well, the distance $$r$$ is the maximum distance from the simulation window that a possibly contributing parent point (outside the simulation window) can exist, while still having daughter points inside the simulation window. This means it is impossible for a hypothetical parent point beyond this distance (outside the extended window) to generate a daughter point that can fall inside the simulation window.

## Steps

##### Number of points

Simulate the underlying or parent Poisson point process on the rectangle with $$N_P$$ points. Then for each point, simulate a Poisson number of points, where each disk now has $$D_i$$ number of points. Then the total number of points is simply $$N=D_1+\dots +D_{P}=\sum_{i=1}^{N_P}D_i$$. The random variables $$P$$ and $$D_i$$ are Poisson random variables with respective means $$\lambda A$$ and $$\mu$$, where $$A$$ is the area of the rectangular simulation window. To simulate these random variables in MATLAB, use the poissrnd function. To do this in R, use the  standard function rpois. In Python, we can use either functions scipy.stats.poisson or numpy.random.poisson from the SciPy or NumPy libraries.

##### Locations of points

The points of the parent point process are randomly positioned 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. The points of all the daughter point process are randomly positioned by using polar  coordinates. For a homogeneous Poisson point process, the $$\theta$$ and $$\rho$$ coordinates of each point are independent  variables,  respectively with uniform and triangle distributions, which was covered in a previous post. Then we  convert coordinates back to Cartesian form, which is easily done in MATLAB with the pol2cart function. In languages without such a function: $$x=\rho\cos(\theta)$$ and $$y=\rho\sin(\theta)$$.

##### Shifting all the points in each cluster disk

In practice (that is, in the code), all the daughter points are simulated in a disk with its centre at the origin. Then for each cluster disk, all the points have to be shifted to the origin is the center of the cluster, which completes the simulation step.

To use vectorization in the code,  the coordinates of each cluster point are repeated by the number of daughters in the corresponding cluster by using the functions repelem in MATLAB, rep in R, and repeat in Python.

## Code

I’ll now give some code in MATLAB, R and Python, which you can see are all very similar. It’s also located here.

##### MATLAB
%Simulation window parameters
xMin=-.5;xMax=.5;
yMin=-.5;yMax=.5;

%Parameters for the parent and daughter point processes
lambdaParent=10;%density of parent Poisson point process
lambdaDautgher=100;%mean number of points in each cluster

%Extended simulation windows parameters
xMinExt=xMin-rExt;
xMaxExt=xMax+rExt;
yMinExt=yMin-rExt;
yMaxExt=yMax+rExt;
%rectangle dimensions
xDeltaExt=xMaxExt-xMinExt;
yDeltaExt=yMaxExt-yMinExt;
areaTotalExt=xDeltaExt*yDeltaExt; %area of extended rectangle

%Simulate Poisson point process for the parents
numbPointsParent=poissrnd(areaTotalExt*lambdaParent,1,1);%Poisson number
%x and y coordinates of Poisson points for the parent
xxParent=xMinExt+xDeltaExt*rand(numbPointsParent,1);
yyParent=yMinExt+yDeltaExt*rand(numbPointsParent,1);

%Simulate Poisson point process for the daughters (ie final poiint process)
numbPoints=sum(numbPointsDaughter); %total number of points

%Generate the (relative) locations in polar coordinates by
%simulating independent variables.
theta=2*pi*rand(numbPoints,1); %angular coordinates

%Convert from polar to Cartesian coordinates
[xx0,yy0]=pol2cart(theta,rho);

%replicate parent points (ie centres of disks/clusters)
xx=repelem(xxParent,numbPointsDaughter);
yy=repelem(yyParent,numbPointsDaughter);
%translate points (ie parents points are the centres of cluster disks)
xx=xx(:)+xx0;
yy=yy(:)+yy0;

%thin points if outside the simulation window
booleInside = ((xx >= xMin) & (xx <= xMax) & (yy >= yMin) & (yy <= yMax));
%retain points inside simulation window
xx=xx(booleInside);
yy=yy(booleInside);

%Plotting
scatter(xx,yy);
shg;


#### 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=-.5;
xMax=.5;
yMin=-.5;
yMax=.5;

#Parameters for the parent and daughter point processes
lambdaParent=10;#density of parent Poisson point process
lambdaDaughter=100;#mean number of points in each cluster

#Extended simulation windows parameters
xMinExt=xMin-rExt;
xMaxExt=xMax+rExt;
yMinExt=yMin-rExt;
yMaxExt=yMax+rExt;
#rectangle dimensions
xDeltaExt=xMaxExt-xMinExt;
yDeltaExt=yMaxExt-yMinExt;
areaTotalExt=xDeltaExt*yDeltaExt; #area of extended rectangle

#Simulate Poisson point process for the parents
numbPointsParent=rpois(1,areaTotalExt*lambdaParent);#Poisson number of points
#x and y coordinates of Poisson points for the parent
xxParent=xMinExt+xDeltaExt*runif(numbPointsParent);
yyParent=yMinExt+yDeltaExt*runif(numbPointsParent);

#Simulate Poisson point process for the daughters (ie final poiint process)
numbPoints=sum(numbPointsDaughter); #total number of points

#Generate the (relative) locations in polar coordinates by
#simulating independent variables.
theta=2*pi*runif(numbPoints); #angular coordinates

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

#replicate parent points (ie centres of disks/clusters)
xx=rep(xxParent,numbPointsDaughter);
yy=rep(yyParent,numbPointsDaughter);

#translate points (ie parents points are the centres of cluster disks)
xx=xx+xx0;
yy=yy+yy0;

#thin points if outside the simulation window
booleInside = ((xx >= xMin) & (xx <= xMax) & (yy >= yMin) & (yy <= yMax));
#retain points inside simulation window
xx=xx[booleInside];
yy=yy[booleInside];

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


Of course, as I have mentioned before, simulating a spatial point processes in R is even easier with the powerful spatial statistics library spatstat. The Matérn cluster point process is simulated by using the function rMatClust, but other cluster point processes, including Neyman-Scott types, are possible.

library(spatstat);
#Simulation window parameters
xMin=-.5;xMax=.5; yMin=-.5;yMax=.5;

#Parameters for the parent and daughter point processes lambdaParent=10;#intensity of parent Poisson point
processlambdaDaughter=100; #mean number of points in each cluster

#Simulate and plot (in one line) Matern cluster point process.

##### Python

Note: in previous posts I used the SciPy functions for random number generation, but now use the NumPy ones, but there is little difference, as SciPy builds off NumPy.

import numpy as np;  # NumPy package for arrays, random number generation, etc
import matplotlib.pyplot as plt  # For plotting
# Simulation window parameters
xMin = -.5;
xMax = .5;
yMin = -.5;
yMax = .5;

# Parameters for the parent and daughter point processes
lambdaParent = 10;  # density of parent Poisson point process
lambdaDaughter = 100;  # mean number of points in each cluster

# Extended simulation windows parameters
xMinExt = xMin - rExt;
xMaxExt = xMax + rExt;
yMinExt = yMin - rExt;
yMaxExt = yMax + rExt;
# rectangle dimensions
xDeltaExt = xMaxExt - xMinExt;
yDeltaExt = yMaxExt - yMinExt;
areaTotalExt = xDeltaExt * yDeltaExt;  # area of extended rectangle

# Simulate Poisson point process for the parents
numbPointsParent = np.random.poisson(areaTotalExt * lambdaParent);  # Poisson number of points
# x and y coordinates of Poisson points for the parent
xxParent = xMinExt + xDeltaExt * np.random.uniform(0, 1, numbPointsParent);
yyParent = yMinExt + yDeltaExt * np.random.uniform(0, 1, numbPointsParent);

# Simulate Poisson point process for the daughters (ie final poiint process)
numbPoints = sum(numbPointsDaughter);  # total number of points

# Generate the (relative) locations in polar coordinates by
# simulating independent variables.
theta = 2 * np.pi * np.random.uniform(0, 1, numbPoints);  # angular coordinates

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

# replicate parent points (ie centres of disks/clusters)
xx = np.repeat(xxParent, numbPointsDaughter);
yy = np.repeat(yyParent, numbPointsDaughter);

# translate points (ie parents points are the centres of cluster disks)
xx = xx + xx0;
yy = yy + yy0;

# thin points if outside the simulation window
booleInside = ((xx >= xMin) & (xx <= xMax) & (yy >= yMin) & (yy <= yMax));
# retain points inside simulation window
xx = xx[booleInside];
yy = yy[booleInside];

# 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.