## Cox point process

In previous posts I have stressed the importance of the Poisson point process, but it can be unsuitable for certain mathematical models.  We can generalize this point process by first considering a non-negative random measure, called a driving or directing measure. Then a Poisson point process, which is independent of the random driving measure, is generated by using the random measure as its intensity or mean measure. This doubly stochastic construction gives what is called a Cox point process.

Typically, we don’t observe the driving measure, so it’s impossible to distinguish a Cox point process from a Poisson point if there’s only one realization available. The properties of a Cox point process are derived by conditioning on the random driving measure, and then using the properties of the Poisson point process.

Cox point processes are also known as doubly stochastic Poisson point processes. Guttorp and Thorarinsdottir argue that it should be called the Quenouille point process, as Maurice Quenouille introduced an example of it before Sir David Cox, but I opt for the more common name.

In this post I’ll cover a couple examples of Cox point processes, but first I’ll give a more precise mathematical definition.

## Definition

We consider a point process defined on some underlying mathematical space $$\mathbb{S}$$, which is sometimes called the carrier space or state space.  For applications, the underlying space is often the real line $$\mathbb{R}$$, the plane $$\mathbb{R}^2$$, or some other familiar mathematical space like a square lattice.

For the first definition, we use the concept of a random measure.

Let $$M$$ be a non-negative random measure on $$\mathbb{S}$$. Then a point process $$\Phi$$ defined on some underlying space $$\mathbb{S}$$ is a Cox point process driven by the intensity measure $$M$$ if the conditional distribution of $$\Phi$$ is a Poisson point process with intensity function $$M$$.

Alternatively, we can give a slightly less general definition of a Cox  point process by using a random intensity function.

Let $$Z=\{Z(x):x\in\mathbb{S} \}$$ be a non-negative random field such that with probability one, $$x\rightarrow Z(x)$$ is a locally integrable function. Then a point process $$\Phi$$ defined on some underlying space $$\mathbb{S}$$ is a Cox point process driven by $$Z$$ if the conditional distribution of $$\Phi$$ is a Poisson point process with intensity function $$Z$$.

The random driving measure $$M$$ is then simply the integral
$$M(B)=\int_B Z(x)\, dx , \quad B\subseteq S.$$

## Over-dispersion

We will soon see that the random driving measures take different forms, giving different Cox point processes. But there is a general observation that can be made for all Cox point processes. For any region $$B \subseteq S$$, it can be shown that the number of points $$\Phi (B)$$ adheres to the inequality
$$\mathbb{Var} [\Phi (B)] \geq \mathbb{E} [\Phi (B)],$$

where $$\mathbb{Var} [\Phi (B)]$$ is the variance of the random variable $$\Phi (B)$$.  As a comparison, for a Poisson point process $$\Phi’$$, the variance of $$\Phi’ (B)$$ is simply $$\mathbb{Var} [\Phi’ (B)] =\mathbb{E} [\Phi’ (B)]$$.  Due to its greater variance, the Cox point process is said to be over-dispersed compared to the Poisson point process.

## Special cases

There is an virtually unlimited number of ways to define a random driving measure, where each one yields a different a Cox point process. But in general mathematicians are and practitioners are restricted by examining only tractable and interesting Cox point processes. I will give some common examples of such Cox point processes, but I stress that the Cox point process family is very large.

### Mixed Poisson point process

For the random driving measure $$M$$, an obvious example is the product form $$M= Y \mu$$,  where $$Y$$ is some independent non-negative random variable and $$\mu$$ is the Lebesgue measure on $$\mathbb{S}$$, meaning $$Y$$ is the only source of randomness. This driving measure gives the mixed Poisson point process.

### Log-Gaussian Cox point process

Instead of a random variable, we can use a non-negative random field to define a random driving measure.  We then have the product form $$M= Y \mu$$, where $$Y$$ is now some independent non-negative random field. (A random field is a collection of random variables indexed by some set, which in this case is the underlying space $$\mathbb{S}$$.)

Arguably the most tractable and used random field is the Gaussian random field. This random field, like Gaussian or normal random variables, takes negative and positive values, but if we define the random field such that its logarithm is a Gaussian field $$Z$$, then we obtain the non-negative random driving measure $$M=\mu e^Z$$, which gives the log-Gaussian Cox point process.

This point process has found applications in spatial statistics.

### Cox-Poisson line-point process

To construct this Cox point process, we first consider a Poisson line process, which I discussed previously.  Given a Poisson line process, we then place an independent one-dimensional Poisson point process on each line, then we obtain an example of a Cox point process, which we could call a Cox line-point process orCox-Poisson line-point process, but I am not sure of the name.

Researchers have recently used this point process to study wireless communication networks in cities, where the streets correspond to Poisson lines. For example, see these two preprints:

### Shot-noise Cox point process

We construct the next Cox point process by first considering a Poisson point process on the space $$\mathbb{S}$$ to create a shot noise term, which we then use as the driving measure of the Cox point process.

More specifically, we first introduce a kernel function $$k(\cdot,\cdot)$$ on $$\mathbb{S}$$, where $$k(x,\cdot)$$ is a probability density function for all points $$x\in \mathbb{S}$$. We then consider a Poisson point process $$\Phi’$$ on $$\mathbb{S}\times (0,\infty)$$, which we assume has a locally integrable intensity function $$\mu$$. (We can interpret  the point process $$\Phi’$$ as a spatially-dependent marked Poisson point process,  where the unmarked Poisson point process is defined on $$\mathbb{S}$$, and each point $$X$$ of this unmarked point process has a mark  $$T \in (0,\infty)$$ with probability density $$\mu(X,t)$$.) The resulting shot noise

$$Z(x)= \sum_{(Y,T)\in \Phi’} T \, k(Y,x)\,.$$

gives the random field, which we then use as the random intensity function to drive the shot-noise Cox point process.

In previous posts, I have detailed how to simulate non-Poisson point processes such as the Matérn and Thomas cluster point processes. These are, more specifically, examples of a Neyman-Scott point process, which is a special case of a shot noise Cox point process. All these point processes find applications in spatial statistics.

## Simulation

Unfortunately, there is no universal way to simulate all Cox point processes — and even if there were one, it would not be the most optimal way for every Cox point process. The simulation method depends on how the Cox point process is constructed, which usually means how its directing or driving measure is defined.

In previous posts I have presented ways (with code) to simulate the following Cox point processes:

In addition to the Matérn and Thomas point processes, there are ways to simulate to more general shot noise Cox point processes, which I will cover in another post.

For general Cox point processes, I suggest the following: Chapter 6 in the monograph Poisson processes by Kingman; Chapter 5 in Statistical Inference and Simulation for Spatial Point Processes by Møller and Waagepetersen; and Section 5.2 in Stochastic Geometry and its Applications by Chiu, Stoyan, Kendall and Mecke. For a much more mathematical treatment, see Chapter 13 in Lectures on the Poisson Process by Last and Penrose. Grandell wrote two detailed monographs titled Mixed Poisson Process and Doubly Stochastic Poisson Processes.

Motivated by applications in spatial statistics, Jesper Møller has (co)-written papers on specific Cox point processes. For example:

• 2001, Møller, Syversveen, and Waagepetersen, Log Gaussian Cox Processes;
• 2003, Møller, Shot noise Cox Processes;
• 2005, Møller and Torrisi,Generalised shot noise Cox processes.

I also suggest the survey article:

• 2003, Møller and Waagepetersen, Modern statistics for spatial point processes.

## Signal strengths of a wireless network

In two previous posts, here and here, I discussed the importance of the quantity called the signal-to-interference ratio, which is usually abbreviated as SIR, for studying communication in wireless networks. In everyday terms, for a listener to hear a certain speaker in a room full of people speaking, the ratio of the speaker’s volume to the sum of the volumes of everyone else heard by the listener. The SIR is the communication bottleneck for any receiver and transmitter pair in a wireless network.

But the strengths (or power values) of the signals are of course also important. In this post I will detail how we can model them using a a simple network model with a single observer.

### Propagation model

For a transmitter located at $$X_i\in \mathbb{R}^2$$, researchers usually attempt to represent the received power of the signal $$P_i$$ with a propagation model. Assuming the power is received at $$x\in \mathbb{R}^2$$, this mathematical model consists of a random and a deterministic component taking the general form
$$P_i(x)=F_i\,\ell(|X_i-x|) ,$$
where $$\ell(r)$$ is a non-negative function in $$r>0$$ and $$F_i$$ is a non-negative random variable.

The function $$\ell(r)$$ is called the pathloss function, and common choices include $$\ell(r)=(\kappa r)^{-\beta}$$ and $$\ell(r)=\kappa e^{-\beta r}$$, where $$\beta>0$$ and $$\kappa>0$$ are model constants.

The random variables $$F_i$$ represent signal phenomena such as multi-path fading and shadowing (also called shadow fading), caused by the signal interacting with the physical environment such as buildings. It is often called fading or shadowing variables.

We assume the transmitters locations $$X_1,\dots,X_n$$ are on the plane $$\mathbb{R}^2$$. Researchers typically assume they form a random point process or, more precisely, the realization of a random point process.

## From two dimensions to one dimension

For studying wireless networks, a popular technique is to consider a wireless network from the perspective of a single observer or user. Researchers then consider the incoming or received signals from the entire network at the location of this observer or user. They do this by considering the inverses of the signal strengths, namely

$$L_i(x): = \frac{1}{P_i}=\frac{1}{F_i \,\ell(|X_i-x|) }.$$

