import mir.math.common: approxEqual, exp;

assert(3.poissonPMF(6.0).approxEqual(exp(-6.0) * 216 / 6));
// Can compute directly with differences of upper incomplete gamma function
assert(3.poissonPMF!"gamma"(6.0).approxEqual(poissonPMF(3, 6.0)));
// For large values of k or lambda, can approximate with normal distribution
assert(1_000_000.poissonPMF!"approxNormal"(1_000_000.0).approxEqual(poissonPMF!"gamma"(1_000_000, 1_000_000.0), 10e-3));
// Or closer with continuity correction
assert(1_000_000.poissonPMF!"approxNormalContinuityCorrection"(1_000_000.0).approxEqual(poissonPMF!"gamma"(1_000_000, 1_000_000.0), 10e-3));