mir.stat.distribution.poisson source code

1 /++
2 This module contains algorithms for the $(LINK2 https://en.wikipedia.org/wiki/Poisson_distribution, Poisson Distribution).
3 
4 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
5 
6 Authors: John Michael Hall
7 
8 Copyright: 2022-3 Mir Stat Authors.
9 
10 +/
11 
12 module mir.stat.distribution.poisson;
13 
14 import mir.bignum.fp: Fp;
15 import mir.internal.utility: isFloatingPoint;
16 
17 /++
18 Algorithms used to calculate poisson dstribution.
19 
20 `PoissonAlgo.direct` can be more time-consuming for large values of the number
21 of events (`k`) or the rate of occurences (`lambda`). Additional algorithms are
22 provided to the user to choose the trade-off between running time and accuracy.
23 
24 See_also:
25     $(LINK2 https://en.wikipedia.org/wiki/Poisson_distribution, Poisson Distribution)
26 +/
27 enum PoissonAlgo {
28     /++
29     Direct
30     +/
31     direct,
32     /++
33     Gamma Incomplete function
34     +/
35     gamma,
36     /++
37     Approximates poisson distribution with normal distribution. Generally a better approximation when
38     `lambda > 1000`.
39     +/
40     approxNormal,
41     /++
42     Approximates poisson distribution with normal distribution (including continuity correction). More 
43     accurate than `PoissonAlgo.approxNormal`. Generally a better approximation when `lambda > 10`.
44     +/
45     approxNormalContinuityCorrection
46 }
47 
48 private
49 @safe pure nothrow @nogc
50 T poissonPMFImpl(T, PoissonAlgo poissonAlgo)(const size_t k, const T lambda)
51     if (isFloatingPoint!T && poissonAlgo == PoissonAlgo.direct)
52     in (lambda > 0, "lambda must be greater than or equal to 0")
53 {
54     import mir.math.common: exp, pow;
55     import mir.math.numeric: factorial;
56 
57     return exp(-lambda) * pow(lambda, k) / (cast(T) factorial(k));
58 }
59 
60 private
61 @safe pure nothrow @nogc
62 T poissonPMFImpl(T, PoissonAlgo poissonAlgo)(const size_t k, const T lambda)
63     if (isFloatingPoint!T && poissonAlgo == PoissonAlgo.gamma)
64     in (lambda > 0, "lambda must be greater than or equal to 0")
65 {
66     import std.mathspecial: gammaIncompleteCompl;
67 
68     return gammaIncompleteCompl(k + 1, lambda) - gammaIncompleteCompl(k, lambda);
69 }
70 
71 private
72 @safe pure nothrow @nogc
73 T poissonPMFImpl(T, PoissonAlgo poissonAlgo)(const size_t k, const T lambda)
74     if (isFloatingPoint!T && poissonAlgo == PoissonAlgo.approxNormal)
75     in (lambda > 0, "lambda must be greater than or equal to 0")
76 {
77     import mir.math.common: sqrt;
78     import mir.stat.distribution.normal: normalPDF;
79 
80     return normalPDF(k, lambda, sqrt(lambda));
81 }
82 
83 private
84 @safe pure nothrow @nogc
85 T poissonPMFImpl(T, PoissonAlgo poissonAlgo)(const size_t k, const T lambda)
86     if (isFloatingPoint!T && poissonAlgo == PoissonAlgo.approxNormalContinuityCorrection)
87     in (lambda > 0, "lambda must be greater than or equal to 0")
88 {
89     import mir.math.common: sqrt;
90     import mir.stat.distribution.normal: normalCDF;
91 
92     return normalCDF(cast(T) k + 0.5, lambda, sqrt(lambda)) - normalCDF(cast(T) k - 0.5, lambda, sqrt(lambda));
93 }
94 
95 /++
96 Computes the poisson probability mass function (PMF).
97 
98 Params:
99     poissonAlgo = algorithm for calculating PMF (default: PoissonAlgo.direct)
100 
101 See_also:
102     $(LINK2 https://en.wikipedia.org/wiki/Poisson_distribution, Poisson Distribution)
103 +/
104 @safe pure nothrow @nogc
105 template poissonPMF(PoissonAlgo poissonAlgo = PoissonAlgo.direct)
106 {
107     /++
108     Params:
109         k = value to evaluate PMF (e.g. number of events)
110         lambda = expected rate of occurence
111     +/
112     T poissonPMF(T)(const size_t k, const T lambda)
113         if (isFloatingPoint!T)
114         in (lambda > 0, "lambda must be greater than or equal to 0")
115     {
116         return poissonPMFImpl!(T, poissonAlgo)(k, lambda);
117     }
118 }
119 
120 /// ditto
121 @safe pure nothrow @nogc
122 template poissonPMF(string poissonAlgo)
123 {
124     mixin("alias poissonPMF = .poissonPMF!(PoissonAlgo." ~ poissonAlgo ~ ");");
125 }
126 
127 ///
128 version(mir_stat_test)
129 @safe pure nothrow @nogc
130 unittest {
131     import mir.math.common: approxEqual, exp;
132     
133     assert(3.poissonPMF(6.0).approxEqual(exp(-6.0) * 216 / 6));
134     // Can compute directly with differences of upper incomplete gamma function
135     assert(3.poissonPMF!"gamma"(6.0).approxEqual(poissonPMF(3, 6.0)));
136     // For large values of k or lambda, can approximate with normal distribution
137     assert(1_000_000.poissonPMF!"approxNormal"(1_000_000.0).approxEqual(poissonPMF!"gamma"(1_000_000, 1_000_000.0), 10e-3));
138     // Or closer with continuity correction
139     assert(1_000_000.poissonPMF!"approxNormalContinuityCorrection"(1_000_000.0).approxEqual(poissonPMF!"gamma"(1_000_000, 1_000_000.0), 10e-3));
140 }
141 
142 // test PoissonAlgo.direct
143 version(mir_stat_test)
144 @safe pure nothrow @nogc
145 unittest {
146     import mir.math.common: approxEqual, exp;
147 
148     assert(0.poissonPMF(5.0).approxEqual(exp(-5.0)));
149     assert(1.poissonPMF(5.0).approxEqual(exp(-5.0) * 5));
150     assert(2.poissonPMF(5.0).approxEqual(exp(-5.0) * 25 / 2));
151     assert(3.poissonPMF(5.0).approxEqual(exp(-5.0) * 125 / 6));
152     assert(4.poissonPMF(5.0).approxEqual(exp(-5.0) * 625 / 24));
153     assert(5.poissonPMF(5.0).approxEqual(exp(-5.0) * 3125 / 120));
154     assert(6.poissonPMF(5.0).approxEqual(exp(-5.0) * 15625 / 720));
155     assert(7.poissonPMF(5.0).approxEqual(exp(-5.0) * 78125 / 5040));
156     assert(8.poissonPMF(5.0).approxEqual(exp(-5.0) * 390625 / 40320));
157     assert(9.poissonPMF(5.0).approxEqual(exp(-5.0) * 1953125 / 362880));
158     assert(10.poissonPMF(5.0).approxEqual(exp(-5.0) * 9765625 / 3628800));
159 }
160 
161 // test PoissonAlgo.gamma
162 version(mir_stat_test)
163 @safe pure nothrow @nogc
164 unittest {
165     import mir.math.common: approxEqual;
166     for (size_t i; i < 20; i++) {
167         assert(i.poissonPMF!"gamma"(5.0).approxEqual(poissonPMF(i, 5.0)));
168     }
169 }
170 
171 // test PoissonAlgo.approxNormal / approxNormalContinuityCorrection
172 version(mir_stat_test)
173 @safe pure nothrow @nogc
174 unittest {
175     import mir.math.common: approxEqual, sqrt;
176     import mir.stat.distribution.normal: normalCDF, normalPDF;
177     for (size_t i; i < 20; i++) {
178         assert(i.poissonPMF!"approxNormal"(5.0).approxEqual(normalPDF(i, 5.0, sqrt(5.0))));
179         assert(i.poissonPMF!"approxNormalContinuityCorrection"(5.0).approxEqual(normalCDF(i + 0.5, 5.0, sqrt(5.0)) - normalCDF(i - 0.5, 5.0, sqrt(5.0))));
180     }
181 }
182 
183 /++
184 Computes the poisson probability mass function (PMF) directly with extended 
185 floating point types (e.g. `Fp!128`), which provides additional accuracy for
186 large values of `lambda` or `k`. 
187 
188 Params:
189     k = value to evaluate PMF (e.g. number of "heads")
190     lambda = expected rate of occurence
191 
192 See_also:
193     $(LINK2 https://en.wikipedia.org/wiki/Poisson_distribution, Poisson Distribution)
194 +/
195 @safe pure nothrow @nogc
196 T fp_poissonPMF(T)(const size_t k, const T lambda)
197     if (is(T == Fp!size, size_t size))
198     in (cast(double) lambda > 0, "lambda must be greater than or equal to 0")
199 {
200     import mir.math.common: exp;
201     import mir.math.internal.fp_powi: fp_powi;
202     import mir.math.numeric: factorial;
203 
204     return T(exp(-cast(double) lambda)) * fp_powi(lambda, k) / factorial(k);
205 }
206 
207 ///
208 version(mir_stat_test)
209 @safe pure nothrow @nogc
210 unittest {
211     import mir.bignum.fp: Fp;
212     import mir.conv: to;
213     import mir.math.common: approxEqual, exp;
214 
215     for (size_t i; i <= 10; i++) {
216         assert(i.fp_poissonPMF(Fp!128(5.0)).to!double.approxEqual(poissonPMF(i, 5.0)));
217     }
218 }
219 
220 private
221 @safe pure nothrow @nogc
222 T poissonCDFImpl(T, PoissonAlgo poissonAlgo)(const size_t k, const T lambda)
223     if (isFloatingPoint!T && poissonAlgo == PoissonAlgo.direct)
224     in (lambda > 0, "lambda must be greater than or equal to 0")
225 {
226     import mir.math.common: exp;
227     import mir.math.numeric: factorial;
228 
229     T output = 1;
230     if (k > 0) {
231         T multiplier = 1;
232         for (size_t i = 1; i < (k + 1); i++) {
233             multiplier *= (lambda / i);
234             output += multiplier;
235         }
236     }
237     return output * exp(-lambda);
238 }
239 
240 private
241 @safe pure nothrow @nogc
242 T poissonCDFImpl(T, PoissonAlgo poissonAlgo)(const size_t k, const T lambda)
243     if (isFloatingPoint!T && poissonAlgo == PoissonAlgo.gamma)
244     in (lambda > 0, "lambda must be greater than or equal to 0")
245 {
246     import std.mathspecial: gammaIncompleteCompl;
247     return cast(T) gammaIncompleteCompl(k + 1, lambda); 
248 }
249 
250 private
251 @safe pure nothrow @nogc
252 T poissonCDFImpl(T, PoissonAlgo poissonAlgo)(const size_t k, const T lambda)
253     if (isFloatingPoint!T && 
254         (poissonAlgo == PoissonAlgo.approxNormal || 
255          poissonAlgo == PoissonAlgo.approxNormalContinuityCorrection))
256     in (lambda > 0, "lambda must be greater than or equal to 0")
257 {
258     import mir.math.common: sqrt;
259     import mir.stat.distribution.normal: normalCDF;
260 
261     T l = k;
262     static if (poissonAlgo == PoissonAlgo.approxNormalContinuityCorrection) {
263         l += 0.5;
264     }
265     return normalCDF(l, lambda, sqrt(lambda));
266 }
267 
268 /++
269 Computes the poisson cumulative distrivution function (CDF).
270 
271 Params:
272     poissonAlgo = algorithm for calculating CDF (default: PoissonAlgo.direct)
273 
274 See_also:
275     $(LINK2 https://en.wikipedia.org/wiki/Poisson_distribution, Poisson Distribution)
276 +/
277 @safe pure nothrow @nogc
278 template poissonCDF(PoissonAlgo poissonAlgo = PoissonAlgo.direct)
279 {
280     /++
281     Params:
282         k = value to evaluate CDF (e.g. number of events)
283         lambda = expected rate of occurence
284     +/
285     T poissonCDF(T)(const size_t k, const T lambda)
286         if (isFloatingPoint!T)
287         in (lambda > 0, "lambda must be greater than or equal to 0")
288     {
289         return poissonCDFImpl!(T, poissonAlgo)(k, lambda);
290     }
291 }
292 
293 /// ditto
294 @safe pure nothrow @nogc
295 template poissonCDF(string poissonAlgo)
296 {
297     mixin("alias poissonCDF = .poissonCDF!(PoissonAlgo." ~ poissonAlgo ~ ");");
298 }
299 
300 ///
301 version(mir_stat_test)
302 @safe pure nothrow @nogc
303 unittest {
304     import mir.math.common: approxEqual;
305     
306     assert(3.poissonCDF(6.0).approxEqual(poissonPMF(0, 6.0) + poissonPMF(1, 6.0) + poissonPMF(2, 6.0) + poissonPMF(3, 6.0)));
307     // Can compute directly with upper incomplete gamma function
308     assert(3.poissonCDF!"gamma"(6.0).approxEqual(poissonCDF(3, 6.0)));
309     // For large values of k or lambda, can approximate with normal distribution
310     assert(1_000_000.poissonCDF!"approxNormal"(1_000_000.0).approxEqual(poissonCDF!"gamma"(1_000_000, 1_000_000.0), 10e-3));
311     // Or closer with continuity correction
312     assert(1_000_000.poissonCDF!"approxNormalContinuityCorrection"(1_000_000.0).approxEqual(poissonCDF!"gamma"(1_000_000, 1_000_000.0), 10e-3));
313 }
314 
315 // test PoissonAlgo.direct
316 version(mir_stat_test)
317 @safe pure nothrow @nogc
318 unittest {
319     import mir.math.common: approxEqual;
320     
321     static double sumOfPoissonPMFs(size_t k, double lambda) {
322         double output = 0;
323         for (size_t i; i < (k + 1); i++) {
324             output += poissonPMF(i, lambda);
325         }
326         return output;
327     }
328     
329     for (size_t i; i < 20; i++) {
330         assert(i.poissonCDF(5.0).approxEqual(sumOfPoissonPMFs(i, 5.0)));
331     }
332 }
333 
334 // test PoissonAlgo.gamma
335 version(mir_stat_test)
336 @safe pure nothrow @nogc
337 unittest {
338     import mir.math.common: approxEqual;
339     for (size_t i; i < 20; i++) {
340         assert(i.poissonCDF!"gamma"(5.0).approxEqual(poissonCDF(i, 5.0)));
341     }
342 }
343 
344 // test PoissonAlgo.approxNormal / approxNormalContinuityCorrection
345 version(mir_stat_test)
346 @safe pure nothrow @nogc
347 unittest {
348     import mir.math.common: approxEqual, sqrt;
349     import mir.stat.distribution.normal: normalCDF;
350     for (size_t i; i < 20; i++) {
351         assert(i.poissonCDF!"approxNormal"(5.0).approxEqual(normalCDF(i, 5.0, sqrt(5.0))));
352         assert(i.poissonCDF!"approxNormalContinuityCorrection"(5.0).approxEqual(normalCDF(i + 0.5, 5.0, sqrt(5.0))));
353     }
354 }
355 
356 private
357 @safe pure nothrow @nogc
358 T poissonCCDFImpl(T, PoissonAlgo poissonAlgo)(const size_t k, const T lambda)
359     if (isFloatingPoint!T && poissonAlgo == PoissonAlgo.direct)
360     in (lambda > 0, "lambda must be greater than or equal to 0")
361 {
362     return T(1) - poissonCDFImpl!(T, poissonAlgo)(k, lambda);
363 }
364 
365 private
366 @safe pure nothrow @nogc
367 T poissonCCDFImpl(T, PoissonAlgo poissonAlgo)(const size_t k, const T lambda)
368     if (isFloatingPoint!T && poissonAlgo == PoissonAlgo.gamma)
369     in (lambda > 0, "lambda must be greater than or equal to 0")
370 {
371     import std.mathspecial: gammaIncomplete;
372     return cast(T) gammaIncomplete(k + 1, lambda); 
373 }
374 
375 private
376 @safe pure nothrow @nogc
377 T poissonCCDFImpl(T, PoissonAlgo poissonAlgo)(const size_t k, const T lambda)
378     if (isFloatingPoint!T && 
379         (poissonAlgo == PoissonAlgo.approxNormal || 
380          poissonAlgo == PoissonAlgo.approxNormalContinuityCorrection))
381     in (lambda > 0, "lambda must be greater than or equal to 0")
382 {
383     import mir.math.common: sqrt;
384     import mir.stat.distribution.normal: normalCCDF;
385 
386     T l = k;
387     static if (poissonAlgo == PoissonAlgo.approxNormalContinuityCorrection) {
388         l += 0.5;
389     }
390     return normalCCDF(l, lambda, sqrt(lambda));
391 }
392 
393 /++
394 Computes the poisson complementary cumulative distrivution function (CCDF).
395 
396 Params:
397     poissonAlgo = algorithm for calculating CCDF (default: PoissonAlgo.direct)
398 
399 See_also:
400     $(LINK2 https://en.wikipedia.org/wiki/Poisson_distribution, Poisson Distribution)
401 +/
402 @safe pure nothrow @nogc
403 template poissonCCDF(PoissonAlgo poissonAlgo = PoissonAlgo.direct)
404 {
405     /++
406     Params:
407         k = value to evaluate CCDF (e.g. number of events)
408         lambda = expected rate of occurence
409     +/
410     T poissonCCDF(T)(const size_t k, const T lambda)
411         if (isFloatingPoint!T)
412         in (lambda > 0, "lambda must be greater than or equal to 0")
413     {
414         return poissonCCDFImpl!(T, poissonAlgo)(k, lambda);
415     }
416 }
417 
418 /// ditto
419 @safe pure nothrow @nogc
420 template poissonCCDF(string poissonAlgo)
421 {
422     mixin("alias poissonCCDF = .poissonCCDF!(PoissonAlgo." ~ poissonAlgo ~ ");");
423 }
424 
425 ///
426 version(mir_stat_test)
427 @safe pure nothrow @nogc
428 unittest {
429     import mir.math.common: approxEqual;
430 
431     assert(3.poissonCCDF(6.0).approxEqual(1.0 - (poissonPMF(0, 6.0) + poissonPMF(1, 6.0) + poissonPMF(2, 6.0) + poissonPMF(3, 6.0))));
432     // Can compute directly with upper incomplete gamma function
433     assert(3.poissonCCDF!"gamma"(6.0).approxEqual(poissonCCDF(3, 6.0)));
434     // For large values of k or lambda, can approximate with normal distribution
435     assert(1_000_000.poissonCCDF!"approxNormal"(1_000_000.0).approxEqual(poissonCCDF!"gamma"(1_000_000, 1_000_000.0), 10e-3));
436     // Or closer with continuity correction
437     assert(1_000_000.poissonCCDF!"approxNormalContinuityCorrection"(1_000_000.0).approxEqual(poissonCCDF!"gamma"(1_000_000, 1_000_000.0), 10e-3));
438 }
439 
440 // test PoissonAlgo.direct
441 version(mir_stat_test)
442 @safe pure nothrow @nogc
443 unittest {
444     import mir.math.common: approxEqual;
445 
446     for (size_t i; i < 20; i++) {
447         assert(i.poissonCCDF(5.0).approxEqual(1.0 - poissonCDF(i, 5.0)));
448     }
449 }
450 
451 // test PoissonAlgo.gamma
452 version(mir_stat_test)
453 @safe pure nothrow @nogc
454 unittest {
455     import mir.math.common: approxEqual;
456     for (size_t i; i < 20; i++) {
457         assert(i.poissonCCDF!"gamma"(5.0).approxEqual(poissonCCDF(i, 5.0)));
458     }
459 }
460 
461 // test PoissonAlgo.approxNormal / approxNormalContinuityCorrection
462 version(mir_stat_test)
463 @safe pure nothrow @nogc
464 unittest {
465     import mir.math.common: approxEqual, sqrt;
466     import mir.stat.distribution.normal: normalCCDF;
467     for (size_t i; i < 20; i++) {
468         assert(i.poissonCCDF!"approxNormal"(5.0).approxEqual(normalCCDF(i, 5.0, sqrt(5.0))));
469         assert(i.poissonCCDF!"approxNormalContinuityCorrection"(5.0).approxEqual(normalCCDF(i + 0.5, 5.0, sqrt(5.0))));
470     }
471 }
472 
473 private
474 @safe pure nothrow @nogc
475 size_t poissonInvCDFImpl(T, PoissonAlgo poissonAlgo)(const T p, const T lambda)
476     if (isFloatingPoint!T && poissonAlgo == PoissonAlgo.direct)
477     in (p >= 0, "p must be greater than or equal to 0")
478     in (p < 1, "p must be less than 1")
479     in (lambda > 0, "lambda must be greater than or equal to 0")
480 {
481     if (p == 0) {
482         return 0;
483     }
484 
485     size_t guess = 0;
486     if (lambda > 16) {
487         guess = poissonInvCDFImpl!(T, PoissonAlgo.approxNormalContinuityCorrection)(p, lambda);
488     }
489     T cdfGuess = poissonCDF!(poissonAlgo)(guess, lambda);
490 
491     if (p <= cdfGuess) {
492         if (guess == 0) {
493             return guess;
494         }
495         for (size_t i = (guess - 1); guess >= 0; i--) {
496             cdfGuess -= poissonPMF!(poissonAlgo)(i + 1, lambda);
497             if (p > cdfGuess) {
498                 guess = i + 1;
499                 break;
500             }
501         }
502     } else {
503         while(p > cdfGuess) {
504             guess++;
505             cdfGuess += poissonPMF!(poissonAlgo)(guess, lambda);
506         }
507     }
508     return guess;
509 }
510 
511 private
512 @safe pure nothrow @nogc
513 T poissonInvCDFImpl(T, PoissonAlgo poissonAlgo)(const T p, const size_t k)
514     if (isFloatingPoint!T && poissonAlgo == PoissonAlgo.gamma)
515     in (p >= 0, "p must be greater than or equal to 0")
516     in (p < 1, "p must be less than 1")
517 {
518     import std.mathspecial: gammaIncompleteComplInverse;
519 
520     if (p == 0) {
521         return T.infinity;
522     }
523     return gammaIncompleteComplInverse(k + 1, p); 
524 }
525 
526 private
527 @safe pure nothrow @nogc
528 size_t poissonInvCDFImpl(T, PoissonAlgo poissonAlgo)(const T p, const T lambda)
529     if (isFloatingPoint!T && 
530         (poissonAlgo == PoissonAlgo.approxNormal || 
531          poissonAlgo == PoissonAlgo.approxNormalContinuityCorrection))
532     in (p >= 0, "p must be greater than or equal to 0")
533     in (p < 1, "p must be less than 1")
534     in (lambda > 0, "lambda must be greater than or equal to 0")
535 {
536     import mir.math.common: floor, sqrt;
537     import mir.stat.distribution.normal: normalCDF, normalInvCDF;
538 
539     T std = sqrt(lambda);
540 
541     // Handles case where p is small better than just using pLowerBound = 0
542     T pLowerBound = 0;
543     T lowerValue = 0;
544     static if (poissonAlgo == PoissonAlgo.approxNormalContinuityCorrection)
545         lowerValue += 0.5;
546     pLowerBound = normalCDF(lowerValue, lambda, std);
547     if (p <= pLowerBound) { 
548         return 0;
549     }
550 
551     auto result = normalInvCDF(p, lambda, std);
552     static if (poissonAlgo == PoissonAlgo.approxNormalContinuityCorrection) {
553         result = result - 0.5;
554     }
555     return cast(size_t) floor(result);
556 }
557 
558 /++
559 Computes the poisson inverse cumulative distrivution function (InvCDF).
560 
561 For algorithms `PoissonAlgo.direct`, `PoissonAlgo.approxNormal`, and 
562 `PoissonAlgo.approxNormalContinuityCorrection`, the inverse CDF returns the 
563 number of events (`k`) given the probability (`p`) and rate of occurence
564 (`lambda`). For the `Poisson.gamma` algorithm, the inverse CDF returns the rate
565 of occurence (`lambda`) given the probability (`p`) and the number of events (`k`).
566 
567 For `PoissonAlgo.direct`, if the value of `lambda` is larger than 16, then an
568 initial guess is made based on `PoissonAlgo.approxNormalContinuityCorrection`.
569 
570 Params:
571     poissonAlgo = algorithm for calculating InvCDF (default: PoissonAlgo.direct)
572 
573 See_also:
574     $(LINK2 https://en.wikipedia.org/wiki/Poisson_distribution, Poisson Distribution)
575 +/
576 @safe pure nothrow @nogc
577 template poissonInvCDF(PoissonAlgo poissonAlgo = PoissonAlgo.direct)
578     if (poissonAlgo == PoissonAlgo.direct ||
579         poissonAlgo == PoissonAlgo.approxNormal ||
580         poissonAlgo == PoissonAlgo.approxNormalContinuityCorrection)
581 {
582     /++
583     Params:
584         p = value to evaluate InvCDF
585         lambda = expected rate of occurence
586     +/
587     size_t poissonInvCDF(T)(const T p, const T lambda)
588         if (isFloatingPoint!T)
589         in (p >= 0, "p must be greater than or equal to 0")
590         in (p <= 1, "p must be less than or equal to 1")
591         in (lambda > 0, "lambda must be greater than or equal to 0")
592     {
593         return poissonInvCDFImpl!(T, poissonAlgo)(p, lambda);
594     }
595 }
596 
597 /// ditto
598 @safe pure nothrow @nogc
599 template poissonInvCDF(PoissonAlgo poissonAlgo)
600     if (poissonAlgo == PoissonAlgo.gamma)
601 {
602     /++
603     Params:
604         p = value to evaluate InvCDF
605         k = number of events
606     +/
607     T poissonInvCDF(T)(const T p, const size_t k)
608         if (isFloatingPoint!T)
609         in (p >= 0, "p must be greater than or equal to 0")
610         in (p <= 1, "p must be less than or equal to 1")
611     {
612         return poissonInvCDFImpl!(T, poissonAlgo)(p, k);
613     }
614 }
615 
616 /// ditto
617 @safe pure nothrow @nogc
618 template poissonInvCDF(string poissonAlgo)
619 {
620     mixin("alias poissonInvCDF = .poissonInvCDF!(PoissonAlgo." ~ poissonAlgo ~ ");");
621 }
622 
623 ///
624 version(mir_stat_test)
625 @safe pure nothrow @nogc
626 unittest {
627     import mir.math.common: approxEqual;
628     
629     assert(0.15.poissonInvCDF(6.0) == 3);
630     // Passing `gamma` returns the rate of occurnece
631     assert(0.151204.poissonInvCDF!"gamma"(3).approxEqual(6));
632     // For large values of k or lambda, can approximate with normal distribution
633     assert(0.5.poissonInvCDF!"approxNormal"(1_000_000.0) == 1_000_000);
634     // Or closer with continuity correction
635     assert(0.5009974.poissonInvCDF!"approxNormalContinuityCorrection"(1_000_000.0) == 1_000_002);
636 }
637 
638 // test PoissonAlgo.direct
639 version(mir_stat_test)
640 @safe pure nothrow @nogc
641 unittest {
642     assert(0.poissonInvCDF(5.0) == 0);
643     for (double x = 0.05; x < 1; x = x + 0.05) {
644         size_t value = x.poissonInvCDF(5.0);
645         assert(value.poissonCDF(5.0) >= x);
646         assert((value - 1).poissonCDF(5.0) < x);
647     }
648 }
649 
650 // test PoissonAlgo.direct, large lambda branch, check small p
651 version(mir_stat_test)
652 @safe pure nothrow @nogc
653 unittest {
654     double x = 1.0e-9;
655     assert(x.poissonInvCDF(20.0) == 0);
656 }
657 
658 // test PoissonAlgo.direct, large lambda branch
659 version(mir_stat_test)
660 @safe pure nothrow @nogc
661 unittest {
662     assert(0.poissonInvCDF(25.0) == 0);
663     for (double x = 0.01; x < 1; x = x + 0.01) {
664         size_t value = x.poissonInvCDF(25.0);
665         assert(value.poissonCDF(25.0) >= x);
666         assert((value - 1).poissonCDF(25.0) < x);
667     }
668 }
669 
670 // test PoissonAlgo.gamma (note the difference in how it is tested)
671 version(mir_stat_test)
672 @safe pure nothrow @nogc
673 unittest {
674     import mir.math.common: approxEqual;
675 
676     assert(0.0.poissonInvCDF!"gamma"(5) == double.infinity);
677     for (double x = 0.05; x < 1; x = x + 0.05) {
678         for (size_t i; i < 10; i++) {
679             assert(poissonCDF!"gamma"(i, poissonInvCDF!"gamma"(x, i)).approxEqual(x));
680         }
681     }
682 }
683 
684 // test PoissonAlgo.approxNormal / approxNormalContinuityCorrection
685 version(mir_stat_test)
686 @safe pure nothrow @nogc
687 unittest {
688     import mir.math.common: floor, sqrt;
689     import mir.stat.distribution.normal: normalInvCDF;
690 
691     assert(0.poissonInvCDF!"approxNormal"(5.0) == 0);
692     assert(0.poissonInvCDF!"approxNormalContinuityCorrection"(5.0) == 0);
693     double checkValue;
694     for (double x = 0.05; x < 1; x = x + 0.05) {
695         checkValue = normalInvCDF(x, 5.0, sqrt(5.0));
696         assert(x.poissonInvCDF!"approxNormal"(5.0) == floor(checkValue));
697         assert(x.poissonInvCDF!"approxNormalContinuityCorrection"(5.0) == floor(checkValue - 0.5));
698     }
699 }
700 
701 /++
702 Computes the poisson log probability mass function (LPMF).
703 
704 Params:
705     k = value to evaluate LPMF (e.g. number of "heads")
706     lambda = expected rate of occurence
707 
708 See_also:
709     $(LINK2 https://en.wikipedia.org/wiki/Poisson_distribution, Poisson Distribution)
710 +/
711 @safe pure nothrow @nogc
712 T poissonLPMF(T)(const size_t k, const T lambda)
713     if (isFloatingPoint!T)
714     in (lambda > 0, "lambda must be greater than or equal to 0")
715 {
716     import mir.math.common: log;
717     import mir.math.internal.log_binomial: logFactorial;
718 
719     return k * log(lambda) - (logFactorial!T(k) + lambda);
720 }
721 
722 ///
723 version(mir_stat_test)
724 @safe pure nothrow @nogc
725 unittest {
726     import mir.math.common: approxEqual, exp;
727 
728     for (size_t i; i <= 10; i++) {
729         assert(i.poissonLPMF(5.0).exp.approxEqual(poissonPMF(i, 5.0)));
730     }
731 }
The OpenD Programming Language