Mathematically, this random function is simply a mapping from the two-dimensional plane $$\mathbb{R}^2$$ to the one-dimensional non-negative real line $$\mathbb{R}_0^+=[0,\infty)$$.

If the transmitters are located according to a non-random point pattern or a random point process, this random mapping generates a random point process on the non-negative real line. The resulting one-dimensional point process of the values $$L_1,L_2,\dots,$$ has been called (independently) propagation (loss) process or path loss (with fading) process. More recently, my co-authors and I decided to call it a projection process, but of course the precise name doesn’t mattter

## Intensity measure of signal strengths

Assuming a continuous monotonic path loss function $$\ell$$ and the fading variables $$F_1, F_2\dots$$ are iid, if the transmitters form a stationary random point process with intensity $$\lambda$$, then the inverse signal strengths $$L_1,L_2,\dots$$ form a random point process on the non-negative real line with the intensity measure $$M$$.

$$M(t) =\lambda \pi \mathbb{E}( [\ell(t F)^{-1} ]^2)\,,$$

where $$\ell^{-1}$$ is the generalized inverse of the function $$\ell$$. This expression can be generalized for a non-stationary point process with general intensity measure $$\Lambda$$.

The inverses $$1/L_1,1/L_2,\dots$$, which are the signal strengths, forprocess with intensity measure

$$\bar{M}(s) =\lambda \pi \mathbb{E}( [\ell( F/s)^{-1} ]^2).$$

## Poisson transmitters gives Poisson signal strengths

Assuming a continuous monotonic path loss function $$\ell$$ and the fading variables $$F_1, F_2\dots$$ are iid, if the transmitters form a Poisson point process with intensity $$\lambda$$, then the inverse signal strengths $$L_1,L_2,\dots$$ form a Poisson point process on the non-negative real line with the intensity measure $$M$$.

If $$L_1,L_2,\dots$$ form a homogeneous Poisson point process, then the inverses $$1/L_1,1/L_2,\dots$$ will also form a Poisson point process with intensity measure $$\bar{M}(s) =\lambda \pi \mathbb{E}( [\ell( F/s)^{-1} ]^2).$$

## Propagation invariance

For $$\ell(r)=(\kappa r)^{-\beta}$$ , the expression for the intensity measure $$M$$ reduces to
$$M(t) = \lambda \pi t^{-2/\beta} \mathbb{E}( F^{-2/\beta})/\kappa^2.$$

What’s striking here is that information of the fading variable $$F$$ is captured simply by one moment $$\mathbb{E}( F^{-2/\beta})$$. This means that two different distributions will give the same results as long as this moment is matching. My co-authors and I have been called this observation propagation invariance.

## Some history

To study just the (inverse) signal strengths as a point process on the non-negative real line was a very useful insight. It was made independently in these two papers:

• 2008, Haenggi, A geometric interpretation of fading in wireless
networks: Theory and applications;
• 2010, Błaszczyszyn, Karray, and Klepper, Impact of the geometry, path-loss exponent and random shadowing on the mean interference factor in wireless cellular networks.

My co-authors and I presented a general expression for the intensity measure $$M$$ in the paper:

• 2018, Keeler, Ross and Xia, When do wireless network signals appear Poisson?.

This paper is also contains examples of various network models.

A good starting point on this topic is the Wikipedia article Stochastic geometry models of wireless networks. The paper that my co-authors and I wrote has details on the projection process.

With Bartek Błaszczyszyn, Sayan Mukherjee, and Martin Haenggi, I co-wrote a short monograph on SINR models called Stochastic Geometry Analysis of Cellular Networks, which is written at a slightly more advanced level. The book put an emphasis on studying the point process formed from inverse signal strengths, we call the projection process.

## The Standard Model of wireless networks

In the previous post I discussed the signal-to-interference-plus ratio or SIR in wireless networks. If noise is included, then then signal-to-interference-plus-noise ratio or just SINR. But I will just write about SIR, as most results that hold for SIR, will also hold for SINR without any great mathematical difficulty.

The SIR is an important quantity due to reasons coming from information theory.  If you’re unfamiliar  with it, I suggest reading the previous post.

In this post, I will describe a very popular mathematical model of the SIR, which I like to call the standard model. (This is not a term used in the literature as I have borrowed it from physics.)

## Definition of SIR

To define the SIR, we consider a wireless network of $$n$$ transmitters with positions located at $$X_1,\dots,X_n$$ in some region of space. At some location $$x$$, we write $$P_i(x)$$ to denote the power value of a signal received at $$x$$ from transmitter  $$X_i$$. Then at location $$x$$, the SIR with respect to transmitter $$X_i$$ is
$$\text{SIR}(x,X_i) := \frac{P_i(x)}{\sum\limits_{j\neq i} P_j(x)} .$$

Researchers usually attempt to represent the received power of the signal $$P_i(x)$$ with a propagation model. This mathematical model  consists of a random and a deterministic component given by
$$P_i(x)=F_i\ell(|X_i-x|) ,$$
where $$\ell(r)$$ is a non-negative function in $$r\geq 0$$ and $$F_i$$ is a non-negative random variable. The function $$\ell(r)$$  is often called the path loss function. The random variables represent random fading or shadowing.

## Standard model

Based on the three model components of fading, path loss, and transmitter locations, there are many combinations possible. That said, researchers generally (I would guess, say, 90 percent or more) use a single combination, which I call the standard model.

The three standard model assumptions are:

1. Singular power law path loss $$\ell(r)=(\kappa r)^{-\beta}$$.
2. Exponential distribution for fading variables, which are independent and identically distributed (iid).
3. Poisson point process for transmitter locations.

Why these three? Well, in short, because they work very well together. Incredibly, it’s sometimes possible to get relatively a simple  mathematical expression for, say, the coverage probability $$\mathbb{P}[\text{SIR}(x,X_i)>\tau ]$$, where $$\tau>0$$.

I’ll now detail the reasons more specifically.

### Path loss

The $$\ell(r)=(\kappa r)^{-\beta}$$ is very simple, despite having a singularity at $$r=0$$. This allows simple algebraic manipulation of equations.

Some, such as myself, are initially skeptical of this function as it gives an infinitely strong signal at the transmitter due to the singularity in the function $$\ell(r)=(\kappa r)^{-\beta}$$. More specifically, the path loss of the signal from transmitter $$X_i$$ approaches infinity as $$x$$ approaches $$X_i$$ .

But apparently, overall, the singularity does not have a significant impact on most mathematical results, at least qualitatively. That said, one still observe consequences of this somewhat physically unrealistic model assumption. And I strongly doubt enough care is taken by researchers to observe and note this.

Interestingly, the original reason why exponential variables were used is because it allowed the SIR problem to be reformulated into a problem of a Laplace transform of a random variable, which for a random variable $$Y$$ is defined as

$$\mathcal{L}_Y(t)=\mathbb{E}(e^{- Y t}) \, .$$

where $$t\geq 0$$. (This is essentially the moment-generating function with $$-t$$ instead of $$t$$.)

The reason for this connection is that the tail distribution of an exponential variable $$F$$ with mean $$\mu$$  is simply $$\mathbb{P}(F>t)= e^{-t/\mu}$$.  In short, with the exponential assumption, various conditioning arguments eventually lead to Laplace transforms of random variables.

### Transmitters locations

No prizes for guessing that researcher overwhelmingly use a (homogeneous) Poisson point process for the transmitter (or receiver) locations. When developing mathematical models with point processes, if you can’t get any results with the Poisson point process, then abandon all hope.

It’s the easier to work with this point process due to its independence property, which leads to another useful property. For Poisson point process, the Palm distribution is known, which is the distribution of a point process conditioned on a point (or collection of points) existing in a specific location of the underlying space on which the point process is defined.  In general, the Palm distribution is not known for many point processes.

## Random propagation effects can lead to Poisson

A lesser known reason why researchers would use the Poisson point process is that, from the perspective of a single observer in the network, it can be used to capture the randomness in the signal strengths.  Poisson approximation results in probability imply that randomly perturbing the signal strengths can make signals appear more Poisson, by which I mean  the signal strengths behave stochastically or statistically as though they were created by a Poisson network of transmitters.

The end result is that a non-Poisson network can appear more Poisson, even if the transmitters do not resemble (the realization of) a Poisson point process. The source of randomness that makes a non-Poisson network appear more Poisson is the random propagation effects of fading, shadowing, randomly varying antenna gains, and so on, or some combination of these.

A good starting point on this topic is the Wikipedia article Stochastic geometry models of wireless networks. This paper is also good:

• 2009, Haenggi, Andrews, Baccelli, Dousse, Franceschetti, Stochastic Geometry and Random Graphs for the Analysis and Design of Wireless Networks.

This paper by my co-authors and I has some details on standard model and why a general network model behaving Poisson in terms of the signal strengths:

• 2018, Keeler, Ross and Xia, When do wireless network signals appear Poisson?.

Early books on the subject include the two-volume textbooks  Stochastic Geometry and Wireless Networks by François Baccelli and Bartek Błaszczyszyn, where the first volume is on theory and the second volume is on applications.  Martin Haenggi wrote a very readable introductory book called Stochastic Geometry for Wireless networks.

Finally, I co-wrote with Bartek Błaszczyszyn, Sayan Mukherjee, and Martin Haenggi a short monograph on SINR models called Stochastic Geometry Analysis of Cellular Networks, which is written at a slightly more advanced level. This book has a section on why signal strengths appear Poisson.

