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.
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
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)
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
Now let's find the mean of the number of nearby pairs by sampling from the distribution a lot of times:
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}"