feat: Allow zero rate Poisson

setup.py

@@ -6,7 +6,7 @@
         'tensorflow~=2.2,>=2.2.1,!=2.3.0',  # c.f. https://github.com/tensorflow/tensorflow/pull/40789
         'tensorflow-probability~=0.10,>=0.10.1',
     ],
-    'torch': ['torch~=1.8'],
+    'torch': ['torch~=1.10'],
     'jax': ['jax~=0.2.8', 'jaxlib~=0.1.58,!=0.1.68'],  # c.f. Issue 1501
     'xmlio': ['uproot>=4.1.1'],
     'minuit': ['iminuit>=2.4'],

src/pyhf/tensor/jax_backend.py

@@ -3,7 +3,7 @@
 config.update('jax_enable_x64', True)
 
 import jax.numpy as jnp
-from jax.scipy.special import gammaln
+from jax.scipy.special import gammaln, xlogy
 from jax.scipy import special
 from jax.scipy.stats import norm
 import numpy as np
@@ -368,14 +368,20 @@ def einsum(self, subscripts, *operands):
     def poisson_logpdf(self, n, lam):
         n = jnp.asarray(n)
         lam = jnp.asarray(lam)
-        return n * jnp.log(lam) - lam - gammaln(n + 1.0)
+        return xlogy(n, lam) - lam - gammaln(n + 1.0)
 
     def poisson(self, n, lam):
         r"""
         The continuous approximation, using :math:`n! = \Gamma\left(n+1\right)`,
         to the probability mass function of the Poisson distribution evaluated
         at :code:`n` given the parameter :code:`lam`.
 
+        .. note::
+
+            Though the p.m.f of the Poisson distribution is not defined for
+            :math:`\lambda = 0`, the limit :math:`\lambda \to 0` is still
+            defined and returns ``1.0`` in the case where ``n=0`` and ``lam=0``.
+
         Example:
 
             >>> import pyhf
@@ -398,7 +404,7 @@ def poisson(self, n, lam):
         """
         n = jnp.asarray(n)
         lam = jnp.asarray(lam)
-        return jnp.exp(n * jnp.log(lam) - lam - gammaln(n + 1.0))
+        return jnp.exp(xlogy(n, lam) - lam - gammaln(n + 1.0))
 
     def normal_logpdf(self, x, mu, sigma):
         # this is much faster than

src/pyhf/tensor/numpy_backend.py

@@ -1,7 +1,7 @@
 """NumPy Tensor Library Module."""
 import numpy as np
 import logging
-from scipy.special import gammaln
+from scipy.special import gammaln, xlogy
 from scipy import special
 from scipy.stats import norm, poisson
 
@@ -349,14 +349,20 @@ def einsum(self, subscripts, *operands):
         return np.einsum(subscripts, *operands)
 
     def poisson_logpdf(self, n, lam):
-        return n * np.log(lam) - lam - gammaln(n + 1.0)
+        return xlogy(n, lam) - lam - gammaln(n + 1.0)
 
     def poisson(self, n, lam):
         r"""
         The continuous approximation, using :math:`n! = \Gamma\left(n+1\right)`,
         to the probability mass function of the Poisson distribution evaluated
         at :code:`n` given the parameter :code:`lam`.
 
+        .. note::
+
+            Though the p.m.f of the Poisson distribution is not defined for
+            :math:`\lambda = 0`, the limit :math:`\lambda \to 0` is still
+            defined and returns ``1.0`` in the case where ``n=0`` and ``lam=0``.
+
         Example:
 
             >>> import pyhf
@@ -379,7 +385,8 @@ def poisson(self, n, lam):
         """
         n = np.asarray(n)
         lam = np.asarray(lam)
-        return np.exp(n * np.log(lam) - lam - gammaln(n + 1.0))
+
+        return np.exp(xlogy(n, lam) - lam - gammaln(n + 1.0))
 
     def normal_logpdf(self, x, mu, sigma):
         # this is much faster than

src/pyhf/tensor/pytorch_backend.py

@@ -367,6 +367,12 @@ def poisson(self, n, lam):
         to the probability mass function of the Poisson distribution evaluated
         at :code:`n` given the parameter :code:`lam`.
 
+        .. note::
+
+            Though the p.m.f of the Poisson distribution is not defined for
+            :math:`\lambda = 0`, the limit :math:`\lambda \to 0` is still
+            defined and returns ``1.0`` in the case where ``n=0`` and ``lam=0``.
+
         Example:
 
             >>> import pyhf

src/pyhf/tensor/tensorflow_backend.py