## Signal-to-interference ratio in wireless networks

The fundamentals of information theory say that to successfully communicate across any potential communication link the signal strength of the communication must be stronger than that of the back ground noise, which leads to the fundamental quantity known as signal-to-noise ratio. Information theory holds in very general (or, in mathematical speak, abstract) settings. The communication could be, for example, a phone call on an old wired landline, two people talking in a bar, or a hand-written letter, for which the respective signals in these examples are the electrical current, speaker’s voice, and the writing. (Respective examples of noise could be, for example, thermal noise in the wires, loud music, or coffee stains on the letter.)

In wireless networks, it’s possible for a receiver to simultaneously detect signals from multiple transmitters, but the receiver typically only wants to receive one signal. The other unwanted or interfering signals form a type of noise, which is usually called interference, and the other (interfering) transmitters are called interferers. Consequently, researchers working on wireless networks study the signal-to-interference ratio, which is usually abbreviated as SIR. Another name for the SIR is carrier-to-interference ratio.

If we also include background noise, which is coming not from the interferers, then the quantity becomes the signal-to-interference-plus-noise ratio or just SINR. But I will just write about SIR, though jumping from SIR to SINR is usually not difficult mathematically.

The concept of SIR makes successful communication more difficult to model and predict, as it just doesn’t depend on the distance of the communication link. Putting the concept in everyday terms, for a listener to hear a certain speaker in a room full of people all speaking to the listener, it is not simply the distance to the speaker, but rather the ratio of the speaker’s volume to the sum of the volumes of everyone else heard by the listener. The SIR is the communication bottleneck for any receiver and transmitter pair in a wireless network.

In wireless network research, much work has been done to examine and understand communication success in terms of interference and SIR, which has led to a popular mathematical model that incorporates how signals propagate and the locations of transmitters and receivers.

## Definition

To define the SIR, we consider a wireless network of transmitters with positions located at $$X_1,\dots,X_n$$ in some region of space. At some location $$x$$, we write $$P_i(x)$$ to denote the power value of a signal received at $$x$$ from transmitter $$X_i$$. Then at location $$x$$, the SIR with respect to transmitter $$X_i$$ is
$$\text{SIR}(x,X_i) :=\frac{P_i(x)}{\sum\limits_{j\neq i} P_j(x)} =\frac{P_i(x)}{\sum\limits_{j=1}^{n} P_j(x)-P_i(x)} .$$

The numerator is the signal and the denominator is the interference.  This ratio tells us that increasing the number of transmitters $$n$$ decreases the original SIR values. But then, in exchange, there is a greater number of transmitters for the receiver to connect to, some of which may have larger $$P_i(x)$$ values and, subsequently, SIR values. This delicate trade-off makes it challenging and interesting to mathematically analyze and design networks that deliver high SIR values.

Researchers usually assume that the SIR is random. A quantity of interest is the tail distribution of the SIR, namely

$$\mathbb{P}[\text{SIR}(x,X_i)>\tau ] := \frac{P_i(x)}{\sum\limits_{j\neq i} P_j(x)} \,,$$

where $$\tau>0$$ is some parameter, sometimes called the SIR threshold. For a given value of $$\tau$$, the probability $$\mathbb{P}[\text{SIR}(x,X_i)>\tau]$$ is sometimes called the coverage probability, which is simply the probability that a signal coming from $$X_i$$ can be received successfully at location $$x$$.

## Mathematical models

### Propagation

Researchers usually attempt to represent the received power of the signal $$P_i(x)$$ with a propagation model. This mathematical model consists of a random and a deterministic component taking the general form
$$P_i(x)=F_i\ell(|X_i-x|) ,$$
where $$F_i$$ is a non-negative random variable and $$\ell(r)$$ is a non-negative function in $$r \geq 0$$.

##### Path loss

The function $$\ell(r)$$ is called the path loss function, and common choices include $$\ell(r)=(\kappa r)^{-\beta}$$ and $$\ell(r)=\kappa e^{-\beta r}$$, where $$\beta>0$$ and $$\kappa>0$$ are model constants, which need to be fitted to (or estimated with) real world data.

Researchers generally assume that the so-called path loss function $$\ell(r)$$ is decreasing in $$r$$, but actual path loss (that is, the change in signal strength over a path travelled) typically increases with distance $$r$$. Researchers originally assumed path loss functions to be increasing, not decreasing, giving the alternative (but equivalent) propagation model
$$P_i(x)= F_i/\ell(|X_i-x|).$$

But nowadays researchers assume that the function $$\ell(r)$$ is decreasing in $$r$$. (Although, based on personal experience, there is still some disagreement on the convention.)

With the random variable $$F_i$$, researchers seek to represent signal phenomena such as multi-path fading and shadowing (also called shadow fading), caused by the signal interacting with the physical environment such as buildings. These variables are often called fading or shadowing variables, depending on what physical phenomena they are representing.

Typical distributions for fading variables include the exponential and gamma distributions, while the log-normal distribution is usually used for shadowing. The entire collection of fading or shadowing variables is nearly always assumed to be independent and identically distributed (iid), but very occasionally random fields are used to include a degree of statistical dependence between variables.

### Transmitters locations

In general, we assume the transmitters locations $$X_1,\dots,X_n$$ are on the plane $$\mathbb{R}^2$$. To model interference, researchers initially proposed non-random models, but they were considered inaccurate and intractable. Now researchers typically use random point processes or, more precisely, the realizations of random point processes for the transmitter locations.

Not surprisingly, the first natural choice is the Poisson point process. Other point processes have been used such as Matérn and Thomas cluster point processes, and Matérn hard-core point processes, as well as determinantal point processes, which I’ll discuss in another post.

## Some history

Early random models of wireless networks go back to the 60s and 70s, but these were based simply on geometry: meaning a transmitter could communicate successfully to a receiver if they were closer than some fixed distance. Edgar Gilbert created the field of continuum percolation with this significant paper:

• 1961, Gilbert, Random plane networks.

Interest in random geometrical models of wireless networks continued into the 70s and 80s. But there was no SIR in these models.

Motivated by understanding SIR, researchers in the late 1990s and early 2000s started tackling SIR problems by using a random model based on techniques from stochastic geometry and point processes. Early papers include:

• 1997, Baccelli, Klein, Lebourges ,and Zuyev, Stochastic geometry and architecture of communication networks;
• 2003, Baccelli and Błaszczyszyn , On a coverage process ranging from the Boolean model to the Poisson Voronoi tessellation, with applications to wireless communications;
• 2006, Baccelli, Mühlethaler, and Błaszczyszyn, An Aloha protocol for multihop mobile wireless networks.

But they didn’t know that some of their results had already been discovered independently by researchers working on wireless networks in the early 1990s. These papers include:

• 1994, Pupolin and Zorzi, Outage probability in multiple access packet radio networks in the presence of fading;
• 1990, Sousa and Silvester, Optimum transmission ranges in a direct-sequence spread-spectrum multihop packet radio network.

The early work focused more on small-scale networks like wireless ad hoc networks. Then the focus shifted dramatically to mobile or cellular phone networks with the publication of the paper:

• 2011, Andrews, Baccelli, Ganti, A tractable approach to coverage and rate in cellular networks.

It’s can be said with confidence that this paper inspired much of the interest in using point processes to develop models of wireless networks. The work generally considers the SINR in the downlink channel for which the incoming signals originate from the phone base stations.

A good starting point on this topic is the Wikipedia article Stochastic geometry models of wireless networks. This paper is also good:

• 2009, Haenggi, Andrews, Baccelli, Dousse, Franceschetti, Stochastic Geometry and Random Graphs for the Analysis and Design of Wireless Networks.

Early books on the subject include the two-volume textbooks Stochastic Geometry and Wireless Networks by François Baccelli and Bartek Błaszczyszyn, where the first volume is on theory and the second volume is on applications. Martin Haenggi wrote a very readable introductory book called Stochastic Geometry for Wireless networks.

Finally, Bartek Błaszczyszyn, Sayan Mukherjee, Martin Haenggi, and I wrote a short book on SINR models called Stochastic Geometry Analysis of Cellular Networks, which is written at a slightly more advanced level. The book put an emphasis on studying the point process formed from inverse signal strengths, we call the projection process.

## Simulating Matérn hard-core point processes

If you wanted to create a point process with repulsion, a reasonable first attempt would be to build off a Poisson point process by removing points according to some rule to ensure that no two points were within a certain distance of each other. Using this natural idea, Bertril Matérn proposed a family of repulsive point processes called Matérn hard-core point processes.

More specifically, Matérn proposed several points processes, including two types of hard-core point processes now called Type I and Type II. (Matérn proposed a third type, called Type III, but it’s considerably harder to simulate on a computer, as detailed in this article.) These types of hard-core point processes are completely different to the Matérn cluster point process.

As I discussed in a previous post, the Poisson point process may not be adequate for representing point phenomena whose points exhibit large degrees of repulsion or clustering. I already covered the Matérn and Thomas cluster point processes, which show distinct clustering in their configurations. In this post, I’ll cover Matérn hard-core point processes. The Type I point processes is the easier of the two, so I’ll start with that one.

## Overview

Simulating Matérn hard-core point processes 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. I have already written about simulating the homogeneous Poisson point processes on a rectangle and a disk, so those posts are good starting points.

