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


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

  • $\begingroup$ I think you might do better on a programming site? $\endgroup$
    – mdewey
    Oct 24, 2017 at 11:57
  • $\begingroup$ I did try first there...but, couldn't find any help there. $\endgroup$
    – user11
    Oct 24, 2017 at 12:05

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

  • $\begingroup$ This is a good approximation. A way to do it "exactly" is thinning a homogenous process. $\endgroup$ Oct 25, 2017 at 10:48
  • $\begingroup$ Depending on what you mean by an inhomogeneous Poisson process there is no approximation here: If you have an inhomogeneous Poisson point process in space with intensity function lambda(x,y) and you then define a grid of cells and count the number of points in each grid then these numbers are independent Poisson variables with the integral of lambda over the cell as the mean value. No approximation is done. Of course if you want to simulate the continuous point process you need to thin a homogeneous process and you have no need for the grid at all. $\endgroup$
    – Ege Rubak
    Oct 25, 2017 at 12:46
  • $\begingroup$ Yes, I was assuming that the underlying point process was needed but you're right that OP doesn't seem to want that. $\endgroup$ Oct 25, 2017 at 13:41
  • $\begingroup$ @EgeRubak: I may not be correct in what I am doing, but my objective is to randomly generate a location according to a Poisson process that is proportional to regional (spatial) intensity map. Can you suggest a solution? $\endgroup$
    – user11
    Oct 25, 2017 at 18:22
  • $\begingroup$ Are the random locations continuously distributed on the map (events possible anywhere in the region) or do you want to work on some grid? And what is the "intensity map"? A spatially continuous function or a fine grid/raster of intensity values over the map? Possibly edit your question with some figures of what you have and what you want. $\endgroup$
    – Ege Rubak
    Oct 26, 2017 at 11:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.