A Poisson process from exponential waiting times

"""Building a Poisson process from exponential waiting times.

Arrivals of rate lam = 2 are generated by adding independent exponential
waiting times (seeded). We count the arrivals in [0, 3] over many runs and
compare the frequencies with the Poisson(6) probabilities, check that the
counts in [0, 1] and (1, 3] are uncorrelated (independent increments), test
the memoryless property of the waiting time, and verify the convolution
identity Po(s) * Po(t) = Po(s + t) term by term.
"""

import math
import random

def poisson_pmf(r, k):
    """e^(-r) r^k / k!"""
    if r < 0 or k < 0:
        raise ValueError("need r >= 0 and k >= 0")
    return math.exp(-r) * r ** k / math.factorial(k)

def arrival_times(lam, horizon, rng):
    """Arrival times in [0, horizon] from i.i.d. exponential gaps."""
    times, t = [], rng.expovariate(lam)
    while t <= horizon:
        times.append(t)
        t += rng.expovariate(lam)
    return times

def count_in(times, a, b):
    """Number of arrivals in the interval (a, b]."""
    return sum(1 for t in times if a < t <= b)

def convolution_gap(s, t, k):
    """|sum_j p_j(s) p_(k-j)(t) - p_k(s + t)|."""
    total = math.fsum(poisson_pmf(s, j) * poisson_pmf(t, k - j)
                      for j in range(k + 1))
    return abs(total - poisson_pmf(s + t, k))

def main():
    lam, runs = 2.0, 20000
    rng = random.Random(1909)
    firsts, seconds, totals, first_gap = [], [], [], []
    for _ in range(runs):
        times = arrival_times(lam, 3.0, rng)
        firsts.append(count_in(times, 0.0, 1.0))
        seconds.append(count_in(times, 1.0, 3.0))
        totals.append(len(times))
        first_gap.append(times[0] if times else math.inf)
    print(f"{runs} runs, rate 2, arrivals in [0, 3] versus Poisson(6):")
    print("   k   simulated   Poisson(6)")
    for k in [0, 2, 4, 6, 8, 10, 12]:
        freq = totals.count(k) / runs
        print(f"  {k:>2}   {freq:.4f}      {poisson_pmf(6.0, k):.4f}")
    mean = sum(totals) / runs
    var = sum((x - mean) ** 2 for x in totals) / runs
    print(f"mean {mean:.3f}, variance {var:.3f} (both 6 in theory)")
    m1, m2 = sum(firsts) / runs, sum(seconds) / runs
    cov = sum((a - m1) * (b - m2) for a, b in zip(firsts, seconds)) / runs
    print(f"Cov(N(0,1], N(1,3]) = {cov:.4f} (0 for independent increments)")
    t_wait = [x for x in first_gap if x > 0.5]
    cond = sum(1 for x in t_wait if x > 1.0) / len(t_wait)
    plain = sum(1 for x in first_gap if x > 0.5) / runs
    print(f"P(T > 1 | T > 0.5) = {cond:.4f}, P(T > 0.5) = {plain:.4f},"
          f" exact e^-1 = {math.exp(-1):.4f}")
    worst = max(convolution_gap(1.3, 2.4, k) for k in range(30))
    print(f"max |Po(1.3) * Po(2.4) - Po(3.7)| over k < 30: {worst:.1e}")
    print()
    print("Epistemic status: the frequencies, covariance and conditional")
    print("probability come from seeded simulations (statistical evidence,")
    print("with errors of order 1/sqrt(runs) = 0.007). The convolution check")
    print("is floating point; the identity is proved exactly in the lesson and")
    print("in Lean.")

if __name__ == "__main__":
    main()