Solved – Poisson distribution sampling in 2D space with spatially varying intensity

poisson distributionpythonrandom variablespatial

I am sampling points from a Poisson process in a 2D space (Ny, Nx) grid with a constant intensity (lam) using the code given below.

How can I sample values from Poisson process in a space where the intensity is not constant, but varying spatially/regionally? I want something as shown in section 1.2.1 on page 7 (and section 7.2 on page 37) in this Tutorial using R. 

import numpy as np
from scipy.spatial import cKDTree as kdtree
import matplotlib.pyplot as plt

Nx, Ny, n_cells_reject_criteria = 100, 100, 3
valid = False

while not valid:
    rate_lambda = 0.02
    #===========generate random samples from homogeneous poisson process===========
    mean_poisson = rate_lambda*Nx*Ny
    n_events_pp = np.random.poisson(lam=mean_poisson)
    x_pp = np.round(np.random.uniform(low=0, high=Nx-1, size=n_events_pp)) # generate n uniformly distributed points
    y_pp = np.round(np.random.uniform(low=0, high=Ny-1, size=n_events_pp)) # generate n uniformly distributed points
    coords_random_ji = ([np.int(j) for j in y_pp], [np.int(i) for i in x_pp])

    #===========test there are no adjacent cells===========
    valid = len(kdtree(coords_random_ji).query_pairs(n_cells_reject_criteria)) == 0

#===========plot resuls===========
#------- create an empty mesh
grid = np.zeros((Ny, Nx), dtype=np.bool)

#------- superimpose the results from rejection sampling
grid[coords_random_ji] = True

#------- create empty figure
fig = plt.figure(figsize=(5, 5)) # in inches
#------- plot
plt.imshow(grid)

Example

As an example of what I mean by sampling with spatially varying intensity, the following figure shows the points (represented as black lines) sampled from a Poisson process superimposed on an image representing intensity colors (0 is the lowest intensity and 1 is the highest intensity), such that more points are sampled at locations with higher intensity and fewer points are sampled at locations with lower intensity.
Few things to note:

  1. Sampling from Poisson process is supposed to be allotted to the same
    grid/mesh dimension (=Nx*Ny) as that of intensity (=Nx*Ny) .
  2. Only 1 point per grid is allowed.
  3. The number of points must be less than or equal to the number of
    cells in the 2D grid (i.e. <= Nx*Ny).
  4. One cell in the grid cannot have more than 1 point.
  5. The location within the grid cell is not important, and by default,
    it is presumed the points are assigned to the center of the cell.
  6. I understand those are black lines and not points, but the intent is
    similar, i.e. the points are densely sampled at regions of higher
    intensity.

enter image description here

Best Answer

After several edits of the original question it is now clear what is wanted:

A grid of Bernoulli random variables (0-1 variables) with the probability of 1 proportional to a given intensity value in the grid. This can be seen as a discrete approximation to a Poisson point process in the limit of small pixel size.

I'm not a Python programmer so I can only give a description of the general algorithm/approach:

  • Let lam(i,j) be the intensity value in cell i,j.
  • Let c be a positive constant of proportionality (must be chosen such that c*lam(i,j)<=1 for all i,j).
  • For each i,j generate a Bernoulli random variable with success probability c*lam(i,j).

This will give you absence/presence indicators for each cell. The total expected number of presences is c times the sum lam(i,j) over all the i,j cells.

Related Question