Given the Poisson point process, the points then need to be thinned in such a manner to ensure that for each point, there is no other point within some fixed $$r>0$$ of the point. This distance $$r>0$$ is the radius of the hard core of each point.

I have already covered the point process operation of thinning. But it’s important to note here that in this construction a dependent thinning is being applied. (If I just applied an independent thinning, then the resulting point process will be another Poisson point process with no repulsion between points.)

## Edge effects

The main trick behind sampling this point process is that it’s possible for points inside the simulation window to be thinned due to their closeness to points that are located outside the simulation window. In other words, points outside the simulation window can cause points inside the window to be thinned. (I discussed a very similar issue in the posts on the Matérn and Thomas cluster point processes.)

To remove these edge effects, the underlying Poisson point process must be simulated on an extended version of the simulation window. The points are then thinned according to a dependent thinning, which is covered in the next section. Then only the retained points inside the simulation window are kept and the remaining points are ignored. Consequently, the underling Poisson points are simulated on an extended window, but we only see the final points inside the simulation window.

To create the extended simulation window, we 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 another point (outside the simulation window) can exist, while still causing points inside the simulation window to be thinned. This means it is impossible for a hypothetical point beyond this distance (outside the extended window) to cause a point inside the simulation window to be thinned.

## Dependent thinning rules

##### Type I

For each point inside the simulation window, check if there are any other points (including those in the extended window) within distance $$r$$ of the point. If no, then keep the point. If yes, then remove the point and the points that are within distance $$r$$ of the point. The remaining points inside the simulation window form a Matérn Type I point process.

This is a relatively simple thinning rule, which only requires calculating all the inter-point distances. But it is also a very strong thinning rule, meaning that it removes many points. Depending on the Poisson point process intensity $$\lambda$$ and core radius $$r$$, it is quite possible that all the points are removed, resulting in an empty configuration.

Now we examine the case when the thinning rule is not as strong.

##### Type II

To create Matérn Type II point process, we assign an independent uniform random variable to each point of the underlying Poisson point process defined on the extended window. In point process terminology, these random variables are called marks, resulting in a marked point process. In the the context of the Matérn Type II point process, these random random marks are usually called ages.

Then for each point in the simulation window, we consider all the points within distance $$r$$ of the point. If this point is the youngest (or, equivalently, the oldest) point, then the point is kept. In other words, the point is only kept if its random mark is smaller (or larger) than the random marks of all the other points within distance $$r$$ of the point. The remaining points inside the simulation window form a Matérn Type II point process.

## Intensity expressions

Using point process and probability theory, one can derive mathematical expressions for the intensities (that is, the average density of points per unit area). These closed-form expressions can then be used to check that the correct number of points are being generated on average over many simulations.

##### Type I

The intensity of the Type I point process is given by

$\mu_1=\lambda e^{-\lambda \pi r^2},$

where $$\lambda \pi r^2$$ is simply the area of the core.

##### Type II

The intensity of the Type II point process is given by

$\mu_2=\frac{1}{\pi r^2}(1-e^{-\lambda \pi r^2}),$

which can be written with the intensity of the the Type I point process as

$\mu_2=\frac{1}{\pi r^2}(1-\frac{\mu_1}{\lambda}).$

## Code

I wrote the sampling code in MATLAB and Python, which are, as usual, very similar to each other. The code, which is is located here, simulates both Type I and II Matérn points processes. It also compares the empirical intensity to the the values given by the mathematical expressions in the previous section.

##### MATLAB

The MATLAB code is here.

##### Python

The Python code is here.

## Results

I have plotted single realizations of the Matern Type I and II point processes, as well as the underlying Poisson point process in the same window.

##### Python

Matérn hard-core point processes are covered in standard books on the related fields of spatial statistics, point processes and stochastic geometry, such as the following: Spatial Point Patterns: Methodology and Applications with R by Baddeley, Rubak and Turner, on page 140; Statistical Analysis and Modelling of Spatial Point Patterns Statistics by Illian, Penttinen, Stoyan, amd Stoyan, Section 6.5.2, starting on page 388; and; Stochastic Geometry and its Applications by Chiu, Stoyan, Kendall and Mecke, Section 5.4, starting on page 176. The first two books are particularly good for beginners.

The book Spatial Point Patterns: Methodology and Applications with R is written by spatial statistics experts Baddeley, Rubak and Turner. It covers the spatial statistics (and point process simulation) R-package spatstat., which has the functions rMaternI and rMaternII for simulating the two point processes respectively.

## Simulating Poisson point processes faster

As an experiment, I tried to write code for simulating many realizations of a homogeneous Poisson point process in a very fast fashion. My idea was to simulate all the realizations in two short steps.

In reality, the findings of this experiment and the contents of this post have little practical value, as computers are so fast at generating Poisson point processes. Still, it was an interesting task, which taught me a couple of things. And I did produce faster code.

### MATLAB

I first tried this experiment in MATLAB.

##### Vectorization

In languages like MATLAB, the trick for speed is to use vectorization, which means applying a single operation to an entire vector or matrix (or array) without using an explicit for-loop. Over the years, the people behind MATLAB has advised to use vectorization instead of for-loops, as for-loops considerably slowed down MATLAB code. (But, as as time goes by, it seems using for-loops in MATLAB doesn’t slow the code down as much as it used to.)

Simulating Poisson point processes is particularly amenable to vectorization, due to the independent nature of the point process. One can simulate the number of points in each realization for all realizations in one step. Then all the points across all realizations can also be positioned in one step. In the two-dimensional case, this results in two one-dimensional arrays (or vectors, in MATLAB parlance) for the $$x$$ and $$y$$ coordinates.  (Of course, in my code, I could have used just one two-dimensional array/vector  for the coordinates of the points, but I didn’t.)

After generating the points, the coordinates of the points need to be grouped into the different realizations and stored in appropriate data structures.

##### Data structures

The problem with storing point processes is that usually each realization has a different number of points, so more sophisticated data structures than regular arrays are needed. For MATLAB, each point process realization can be stored in a data object called a cell array.  For a structure array, it’s possible for each element (that is, structure) to be a different size, making them ideal for storing objects like point processes with randomly varying sizes.

In the case of two-dimensional point processes, two cell arrays  can be used to store the $$x$$ and $$y$$ coordinates of all the point process realizations. After randomly positioning all the points, they can be grouped into a cell array, where each cell array element represents a realization of the Poisson point process, by using the inbuilt function MATLAB mat2cell, which converts a matrix (or array) into a cell array.

Alternatively, we could use another MATLAB data object called a structure array. In MATLAB structures have fields, which can be, for example for a point process can be the locations of the points. Given cell arrays of equal size, we can convert them into a single structure array by using the inbuilt MATLAB function struct.

### Python

After successfully simulating Poisson point processes in MATLAB, I tried it in Python with the NumPy library.

##### Vectorization

I basically replicated what I did in MATLB using Python by positioning all the points in a single step. This gives two one-dimensional NumPy arrays for the $$x$$ and $$y$$ coordinates of all the point process realizations. (Again, I could have instead stored the coordinates as a single two-dimensional array, but I didn’t.)

Perhaps surprisingly, the vectorization step doesn’t speed things up much in Python with the NumPy library. This may be due to the fact that the underlying code is actually written in the C language.  That motivated me to see what methods have been implemented for simulating Poisson variables, which is the topic of the next couple posts.

##### Data structures

In Python, the data structure  known as a list is the natural choice. Similarly to cell arrays in MATLAB, two lists can be used for the $$x$$ and $$y$$ coordinates of all the points. Instead of MATLAB’s function mat2cell, I used the NumPy function numpy.split to create two lists from the two NumPy arrays containing the point coordinates.

Python does not seem to have an immediate equivalent to structure arrays in MATLAB. But in Python one can define a new data structure or class with fields, like a structure. Then one can create a list of those data structures with fields, which are called attribute references in Python.

### Code

The code in MATLAB and Python can be found here. For a comparison, I also generated the same number of point process realizations (without using vectorization) by using a trusty for-loop. The code compares the times of taken for implemented the two different approaches, labelled internally as Method A and Method B. There is a some time difference in the MATLAB code, but not much of a difference in the Python case.

I have commented out the sections that create data structures (with fields or attribute references) for storing all the point process realizations, but those sections should also work when uncommented.

## Simulating a Cox point process based on a Poisson line process

In the previous post, I described how to simulate a Poisson line process, which in turn was done by using insight from an earlier post on the Bertrand paradox.

Now, given a Poisson line process, for each line, if we generate an independent one-dimensional Poisson point point process on each line, then we obtain an example of a Cox point process. Cox point processes are also known as doubly stochastic Poisson point processes. On the topic of names, Guttorp and Thorarinsdottir argue that it should be called the Quenouille point process, as Maurice Quenouille introduced an example of it before Sir David Cox, but I opt for the more common name.

## Cox point proceesses

A Cox point process is a generalization of a Poisson point process. It is created by first considering a non-negative random measure, sometimes called a driving measure. Then a Poisson point process, which is independent of the random driving measure, is generated by using the random measure as its intensity or mean measure.

The driving measure of a Cox point process can be, for example, a non-negative random variable or field multiplied by a Lebesgue measure. In our case, the random measure is the underlying Poisson line process coupled with the Lebesgue measure on the line (that is, length).

Cox processes form a very large and general family of point processes, which exhibit clustering. In previous posts, I have covered two special cases of Cox point processes: the Matérn and Thomas cluster point processes. These are, more specifically, examples of a Neyman-Scott point process, which is a special case of a shot noise Cox point process. These two point processes are fairly easy to simulate, but that’s not the case for Cox point processes in general. Some are considerably easier than others.

