1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
use crate::{stirlerr::stirlerr, binodeviance::binodeviance, gammaln};
/// Poisson probability density function (PDF).
///
/// Returns the Poisson probability density function with parameter `lambda` at the value `x`.
///
/// # Mathematical Definition
/// <math display="block">
/// <mi>f</mi><mo>(</mo><mi>x</mi><mo>;</mo><mi>λ</mi><mo>)</mo>
/// <mo>=</mo>
/// <mfrac>
/// <mrow>
/// <msup><mi>λ</mi><mi>x</mi></msup>
/// <msup><mi>e</mi><mrow><mo>-</mo><mi>λ</mi></mrow></msup>
/// </mrow>
/// <mrow>
/// <mi>x</mi><mo>!</mo>
/// </mrow>
/// </mfrac>
/// </math>
/// for <math><mi>x</mi><mo>∈</mo><mo>{</mo><mn>0</mn><mo>,</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mo>...</mo><mo>}</mo></math>.
///
/// # Domain
/// - `x`: Non-negative integer.
/// - `lambda`: Non-negative real number.
/// - Returns `0.0` if `x` is not a non-negative integer.
/// - Returns `NaN` if `lambda < 0` or any input is `NaN`.
///
/// # Examples
/// ```
/// use abax::poisspdf;
///
/// // P(X=2 | lambda=3) = (3^2 * e^-3) / 2! = 9 * 0.04978706836 / 2 = 0.2240418
/// assert!((poisspdf(2.0, 3.0) - 0.2240418076553877).abs() < 1e-15);
///
/// // P(X=0 | lambda=0) = 1
/// assert_eq!(poisspdf(0.0, 0.0), 1.0);
/// ```
pub fn poisspdf(x: f64, lambda: f64) -> f64 {
// 1. Handle NaN propagating cases
if x.is_nan() || lambda.is_nan() {
return f64::NAN;
}
// 2. Invalid parameter space: lambda cannot be negative
if lambda < 0.0 {
return f64::NAN;
}
// 3. Exact evaluation edge case specified by MATLAB: x == 0 and lambda == 0
if x == 0.0 && lambda == 0.0 {
return 1.0;
}
// 4. The density function is strictly zero unless x is a non-negative integer
if x < 0.0 || x != x.round() {
return 0.0;
}
// 5. If lambda is exactly 0, any x > 0 has a probability mass of 0
if lambda == 0.0 {
return 0.0;
}
// --- Core Calculation Block ---
// At this stage, we are guaranteed: x >= 0, x == round(x), and lambda > 0.
// Constant: log(sqrt(2 * pi))
const LN_SR_2PI: f64 = 0.9189385332046727;
// Guard against underflow or division flags via realmin checks
let small_x = x <= lambda * f64::MIN_POSITIVE;
let large_x = lambda < x * f64::MIN_POSITIVE;
if small_x {
// If x is extremely small relative to lambda
(-lambda).exp()
} else if large_x {
// If lambda is extremely small relative to x
(-lambda + x * lambda.ln() - gammaln(x + 1.0)).exp()
} else {
// High-precision standard pathway using Loader's algorithm
(-LN_SR_2PI - 0.5 * x.ln() - stirlerr(x) - binodeviance(x, lambda)).exp()
}
}
#[cfg(test)]
mod tests {
use super::*;
const TOL: f64 = 1e-15;
#[test]
fn test_poisspdf_basic() {
// P(X=2 | lambda=3)
assert!((poisspdf(2.0, 3.0) - 0.2240418076553877).abs() < TOL);
// P(X=0 | lambda=1)
assert!((poisspdf(0.0, 1.0) - 0.36787944117144233).abs() < TOL);
// P(X=1 | lambda=1)
assert!((poisspdf(1.0, 1.0) - 0.36787944117144233).abs() < TOL);
// P(X=5 | lambda=2)
assert!((poisspdf(5.0, 2.0) - 0.0360894088630967).abs() < TOL);
}
#[test]
fn test_poisspdf_edge_cases() {
// x=0, lambda=0
assert_eq!(poisspdf(0.0, 0.0), 1.0);
// x > 0, lambda=0
assert_eq!(poisspdf(1.0, 0.0), 0.0);
// Non-integer x
assert_eq!(poisspdf(1.5, 1.0), 0.0);
// Negative x
assert_eq!(poisspdf(-1.0, 1.0), 0.0);
}
#[test]
fn test_poisspdf_invalid_inputs() {
assert!(poisspdf(f64::NAN, 1.0).is_nan());
assert!(poisspdf(1.0, f64::NAN).is_nan());
assert!(poisspdf(1.0, -1.0).is_nan());
}
#[test]
fn test_poisspdf_large_values() {
// Large x, lambda (triggers Loader's expansion)
assert!((poisspdf(100.0, 50.0) - 1.630319352147727e-10).abs() < 1e-14);
// Small lambda, large x (triggers direct log-probability formula)
assert!((poisspdf(10.0, 0.1) - 2.493489357462383e-17).abs() < 1e-32); // Adjusted tolerance for very small numbers
}
}