In [1]:
import numpy as np

Here is some code to simulate a PPP of pigeons on a wire, and to count the number of pairs that are nearby to each other.

In [2]:
def draw_pigeons(L, rho=1):
    # draw number of pigeons
    N = np.random.poisson(L*rho)
    # draw positions
    X = np.random.uniform(low=0.0, high=L, size=N)
    return X

def count_pairs(X, d):
    # count number of pairs within distance d
    npairs = 0
    Y = sorted(X)
    for i, x in enumerate(Y):
        for j in range(i+1, len(X)):
            y = Y[j]
            if y - x <= d:
                npairs += 1
            else:
                break
    return npairs
In [3]:
X = np.array([0, 2, 2.1, 2.2, 3, 4, 1, 4.1, 10])
# X.sort()
assert(count_pairs(X, 0.5) == 4)
# count_pairs(X, 0.5)
In [4]:
L = 1e5  # length of wire in meters
rho = 1 # pigeon per meter
X = draw_pigeons(L, rho)    
len(X), count_pairs(X, 0.05), 0.1 * L * rho * rho / 2
Out[4]:
(100541, 4898, 5000.0)

Now let's find the mean of the number of nearby pairs by sampling from the distribution a lot of times:

In [6]:
nreps = 100
npairs = [count_pairs(draw_pigeons(L, rho), 0.05) for _ in range(nreps)]
mean_pairs = np.mean(npairs)
sd_pairs = np.std(npairs)
f"Average number of pairs = {mean_pairs}, with SD={sd_pairs}; expected value={0.1*L*rho*rho/2}"
Out[6]:
'Average number of pairs = 5004.24, with SD=74.42474319740714; expected value=5000.0'