## Motivation

I will focus on simulating the Cox point process formed from a Poisson line process with homogeneous Poisson point processes. I do this for two reasons. First, it’s easy to simulate, given we can simulate a Poisson line process. Second, it has been used and studied recently in the mathematics and engineering literature for investigating wireless communication networks in cities, where the streets correspond to Poisson lines; for example, see these two preprints:

Incidentally, I don’t know what to call this particular Cox point process. A Cox line-point process? A Cox-Poisson line-point process? But it doesn’t matter for simulation purposes.

## Method

We will simulate the Cox (-Poisson line-) point process on a disk. Why a disk? I suggest reading the previous posts on the Poisson line process the Bertrand paradox for why the disk is a natural simulation window for line processes.

Provided we can simulate a Poisson line process, the simulation method is quite straightforward, as I have essentially already described it.

##### Line process

First simulate a Poisson line process on a disk. We recall that for each line of the line process, we need to generate two independent random variables $$\Theta$$ and $$P$$ describing the position of the line. The first random variable $$\Theta$$ gives the line orientation, and it is a uniform random variable on the interval $$(0,2\pi)$$.

The second random variables $$P$$ gives the distance from the origin to the disk edge, and it is a uniform random variable on the interval $$(0,r)$$, where $$r$$ is the radius of the disk. The distance from the point $$(\Theta, P)$$ to the disk edge (that is, the circle) along the chord is:

$$Q=\sqrt{r^2-P^2}.$$

The endpoints of the chord (that is, the points on the disk edge) are then:

Point 1: $$X_1=P \cos \Theta+ Q\sin \Theta$$, $$Y_1= P \sin \Theta- Q\cos \Theta$$,

Point 2: $$X_2=P \cos \Theta- Q\sin \Theta$$, $$Y_2= P \sin \Theta+Q \cos \Theta$$.

The length of the line segment is $$2 Q$$. We can say this random line is described by the point $$(\Theta,P)$$.

##### One-dimensional Poisson point process

For each line (segment) in the line process, simulate a one-dimensional Poisson point process on it. Although I have never discussed how to simulate a one-dimensional (homogeneous) Poisson point process, it’s essentially one dimension less than simulating a homogeneous Poisson point process on a rectangle.

More specifically, given a line segment $$(\Theta,P)=(\theta,p)$$, you simulate a homogeneous Poisson point process with intensity $$\mu$$ on a line segment with length $$2 q$$, where $$q=\sqrt{r^2-p^2}$$. (I am now using lowercase letters to stress that the line is no longer random.) To simulate the homogeneous Poisson point process, you generate a Poisson random variable with parameter $$2 \mu q$$.

Now you need to place the points uniformly on the line segment. To do this, consider a single point on a single line. For this point, generate a single uniform variable $$U$$ on the interval $$(-1,1)$$. The tricky part is now getting the Cartesian coordinates right. But the above expressions for the endpoints suggest that the single random point has the Cartesian coordinates:

$$x=p \cos \theta+ U q\sin \theta$$, $$y=p \sin \theta- U q\cos \theta$$.

The two extreme cases of the uniform random variable $$U$$ (that is, $$U=-1$$ and $$U=1$$) correspond to the two endpoints of the line segment. We recall that $$Q$$ is the distance from the midpoint of the line segment to the disk edge along the line segment, so it makes sense that we want to vary this distance uniformly in order to uniformly place a point on the line segment. This uniform placement step is done for all the points of the homogeneous Point process on that line segment.

You repeat this procedure for every line segment. And that’s it: a Cox point process built upon a Poisson line process.

## Code

As always, the code from all my posts is online. For this post, I have written the code in MATLAB, R and Python.

For the first step, the reading material is basically the same as that for the Poisson line process, which overlaps with that of the Bertrand paradox. For the one-dimensional Poisson point process, we can use the reading material on the homogeneous Poisson point process on a rectangle.

For general Cox point processes, I recommend starting with the following: Chapter 6 in the monograph Poisson processes by Kingman; Chapter 5 in Statistical Inference and Simulation for Spatial Point Processes by Møller and Waagepetersen; and Section 5.2 in Stochastic Geometry and its Applications by Chiu, Stoyan, Kendall and Mecke. For a much more mathematical treatment, see Chapter 13 in Lectures on the Poisson Process by Last and Penrose.

For this particularly Cox point process, see the two aforementioned preprints, located here and here.

## Simulating a Poisson line process

Instead of points, we can consider other objects randomly scattered on some underlying mathematical space. If we take a Poisson point process, then we can use (infinitely long) straight lines instead of points, giving a Poisson line process. Researchers have studied and used this random object to model physical phenomena. In this post I’ll cover how to simulate a homogeneous Poisson line process in MATLAB, R and Python. The code which can be downloaded from here

## Overview

For simulating a Poisson line process, the key question is how to randomly position the lines.  This is related to a classic problem in probability theory called the Bertrand paradox.  I discussed this illustration in probability in a previous post, where I also included code for simulating it. I highly recommend reading that post first.

The Bertrand paradox involves three methods for positioning random lines. We use Method 2 to achieve a uniform positioning of lines, meaning the number of lines and orientation of lines is statistically uniform. Then it also makes sense to use this method for simulating a homogeneous (or uniform) Poisson line process.

We can interpret a line process as a point process. For a line process on the plane $$\textbf{R}^2$$, it can be described by a point process on $$(0,\infty)\times (0,2\pi)$$, which is an an infinitely long cylinder. In other words, the Poisson line process can be described as a Poisson point process.

For simulating a Poisson line process, it turns out the disk is the most natural setting. (Again, this goes back to the Bertrand paradox.) More specifically, how the (infinitely long) lines intersect a disk of a fixed radius $$r>0$$. The problem of simulating a Poisson line process reduces to randomly placing chords in a disk. For other simulation windows in the plane, we can always bound any non-circular region with a sufficiently large disk.

## Steps

##### Number of lines

Of course, with all things Poisson, the number of lines will be  a Poisson random variable, but what will its parameter be? This homogeneous (or uniform) Poisson line process forms a one-dimensional homogeneous (or uniform) Poisson point process around the edge of the disk with a circumference $$2 \pi r$$. Then the number of lines is simply a Poisson variable with parameter $$\lambda 2 \pi r$$.

##### Locations of points

To position a single line uniformly in a disk, we need to generate two uniform random variables. One random variable gives the angle describing orientation of the line, so it’s a uniform random variable $$\Theta$$ on the interval $$(0,2\pi)$$.

The other random variable gives the distance from the origin to the disk edge, meaning it’s a uniform random variable $$P$$ on the interval $$(0,r)$$, where $$r$$ is the radius of the disk.  The random radius and its perpendicular chord create a right-angle triangle.  The distance from the point $$(\Theta, P)$$ to the disk edge (that is, the circle) along the chord is:

$$Q=\sqrt{r^2-P^2}.$$

The endpoints of the chord (that is, the points on the disk edge) are then:

Point 1: $$X_1=P \cos \Theta+ Q\sin \Theta$$, $$Y_1= P \sin \Theta- Q\cos \Theta$$,

Point 2: $$X_2=P \cos \Theta- Q\sin \Theta$$, $$Y_2= P \sin \Theta+Q \cos \Theta$$.

## Code

I have implemented the simulation procedure in MATLAB, R and Python, which, as usual, are all very similar. I haven’t put my code here, because the software behind my website keeps mangling it.  As always, I have uploaded my code to a repository; for this post, it’s in this directory.

I have written the code in R, but I wouldn’t use it in general. That’s because if you’re using R, then, as I have said before, I strongly recommend using the powerful spatial statistics library spatstat. For a simulating Poisson line process, there’s the function rpoisline.

The chief author of spatstat, Adrian Baddeley, has written various lectures and books on the related topics of point processes, spatial statistics, and geometric probability. In this post, he answered why the angular coordinates have to be chosen uniformly.

### Results

MATLAB

R

Python

To read about the Poisson line process, it’s best to start with the Bertrand problem, which is covered in many works on geometric probability and related fields. A good and recent introduction is given by Calka in (Section 1.3) of the lectures titled Stochastic Geometry: Modern Research Frontiers, which were edited by Coupier and published by Springer.  Also see, for example, problem 1.2 in Geometrical Probability by Kendall and Moran or page 44 in Geometric Probability by Solomon.

For the Poisson line process, I highly recommend Section 7.2 in the monograph Poisson processes by Kingman. Also see Example 8.2 in the standard textbook Stochastic Geometry and its Applications by Chiu, Stoyan, Kendall and Mecke. The Poisson line process book also appears in Exercise 5.2 in the book Stochastic Simulation – Algorithms and Analysis by  Asmussen and Glynn.

For online resources, this set of lectures by Christian Lantuéjoul
covers the Poisson line process. Wilfrid Kendall covers the Poisson line process in this talk in relation to so-called Poisson cities.

## Testing the Julia language with point process simulations

I started writing these posts (or blog entries) about a year ago. In my first post I remarked how I wanted to learn to write stochastic simulations in a new language. Well, I found one. It’s called Julia. Here’s my code. And here are my thoughts.

### Overview

For scientific programming, the Julia language has arisen as a new contender. Originally started in 2012, its founders and developers have (very) high aspirations, wanting the language to be powerful and accessible, while still having run speeds comparable to C. There’s been excitement about it, and even a Nobel Laureate in economics, Thomas Sargent, has endorsed it. He co-founded the QuantEcon project, whose website has this handy guide or cheat sheet for commands between MATLAB, Python and Julia.

