Poisson Point Process in Python 3 with numpy, without scipy

Question

I need to write a function in Python 3 which returns an array of positions (x,y) on a rectangular field (e.g. 100x100 points) that are scattered according to a homogenous spatial Poisson process.

So far I have found this resource with Python code, but unfortunately, I'm unable to find/install scipy for Python 3:

http://connor-johnson.com/2014/02/25/spatial-point-processes/

It has helped me understand what a Poisson point process actually is and how it works, though.

I have been playing around with numpy.random.poisson for a while now, but I am having a tough time interpreting what it returns.

http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.poisson.html

>>> import numpy as np
>>> np.random.poisson(1, (1, 5, 5))
array([[[0, 2, 0, 1, 0],
        [3, 2, 0, 2, 1],
        [0, 1, 3, 3, 2],
        [0, 1, 2, 0, 2],
        [1, 2, 1, 0, 3]]])

What I think that command does is creating one 5x5 field = (1, 5, 5) and scattering objects with a rate of lambda = 1 over that field. The numbers displayed in the resulting matrix are the probability of an object lying on that specific position.

How can I scatter, say, ten objects over that 5x5 field according to a homogenous spatial Poisson process? My first guess would be to iterate over the whole array and insert an object on every position with a "3", then one on every other position with a "2", and so on, but I'm unsure of the actual probability I should use to determine if an object should be inserted or not.

According to the following resource, I can simulate 10 objects being scattered over a field with a rate of 1 by simply multiplying the rate and the object count (10*1 = 10) and using that value as my lambda, i.e.

>>> np.random.poisson(10, (1, 5, 5))
array([[[12, 12, 10, 16, 16],
        [ 8,  6,  8, 12,  9],
        [12,  4, 10,  3,  8],
        [15, 10, 10, 15,  7],
        [ 8, 13, 12,  9,  7]]])

However, I don't see how that should make things easier. I only increase the rate at which objects appear by 10 that way.

Poisson point process in matlab

To sum it up, my primary question is: How can I use numpy.random.poisson(lam, size) to model a number n of objects being scattered over a 2-dimensional field dx*dy?

Answer 1

That's all correct. You definitely don't need SciPy, though when I first simulated a Poisson point process in Python I also used SciPy. I presented the original code with details in the simulation process in this post:

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

I just use NumPy in the more recent code:

import numpy as np; #NumPy package for arrays, random number generation, etc
import matplotlib.pyplot as plt #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
lambda0=100; #intensity (ie mean density) of the Poisson process

#Simulate a Poisson point process
numbPoints = np.random.poisson(lambda0*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

The code can also be found here:

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

I've also uploaded there more Python (and MATLAB and Julia) code for simulating several points processes, including Poisson point processes on various shapes and cluster point processes.

https://github.com/hpaulkeeler/posts

Answer 2

It seems I've looked at the problem in the wrong way. After more offline research I found out that it actually is sufficient to create a random Poisson value which represents the number of objects, for example n = np.random.poisson(100) and create the same amount of random values between 0 and 1

x = np.random.rand(n)
y = np.random.rand(n)

Now I just need to join the two arrays of x- and y-values to an array of (x,y) tuples. Those are the random positions I was looking for. I can multiply every x and y value by the side length of my field, e.g. 100, to scale the values to the 100x100 field I want to display.

I thought that the "randomness" of those positions should be determined by a random Poisson process, but it seems that just the number of positions needs to be determined by it, not the actual positional values.