Poisson distribution sampling in 2D space with spatially varying intensity?

Question

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.

Answer 1

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.