That guide suggests that Julia’s main syntax inspiration comes from MATLAB. But perhaps its closest (and greatest) competitor in scientific programming languages is Python, which has become a standard language used in scientific programming, particularly in machine learning. Another competitor is the statistics language R, which is popular for data science. But R is not renown for its speed.

As an aside, machine learning is closely related to what many call data science. I consider the two disciplines as largely overlapping with statistics, where their respective emphases are on theory and practice. In these fields, often the languages Python and R are used. There are various websites discussing which language is better, such as this one, which in turn is based on this one. In general, it appears that computer scientists and statisticians respectively prefer using Python and R.

Returning to the Julia language, given its young age, the language is still very much evolving, but I managed to find suitable Julia functions for stochastic simulations. I thought I would try it out by simulating some point processes, which I have done several times before. I successfully ran all my code with Julia Version 1.0.3.

In short, I managed to replicate in (or even translate to) Julia the code that I presented in the following posts:

Simulating a homogeneous Poisson point process on a rectangle

Simulating a Poisson point process on a disk

Simulating a Poisson point process on a triangle

Simulating an inhomogeneous Poisson point process

Simulating a Matérn cluster point process

Simulating a Thomas cluster point process

The Julia code, like all the code I present here, can be found on my Github repository, which for this post is located here.

### Basics

##### Language type and syntax

The Wikipedia article on Julia says:

Julia is a high-level general-purpose dynamic programming language designed for high-performance numerical analysis and computational science.

Scientific programming languages like the popular three MATLAB, R and Python, are interpreted languages. But the people behind Julia say:

it is a flexible dynamic language, appropriate for scientific and numerical computing, with performance comparable to traditional statically-typed languages.

Because Julia’s compiler is different from the interpreters used for languages like Python or R, you may find that Julia’s performance is unintuitive at first.

I already remarked that Julia’s syntax is clearly inspired by MATLAB, as one can see in this guide for MATLAB, Python and Julia. But there are key differences. For example, to access an array entry in Julia, you use square brackets (like in most programming languages), whereas parentheses are used in MATLAB.

##### Packages

Julia requires you to install certain packages or libraries, like most languages. For random simulations and plots, you have to install the Julia packages Distributions and Plots, which is done by running the code.

Pkg.add("Distributions");
Pkg.add("Plots");

After doing that, it’s best to restart Julia. These packages are loaded with the using command:

Using Distributions;
Using Plots;

Also, the first time it takes a while to run any code using those newly installed packages.

I should stress that there are different plotting libraries, but Plots, which contains many plotting libraries, is the only one I could get working. Another is PlotPy, which uses the Python library. As a beginner, it seems to me that the Julia community has not focused too much on developing new plotting functions, and has instead leveraged pre-existing libraries.

For standard scientific and statistical programming, you will usually also need the packages LinearAlgebra and Statistics.

##### Data types

Unlike MATLAB or R, Julia is a language that has different data types for numbers, such as integers and floating-point numbers (or floats). This puts Julia in agreement with the clear majority of languages, making it nothing new for most programmers. This is not a criticism of the language, but this can be troublesome if you’ve grown lazy after years of using MATLAB and R.

##### Simulating random variables

In MATLAB, R and Python, we just need to call a function for simulating uniform, Poisson, and other random variables. There’s usually a function for each type of random variable (or probability distribution).

Julia does simulation of random objects in a more, let’s say, object-oriented way (but I’m told, it’s not an object-oriented language). The probability distributions of random variables are objects, which are created and then sent to a general function for random generation. For example, here’s the code for simulating a Poisson variable with mean $$\mu=10$$.

mu=10;
distPoisson=Poisson(mu);
numbPoisson=rand(distPoisson);

Similarly, here’s how to simulate a normal variable with mean $$\mu=10$$ and standard deviation $$\sigma=1$$.

sigma=1;
distNormal=Normal(mu,sigma);
numbNormal=rand(distNormal);

Of course the last two lines can be collapsed into one.

mu=10;
sigma=1;
numbNormal=rand(Normal(mu,sigma));

But if you just want to create standard uniform variables on the interval (0,1), then the code is like that in MATLAB. For example, this code creates a $$4\times3$$ matrix (or array) $$X$$ whose entries are simulation outcomes of independent uniform random variables:

X=rand(4,3);

The resulting matrix $$X$$ is a Float 64 array.

##### Arrays

The indexing of arrays in Julia starts at one, like MATLAB and R. When you apply a function to an array, you generally need to use the dot notation. For example, if I try to run the code:

Y=sqrt(rand(10,1)); #This line will result in an error.

then on my machine (with Julia Version 1.0.3) I get the error:

ERROR: DimensionMismatch(“matrix is not square: dimensions are (10, 1)”)

But this code works:

Y=sqrt.(rand(10,1));

Also, adding scalars to arrays can catch you in Julia, as you also often need to use the dot notation. This code:

Y=sqrt.(rand(10,1));
Z=Y+1; #This line will result in an error.

gives the error:

ERROR: MethodError: no method matching +(::Array{Float64,2}, ::Int64)

This is fixed by adding a dot:

Y=sqrt.(rand(10,1));
Z=Y.+1; #This line will work.

Note the dot has to be on the left hand side. I ended up just using dot notation every time to be safe.

Other traps exist. For example, with indexing, you need to convert floats to integers if you want to use them as indices.

##### Repeating array elements

There used to be a Julia function called repmat, like the one in MATLAB , but it was merged with a function called repeat. I used such repeating operations to avoid explicit for-loops, which is generally advised in languages like MATLAB and R. For example, I used the repelem function in MATLAB to simulate Matérn and Thomas cluster point processes. To do this in Julia, I had to use this nested construction:

y=vcat(fill.(x, n)...);

This line means that the first value in $$x$$ is repeated $$n[1]$$ times, where $$n[1]$$ is the first entry of $$n$$ (as indexing in Julia starts at one), then the second value of $$x$$ is repeated $$n[2]$$ times, and so on. For example, the vectors $$x=[7,4,9]$$ and $$n=[2,1,3]$$, the answer is $$y=[7,7,4,9,9,9]$$.

To do this in Julia, the construction is not so bad, if you know how, but it’s not entirely obvious. In MATLAB I use this:

y=repelem(x,n);

Similarly in Python:

y=np.repeat(x,n);
##### Different versions of Julia

I found that certain code would work (or not work) and then later the same code would not work (or would work) on machines with different versions of Julia, demonstrating how the language is still being developed. More specifically, I ran code on Julia Version 1.0.3 (Date 2018-12-18) and Julia Version 0.6.4 (Date: 2018-07-09). (Note how there’s only a few months difference in the dates of the two versions.)

Consider the code with the errors (due to the lack of dot operator) in the previous section. The errors occurred on one machine with Julia Version 1.0.3, but the errors didn’t occur on another machine with the older Julia Version 0.6.4. For a specific example, the code:

Y=sqrt.(rand(10,1));
Z=Y+1; #This line will not result in an error on Version 0.6.4.

gives no error with Julia Version 0.6.4, while I have already discussed how it gives an error with Julia Version 1.0.3.

For another example, I copied from this MATLAB-Python-Julia guide the following command:

A = Diagonal([1,2,3]); #This line will (sometimes?) result in an error.

It runs on machine with Julia Version 0.6.4 with no problems. But running it on the machine with Julia Version 1.0.3 gives the error:

ERROR: UndefVarError: Diagonal not defined

That’s because I have not used the LinearAlgebra package. Fixing this, the following code:

using LinearAlgebra; #Package needed for Diagonal command.
A = Diagonal([1,2,3]); #This line should now work.

gives no error with Julia Version 1.0.3.

If you have the time and energy, you can search the internet and find online forums where the Julia developers have discussed why they have changed something, rendering certain code unworkable with the latest versions of Julia.

##### Optimization

It seems that performing optimization on functions is done with the Optim package.

Pkg.add("Optim");

But some functions need the Linesearches package, so it’s best to install that as well.

Pkg.add("Linesearches");

Despite those two optimization packages, I ended up using yet another package called BlackBoxOptim.

Pkg.add("BlackBoxOptim");

In this package, I used a function called bboptimize. This is the first optimziation function that I managed to get working. I do not know how it compares to the functions in the Optim and Linesearches packages.

In a previous post, I used optimization functions to simulate a inhomogeneous or nonhomogeneous Poisson point process on a rectangle. I’ve also written Julia code for this simulation, which is found below. I used bboptimize, but I had some problems when I initially set the search regions to integers, which the package did not like, as the values need to be floats. That’s why I multiple the rectangle dimensions by $$1.0$$ in the following code:

boundSearch=[(1.0xMin,1.0xMax), (1.0yMin, 1.0yMax)]; #bounds for search box
#WARNING: Values of boundSearch cannot be integers!
resultsOpt=bboptimize(fun_Neg;SearchRange = boundSearch);
lambdaNegMin=best_fitness(resultsOpt); #retrieve minimum value found by bboptimize

## Conclusion

In this brief experiment, I found the language Julia good for doing stochastic simulations, but too tricky for doing simple things like plotting. Also, depending on the version of Julia, sometimes my code would work and sometimes it wouldn’t. No doubt things will get better with time.

As I said, Julia is still very much an ongoing project. Here’s a couple of links that helped me learn the basics.

