Poisson Distribution

Compute the probability of exactly k events when they arrive independently at an average rate of lam per interval.

P(k) = lam^k * exp(-lam) / k!

Return 0.0 for negative k.

Written directly this overflows twice over for large k: lam^k goes to infinity while exp(-lam) goes to zero, and k! overflows a float around k = 171. Work in log space instead and exponentiate once at the end:

log P(k) = k * log(lam) - lam - lgamma(k + 1)

math.lgamma(k + 1) is log(k!) computed without ever forming k!. Trading a product of huge and tiny numbers for a sum of moderate ones is the standard fix, and it is the same reason models work with log probabilities throughout.

Handle lam = 0 separately, since log(0) is undefined: all the mass sits on k = 0.

Function signature

def poisson_pmf(k, lam):

Example

Input

poisson_pmf(k=2, lam=3.0)

Output

0.22404180765538786

Reasoning

3^2 * exp(-3) / 2! = 4.5 * 0.0498, which is about 0.2240.

Constraints and notes

• Return 0.0 for negative k.
• Work in log space with math.lgamma.
• When lam = 0, all the probability is on k = 0.