Implementation — reading a minimal Poisson sampler

Transcribe the formula directly into a function

This course's JavaScript is a minimal implementation that moves the math into functions without hiding anything. No external libraries; inputs that violate the preconditions raise an exception instead of being silently accepted.

The minimal JavaScript implementation

Here is the minimal implementation actually used by the Poisson calculations and the simulator.

function assertFiniteNonNegative(name, value) {
  if (!Number.isFinite(value) || value < 0) {
    throw new RangeError(`${name} must be a finite non-negative number`);
  }
}

function assertNonNegativeInteger(name, value) {
  if (!Number.isInteger(value) || value < 0) {
    throw new RangeError(`${name} must be a non-negative integer`);
  }
}

function factorialInt(n) {
  assertNonNegativeInteger("n", n);
  let acc = 1;
  for (let i = 2; i <= n; i += 1) {
    acc *= i;
  }
  return acc;
}

function poissonPmf(lambda, k) {
  assertFiniteNonNegative("lambda", lambda);
  assertNonNegativeInteger("k", k);
  return Math.exp(-lambda) * (lambda ** k) / factorialInt(k);
}

function poissonCdf(lambda, k) {
  assertFiniteNonNegative("lambda", lambda);
  assertNonNegativeInteger("k", k);
  let sum = 0;
  for (let i = 0; i <= k; i += 1) {
    sum += poissonPmf(lambda, i);
  }
  return sum;
}

function samplePoissonKnuth(lambda, rng = Math.random) {
  assertFiniteNonNegative("lambda", lambda);
  if (lambda === 0) return 0;

  const limit = Math.exp(-lambda);
  let product = 1;
  let count = 0;

  while (product > limit) {
    product *= rng();
    count += 1;
  }
  return count - 1;
}

poissonPmf gives a single bar, poissonCdf is the left-to-right cumulative sum, and samplePoissonKnuth draws one sample from a stream of uniform random numbers.

Code-to-formula mapping

As an introduction, being able to read poissonPmf and poissonCdf is plenty. The sampler is your entry point when you later move on to simulation or Monte Carlo work.

Why the Knuth sampler works

samplePoissonKnuth is a classic algorithm (Donald Knuth, 1969) for drawing a single sample from a Poisson distribution. The code looks like a black box at first, but the idea behind it is the well-known fact that "the inter-arrival times in a Poisson process follow an exponential distribution."

• For a Poisson process with rate λ (per unit interval), the gap between consecutive events follows Exponential(λ).
• The number of events X arriving in the interval [0, 1] follows Poisson(λ).
• Therefore "how many exponential gaps fit before we exhaust 1 unit of time" gives a Poisson sample.

Each exponential variate E_i can be produced from a uniform U_i ∈ (0, 1] as E_i = -ln(U_i) / λ. A Poisson(λ) sample is "the largest n such that E_1 + E_2 + ... + E_n ≤ 1." Multiplying both sides by -λ and exponentiating turns the sum into a product:

That equivalence is the heart of the algorithm. The code fixes limit = e^(−λ) up front, multiplies product by uniforms, and counts how many multiplications it takes for product to fall below limit. Concretely:

• product: the running product of uniforms drawn so far, U_1 × U_2 × ... × U_count.
• while (product > limit): keep looping while the product is still above e^(−λ), i.e., while the cumulative exponential sum has not yet exceeded 1.
• count - 1: the loop body increments count immediately after product *= rng(), so the iteration that finally breaks the condition has already incremented count. Subtracting 1 removes that overshoot, giving the desired count.

The early return 0 for lambda === 0 avoids the edge case where e^0 = 1 equals the initial product = 1, which would either stall the loop or behave unstably under floating-point rounding.

This algorithm is short and intuitive when λ is small, but it does O(λ) uniform draws per sample. Once λ gets large (say above 30) it becomes slow, and serious applications switch to algorithms designed for large λ such as PTRS or PA.

Implementation notes

This implementation is intentionally minimal for teaching. Four points to keep in mind if you push it toward production use:

• lambda ** k and ES2016: the exponentiation operator ** was added in ES2016. It is fine on modern Node.js and modern browsers, but if you need to support older environments, Math.pow(lambda, k) is broader. We chose ** here because it visually mirrors the math notation λ^k.
• Overflow in factorialInt: factorialInt is a plain loop. JavaScript Number has 53 bits of precision, so floating-point error creeps in around 22!, and 171! overflows to Infinity. For large λ paired with large k, you should switch to the recurrence P(X = k+1) = P(X = k) × λ / (k+1) or compute in log space. This implementation assumes the small-λ teaching regime.
• The Knuth sampler's count - 1: as explained above, by the time the loop exits, count already includes the iteration that broke the condition, so we return count - 1 to land on a valid Poisson(λ) sample. Starting count from 0 keeps the offset natural.
• Choice of exception: RangeError: RangeError is the right fit when the type is correct but the value is out of range, such as lambda < 0 or a non-integer k. TypeError is conventionally reserved for cases where the argument's type itself is wrong (for example, a string supplied where a number is required). For a teaching implementation, raising explicitly is much easier to trace than silent rounding or returning NaN.

Key takeaways from Chapter 6

• The Poisson pmf / cdf transcribe directly into JavaScript.
• The Knuth sampler rewrites "how many exponential gaps fit in interval 1" as a product of uniforms compared with e^(−λ).
• Raising a RangeError on invalid input makes the code easier to follow, both as a lesson and as an implementation.
• The simulator's numbers come from exactly this minimal implementation.