https://en.wikibooks.org/wiki/Introducing_Julia/Arrays_and_tuples

https://voxeu.org/content/which-numerical-computing-language-best-julia-matlab-python-or-r

Julia, Matlab, and C

https://modelingguru.nasa.gov/docs/DOC-2676

## Code

I’ve only posted here code for some of simulations, but the rest of the code is available on my GitHub repository located here. You can see how the code is comparable to that of MATLAB.

##### Poisson point process on a rectangle

I wrote about this point process here. The code is located here.

using Distributions #for random simulations
using Plots #for plotting

#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=rand(Poisson(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
plot1=scatter(xx,yy,xlabel ="x",ylabel ="y", leg=false);
display(plot1);
##### Inhomogeneous Poisson point process on a rectangle

I wrote about this point process here. The code is located here.

using Distributions #for random simulations
using Plots #for plotting
using BlackBoxOptim #for blackbox optimizing

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

s=0.5; #scale parameter

#Point process parameters
function fun_lambda(x,y)
100*exp.(-(x.^2+y.^2)/s^2); #intensity function
end

###START -- find maximum lambda -- START ###
#For an intensity function lambda, given by function fun_lambda,
#finds the maximum of lambda in a rectangular region given by
#[xMin,xMax,yMin,yMax].
#NOTE: Need xMin, xMax, yMin, yMax to be floats eg xMax=1. See boundSearch

function fun_Neg(x)
-fun_lambda(x[1],x[2]); #negative of lambda
end
xy0=[(xMin+xMax)/2.0,(yMin+yMax)/2.0];#initial value(ie centre)

#Find largest lambda value
boundSearch=[(1.0xMin,1.0xMax), (1.0yMin, 1.0yMax)];
#WARNING: Values of boundSearch cannot be integers!
resultsOpt=bboptimize(fun_Neg;SearchRange = boundSearch);
lambdaNegMin=best_fitness(resultsOpt); #retrieve minimum value found by bboptimize
lambdaMax=-lambdaNegMin;
###END -- find maximum lambda -- END ###

#define thinning probability function
function fun_p(x,y)
fun_lambda(x,y)/lambdaMax;
end

#Simulate a Poisson point process
numbPoints=rand(Poisson(areaTotal*lambdaMax)); #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

#calculate spatially-dependent thinning probabilities
p=fun_p(xx,yy);
#Generate Bernoulli variables (ie coin flips) for thinning
booleRetained=rand(numbPoints,1).<p; #points to be retained
xxRetained=xx[booleRetained]; yyRetained=yy[booleRetained];

#Plotting
plot1=scatter(xxRetained,yyRetained,xlabel ="x",ylabel ="y", leg=false);
display(plot1);
##### Thomas point process on a rectangle

I wrote about this point process here. The code is located here.

using Distributions #for random simulations
using Plots #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=10;#mean number of points in each cluster
sigma=0.05; #sigma for normal variables (ie random locations) of daughters

#Extended simulation windows parameters
rExt=7*sigma; #extension parameter
#for rExt, use factor of deviation sigma eg 6 or 7
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
numbPointsParent=rand(Poisson(areaTotalExt*lambdaParent)); #Poisson number of points

#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 Cartesian coordinates by
#simulating independent normal variables
xx0=rand(Normal(0,sigma),numbPoints);
yy0=rand(Normal(0,sigma),numbPoints);

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

#Shift centre of disk to (xx0,yy0)
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
plot1=scatter(xx,yy,xlabel ="x",ylabel ="y", leg=false);
display(plot1);

## Checking Poisson point process simulations

In previous posts I described how to simulate homogeneous Poisson point processes on a rectangledisk and triangle. Then I covered how to randomly thin a point process in a spatially dependent manner. Building off these posts, I  wrote in my last post how to simulate an inhomogeneous or nonhomogeneous Poisson point process. Now I’ll describe how to verify that the simulation code correctly simulates the desired Poisson point process.

I should add that although this post is focused on the Poisson point process, parts of the material hold for all point processes. Also, not surprisingly, some of the material and the code overlaps with that presented in the post on the inhomogeneous point process.

## Basics

Any Poisson point process is defined with a measure called the intensity measure or mean measure, which I’ll denote by $$\Lambda$$.  For practical purposes, I assume that the intensity measure $$\Lambda$$ has a derivative $$\lambda(x,y)$$, where I have used Cartesian coordinates.  The function $$\lambda(x,y)$$ is often called the intensity function or just intensity, which I  assume is bounded, so $$\lambda(x,y)<\infty$$ for all points in a simulation window $$W$$. Finally, I assume that the simulation window $$W$$ is a rectangle.

Several times before I have mentioned that simulating a Poisson point process requires simulating two random components: the number of points and the locations of points.  Working backwards, to check a Poisson simulation, we must run the Poisson simulation a large number of times (say $$10^3$$ or $$10^4$$), and collect the statistics on these two properties. We’ll start by examining the easiest of the two.

## Number of points

For any Poisson point process, the number of points is a Poisson random variable with a parameter (that is, a mean) $$\Lambda(W)$$, which under our previous assumptions is given by the surface integral

$$\Lambda(W)=\int_W \lambda(x,y)dxdy.$$

Presumably we can evaluate such an integral analytically or numerically in order to simulate the Poisson point process. To check that we correctly simulate the random the number of points, we just need to simulate the point process a large number  of times, and compare the statistics to those given by the analytic expressions.

##### Moments

The definition of the intensity measure of a point process is a measure that gives the average or expected number of points in some region. As the number of simulations increases, the (sample) average number of points will converge to the intensity measure $$\Lambda(W)$$.  I should stress that this is a test for the intensity measure, a type of first moment, which will work for the intensity measure of any point process.

For Poisson point processes, there is another moment test that can be done.  It can be shown mathematically that the variance of the number of points will also converge to the intensity measure $$\Lambda(W)$$, giving a second empirical test based on second moments. There is no pointp process theory here, as this moment result is simply due to the number of points being distributed according to a Poisson distribution. The second moment is very good for checking Poissonness, forming the basis for statistical tests. If this and the first moment test hold, then there’s a strong chance the number of points is a Poisson variable.

##### Empirical distribution

Beyond the first two moments, an even more thorough test is to estimate an empirical distribution for the number of points. That is, we perform a histogram on the number of points, and then we normalize it, so the total mass sums to one, producing an empirical probability distribution for the number of points.

The results should closely match the probability mass function of a Poisson distribution with parameter $$\Lambda(W)$$, namely

$$P(N=n)= \frac{[\Lambda(W)]^n}{n!} e^{-\Lambda(W)} ,$$

where $$N$$ is the random number of points in the window $$W$$, which has the area $$|W|$$. If the empirical distribution is close to results given by the above expression, then we can confidently say that the number of points is a Poisson random variable with the right parameter or mean.

## Locations of points

To check the positioning of points, we can empirically estimate the intensity function $$\lambda(x,y)$$. A simple way to do this is to perform a two-dimensional histogram on the point locations. This is very similar to the one-dimensional histogram, which I suggested to do for testing the number of points. It just involves counting the number of points that fall into two-dimensional non-overlapping subsets called bins.

To estimate the intensity function, each bin count needs to be rescaled by the area of the bin, giving a density value for each bin. This empirical estimate of the intensity function should resemble the true intensity function $$\lambda(x,y)$$. For a visual comparison, we can use a surface plot to illustrate the two sets of results.

This procedure will work for estimating the intensity function of any point process, not just a Poisson one.

In spatial statistics, there are more advanced statistical tests for testing how Poisson a point pattern is.  But these tests are arguably too complicated for checking a simple point process that is rather easy to simulate. Furthermore, researchers usually apply these tests to a small number of point patterns. In this setting, it is not possible to accurately obtain empirical distributions without further assumptions. But with simulations, we can generate many simulations and obtain good empirical distributions. In short, I would never use such tests for just checking that I have properly coded a Poisson simulation.

## Results

I produced the results with ten thousand simulations, which gave good results and took only a few seconds to complete on a standard desktop computer. Clearly increasing the number of simulations increases the accuracy of the statistics, but it also increases the computation time.

For the results, I used the intensity function

$$\lambda(x,y)=\lambda_1(x,y)+\lambda_2(x,y),$$

where

$$\lambda_1(x,y)=80e^{-((x+0.5)^2+(y+0.5)^2)/s^2},$$

$$\lambda_2(x,y)=100e^{-((x-0.5)^2+(y-0.5)^2)/s^2},$$

and $$s>0$$ is a scale parameter. We can see that this function has two maxima or peaks at $$(-0.5,-0.5)$$ and $$(0.5,0.5)$$.

##### Python

I have not covered much new theoretical stuff in this post, so looking at the references in previous posts, such as this one, should help.

For two-dimensional histograms, I recommend going to the respective MATLAB and Python function websites.  Here’s an example of a two-dimensional histogram implemented in Python.

Probably the most difficult part for me was performing the plotting in Python.  I recommend these links:

https://plot.ly/python/3d-surface-plots/

https://chrisalbon.com/python/basics/set_the_color_of_a_matplotlib/

https://matplotlib.org/examples/color/colormaps_reference.html

## Code

All code from my posts, as always, can be found on the my GitHub repository. The code for this post is located here.

The estimating the statistical moments is standard. Performing the histograms is also routine, but when normalizing, you have to choose the option that returns the empirical estimate of the probability density function (pdf).

