from scipy.special import i0
import numpy as np

sample_count = 100000
mean = 10

samples = np.random.poisson(mean, (sample_count, 2))
x = samples[:, 0]
y = samples[:, 1]
proportion = sum(x <= y) / sample_count
predicted = (1 + np.exp(-2*mean) * i0(2*mean)) / 2

print("mean: {}".format(10))
print("proportion of {} samples with x <= y: {}".format(sample_count, proportion))
print("predicted proportion: {}".format(predicted))