@@ -2,7 +2,6 @@
 import logging
 import tensorflow as tf
 import tensorflow_probability as tfp
-from numpy import nan
 
 log = logging.getLogger(__name__)
 
@@ -427,22 +426,20 @@ def poisson_logpdf(self, n, lam):
             TensorFlow Tensor: Value of the continuous approximation to log(Poisson(n|lam))
         """
         lam = self.astensor(lam)
-        # Guard against Poisson(n=0 | lam=0)
-        # c.f. https://github.com/scikit-hep/pyhf/issues/293
-        valid_obs_given_rate = tf.logical_or(
-            tf.math.not_equal(lam, n), tf.math.not_equal(n, 0)
-        )
-
-        return tf.where(
-            valid_obs_given_rate, tfp.distributions.Poisson(lam).log_prob(n), nan
-        )
+        return tfp.distributions.Poisson(lam).log_prob(n)
 
     def poisson(self, n, lam):
         r"""
         The continuous approximation, using :math:`n! = \Gamma\left(n+1\right)`,
         to the probability mass function of the Poisson distribution evaluated
         at :code:`n` given the parameter :code:`lam`.
 
+        .. note::
+
+            Though the p.m.f of the Poisson distribution is not defined for
+            :math:`\lambda = 0`, the limit :math:`\lambda \to 0` is still
+            defined and returns ``1.0`` in the case where ``n=0`` and ``lam=0``.
+
         Example:
             >>> import pyhf
             >>> pyhf.set_backend("tensorflow")
@@ -465,17 +462,7 @@ def poisson(self, n, lam):
             TensorFlow Tensor: Value of the continuous approximation to Poisson(n|lam)
         """
         lam = self.astensor(lam)
-        # Guard against Poisson(n=0 | lam=0)
-        # c.f. https://github.com/scikit-hep/pyhf/issues/293
-        valid_obs_given_rate = tf.logical_or(
-            tf.math.not_equal(lam, n), tf.math.not_equal(n, 0)
-        )
-
-        return tf.where(
-            valid_obs_given_rate,
-            tf.exp(tfp.distributions.Poisson(lam).log_prob(n)),
-            nan,
-        )
+        return tf.exp(tfp.distributions.Poisson(lam).log_prob(n))
 
     def normal_logpdf(self, x, mu, sigma):
         r"""

tests/test_tensor.py

@@ -287,17 +287,14 @@ def test_pdf_calculations(backend):
         ],
         nan_ok=True,
     )
-    # poisson(lambda=0) is not defined, should return NaN
+    # Allow poisson(lambda=0) under limit Poisson(lambda -> 0) = 1
     assert tb.tolist(
         tb.poisson(tb.astensor([0, 0, 1, 1]), tb.astensor([0, 1, 0, 1]))
-    ) == pytest.approx(
-        [np.nan, 0.3678794503211975, 0.0, 0.3678794503211975], nan_ok=True
-    )
+    ) == pytest.approx([1.0, 0.3678794503211975, 0.0, 0.3678794503211975])
     assert tb.tolist(
         tb.poisson_logpdf(tb.astensor([0, 0, 1, 1]), tb.astensor([0, 1, 0, 1]))
     ) == pytest.approx(
-        np.log([np.nan, 0.3678794503211975, 0.0, 0.3678794503211975]).tolist(),
-        nan_ok=True,
+        np.log([1.0, 0.3678794503211975, 0.0, 0.3678794503211975]).tolist()
     )
 
     # Ensure continuous approximation is valid
@@ -333,17 +330,14 @@ def test_pdf_calculations_pytorch(backend):
         ],
     )
 
-    # poisson(lambda=0) is not defined, should return NaN
+    # Allow poisson(lambda=0) under limit Poisson(lambda -> 0) = 1
     assert tb.tolist(
         tb.poisson(tb.astensor([0, 0, 1, 1]), tb.astensor([0, 1, 0, 1]))
-    ) == pytest.approx(
-        [np.nan, 0.3678794503211975, 0.0, 0.3678794503211975], nan_ok=True
-    )
+    ) == pytest.approx([1.0, 0.3678794503211975, 0.0, 0.3678794503211975])
     assert tb.tolist(
         tb.poisson_logpdf(tb.astensor([0, 0, 1, 1]), tb.astensor([0, 1, 0, 1]))
     ) == pytest.approx(
-        np.log([np.nan, 0.3678794503211975, 0.0, 0.3678794503211975]).tolist(),
-        nan_ok=True,
+        np.log([1.0, 0.3678794503211975, 0.0, 0.3678794503211975]).tolist()
     )
 
     # Ensure continuous approximation is valid