Fortunately, scientific programming languages usually have functions for performing two-dimensional histograms. What is  a bit tricky is how to normalize or rescale the bin counts. The histogram functions can,  for example, divide by the number of simulations, the area of each bin, or both. In the end, I chose the pdf option  in both MATLAB and Python to give an empirical estimate of the probability density function, and then multiplied it by the average number of points, which was calculated in the previous check. (Although,  I could have done this in a single step in MATLAB, but not in Python, so I chose to do it in a couple of steps in both languages so the code matches more closely.)

##### MATLAB

I used the surf function to plot the intensity function and its estimate; see below for details on the histograms.

close all;

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

numbSim=10^4; %number of simulations

s=0.5; %scale parameter

%Point process parameters
fun_lambda=@(x,y)(100*exp(-((x).^2+(y).^2)/s^2));%intensity function

%%%START -- find maximum lambda -- START %%%
%For an intensity function lambda, given by function fun_lambda,
%finds the maximum of lambda in a rectangular region given by
%[xMin,xMax,yMin,yMax].
funNeg=@(x)(-fun_lambda(x(1),x(2))); %negative of lambda
%initial value(ie centre)
xy0=[(xMin+xMax)/2,(yMin+yMax)/2];%initial value(ie centre)
%Set up optimization step
options=optimoptions('fmincon','Display','off');
%Find largest lambda value
[~,lambdaNegMin]=fmincon(funNeg,xy0,[],[],[],[],...
[xMin,yMin],[xMax,yMax],'',options);
lambdaMax=-lambdaNegMin;
%%%END -- find maximum lambda -- END%%%

%define thinning probability function
fun_p=@(x,y)(fun_lambda(x,y)/lambdaMax);

%for collecting statistics -- set numbSim=1 for one simulation
numbPointsRetained=zeros(numbSim,1); %vector to record number of points
for ii=1:numbSim
%Simulate Poisson point process
numbPoints=poissrnd(areaTotal*lambdaMax);%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

%calculate spatially-dependent thinning probabilities
p=fun_p(xx,yy);

%Generate Bernoulli variables (ie coin flips) for thinning
booleRetained=rand(numbPoints,1)<p; %points to be retained %x/y locations of retained points xxRetained=xx(booleRetained); yyRetained=yy(booleRetained); %collect number of points simulated numbPointsRetained(ii)=length(xxRetained); end %Plotting plot(xxRetained,yyRetained,'bo'); %plot retained points xlabel('x');ylabel('y'); %run empirical test on number of points generated if numbSim&gt;=10
%total mean measure (average number of points)
LambdaNumerical=integral2(fun_lambda,xMin,xMax,yMin,yMax)
%Test: as numbSim increases, numbPointsMean converges to LambdaNumerical
numbPointsMean=mean(numbPointsRetained)
%Test: as numbSim increases, numbPointsVar converges to LambdaNumerical
numbPointsVar=var(numbPointsRetained)

end


For the histogram section, I used the histcounts and histcounts2 functions respectively to estimate the distribution of the number of points and the intensity function. I used the pdf option.

Number of points

histcounts(numbPointsRetained,binEdges,'Normalization','pdf');


Locations of points

histcounts2(xxVectorAll,yyVectorAll,numbBins,'Normalization','pdf');

##### Python

I used the Matplotlib library to plot the intensity function and its estimate; see below for details on the histograms.

import numpy as np; #NumPy package for arrays, random number generation, etc
import matplotlib.pyplot as plt #for plotting
from matplotlib import cm #for heatmap plotting
from mpl_toolkits import mplot3d #for 3-D plots
from scipy.optimize import minimize #for optimizing
from scipy import integrate #for integrating
from scipy.stats import poisson #for the Poisson probability mass function

plt.close("all"); #close all plots

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

numbSim=10**4; #number of simulations
numbBins=30; #number of bins for histogram

#Point process parameters
s=0.5; #scale parameter

def fun_lambda(x,y):
#intensity function
lambdaValue=80*np.exp(-((x+0.5)**2+(y+0.5)**2)/s**2)+100*np.exp(-((x-0.5)**2+(y-0.5)**2)/s**2);
return lambdaValue;

###START -- find maximum lambda -- START ###
#For an intensity function lambda, given by function fun_lambda,
#finds the maximum of lambda in a rectangular region given by
#[xMin,xMax,yMin,yMax].
def fun_Neg(x):
return -fun_lambda(x[0],x[1]); #negative of lambda

xy0=[(xMin+xMax)/2,(yMin+yMax)/2];#initial value(ie centre)
#Find largest lambda value
resultsOpt=minimize(fun_Neg,xy0,bounds=((xMin, xMax), (yMin, yMax)));
lambdaNegMin=resultsOpt.fun; #retrieve minimum value found by minimize
lambdaMax=-lambdaNegMin;
###END -- find maximum lambda -- END ###

#define thinning probability function
def fun_p(x,y):
return fun_lambda(x,y)/lambdaMax;

#for collecting statistics -- set numbSim=1 for one simulation
numbPointsRetained=np.zeros(numbSim); #vector to record number of points
xxAll=[]; yyAll=[];

### START -- Simulation section -- START ###
for ii in range(numbSim):
#Simulate a Poisson point process
numbPoints = np.random.poisson(lambdaMax*areaTotal);#Poisson number of points
xx = xDelta*np.random.uniform(0,1,numbPoints)+xMin;#x coordinates of Poisson points
yy = yDelta*np.random.uniform(0,1,numbPoints)+yMin;#y coordinates of Poisson points

#calculate spatially-dependent thinning probabilities
p=fun_p(xx,yy);

#Generate Bernoulli variables (ie coin flips) for thinning
booleRetained=np.random.uniform(0,1,numbPoints)<p; #points to be thinned

#x/y locations of retained points
xxRetained=xx[booleRetained]; yyRetained=yy[booleRetained];
numbPointsRetained[ii]=xxRetained.size;
xxAll.extend(xxRetained); yyAll.extend(yyRetained);
### END -- Simulation section -- END ###

#Plotting a simulation
fig1 = plt.figure();
plt.scatter(xxRetained,yyRetained, edgecolor='b', facecolor='none');
plt.xlabel("x"); plt.ylabel("y");
plt.title('A single realization of a Poisson point process');
plt.show();

#run empirical test on number of points generated
###START -- Checking number of points -- START###
#total mean measure (average number of points)
LambdaNumerical=integrate.dblquad(fun_lambda,xMin,xMax,lambda x: yMin,lambda y: yMax)[0];
#Test: as numbSim increases, numbPointsMean converges to LambdaNumerical
numbPointsMean=np.mean(numbPointsRetained);
#Test: as numbSim increases, numbPointsVar converges to LambdaNumerical
numbPointsVar=np.var(numbPointsRetained);
binEdges=np.arange(numbPointsRetained.min(),(numbPointsRetained.max()+1))-0.5;
pdfEmp, binEdges=np.histogram(numbPointsRetained, bins=binEdges,density=True);

nValues=np.arange(numbPointsRetained.min(),numbPointsRetained.max());
#analytic solution of probability density
pdfExact=(poisson.pmf(nValues,LambdaNumerical));

#Plotting
fig2 = plt.figure();
plt.scatter(nValues,pdfExact, color='b', marker='s',facecolor='none',label='Exact');
plt.scatter(nValues,pdfEmp, color='r', marker='+',label='Empirical');
plt.xlabel("n"); plt.ylabel("P(N=n)");
plt.title('Distribution of the number of points');
plt.legend();
plt.show();
###END -- Checking number of points -- END###

###START -- Checking locations -- START###
#2-D Histogram section
p_Estimate, xxEdges, yyEdges = np.histogram2d(xxAll, yyAll,bins=numbBins,density=True);
lambda_Estimate=p_Estimate*numbPointsMean;

xxValues=(xxEdges[1:]+xxEdges[0:xxEdges.size-1])/2;
yyValues=(yyEdges[1:]+yyEdges[0:yyEdges.size-1])/2;
X, Y = np.meshgrid(xxValues,yyValues) #create x/y matrices for plotting

#analytic solution of probability density
lambda_Exact=fun_lambda(X,Y);

#Plot empirical estimate
fig3 = plt.figure();
plt.rc('text', usetex=True);
plt.rc('font', family='serif');
ax=plt.subplot(211,projection='3d');
surf = ax.plot_surface(X, Y,lambda_Estimate,cmap=plt.cm.plasma);
plt.xlabel("x"); plt.ylabel("y");
plt.title('Estimate of $\lambda(x)$');
plt.locator_params(axis='x', nbins=5);
plt.locator_params(axis='y', nbins=5);
plt.locator_params(axis='z', nbins=3);
#Plot exact expression
ax=plt.subplot(212,projection='3d');
surf = ax.plot_surface(X, Y,lambda_Exact,cmap=plt.cm.plasma);
plt.xlabel("x"); plt.ylabel("y");
plt.title('True $\lambda(x)$');
plt.locator_params(axis='x', nbins=5);
plt.locator_params(axis='y', nbins=5);
plt.locator_params(axis='z', nbins=3);
###END -- Checking locations -- END###


For the histogram section, I used the histogram and histogram2d functions respectively to estimate the distribution of the number of points and the intensity function. I used the pdf option. (The SciPy website recommends not using the the normed option.)

Number of points

np.histogram(numbPointsRetained, bins=binEdges,density=True);


Locations of points

 np.histogram2d(xxArrayAll, yyArrayAll,bins=numbBins,density=True);