Merge pull request #2388 from gschwefer/pixel-likelihood-fix

Bug fix in full numeric pixel likelihood calculation

ctapipe/image/pixel_likelihood.py

@@ -29,6 +29,7 @@
 
 import numpy as np
 from scipy.integrate import quad
+from scipy.special import factorial
 from scipy.stats import poisson
 
 __all__ = [
@@ -71,11 +72,8 @@ def neg_log_likelihood_approx(image, prediction, spe_width, pedestal):
 
         - \\ln{P} = \\frac{\\ln{2 π} + \\ln{θ}}{2} + \\frac{(s - μ)^2}{2 θ}
 
-    and since we can remove constants and factors in the minimization:
-
-    .. math::
-
-        - \\ln{P} = \\ln{θ} + \\frac{(s - μ)^2}{θ}
+        We keep the constants in this because the actual value of the likelihood
+        can be used to calculate a goodness-of-fit value
 
 
     Parameters
@@ -95,13 +93,15 @@ def neg_log_likelihood_approx(image, prediction, spe_width, pedestal):
     """
     theta = pedestal**2 + prediction * (1 + spe_width**2)
 
-    neg_log_l = np.log(theta + EPSILON) + (image - prediction) ** 2 / theta
+    neg_log_l = 0.5 * (
+        np.log(2 * np.pi * theta + EPSILON) + (image - prediction) ** 2 / theta
+    )
 
     return np.sum(neg_log_l)
 
 
 def neg_log_likelihood_numeric(
-    image, prediction, spe_width, pedestal, confidence=(0.001, 0.999)
+    image, prediction, spe_width, pedestal, confidence=0.999
 ):
     """
     Calculate likelihood of prediction given the measured signal,
@@ -117,8 +117,9 @@ def neg_log_likelihood_numeric(
         Width of single p.e. peak (:math:`σ_γ`).
     pedestal: ndarray
         Width of pedestal (:math:`σ_p`).
-    confidence: tuple(float, float), 0 < x < 1
-        Confidence interval of poisson integration.
+    confidence: float, 0 < x < 1
+        Upper end of Poisson confidence interval of maximum prediction.
+        Determines upper end of poisson integration.
 
     Returns
     -------
@@ -131,16 +132,17 @@ def neg_log_likelihood_numeric(
 
     likelihood = epsilon
 
-    ns = np.arange(*poisson(np.max(prediction)).ppf(confidence))
+    n_signal = np.arange(poisson(np.max(prediction)).ppf(confidence) + 1)
 
-    ns = ns[ns >= 0]
+    n_signal = n_signal[n_signal >= 0]
 
-    for n in ns:
+    for n in n_signal:
         theta = pedestal**2 + n * spe_width**2
         _l = (
             prediction**n
             * np.exp(-prediction)
-            / theta
+            / np.sqrt(2 * np.pi * theta)
+            / factorial(n)
             * np.exp(-((image - n) ** 2) / (2 * theta))
         )
         likelihood += _l
@@ -222,7 +224,7 @@ def _integral_poisson_likelihood_full(image, prediction, spe_width, ped):
     image = np.asarray(image)
     prediction = np.asarray(prediction)
     like = neg_log_likelihood(image, prediction, spe_width, ped)
-    return like * np.exp(-0.5 * like)
+    return 2 * like * np.exp(-like)
 
 
 def mean_poisson_likelihood_full(prediction, spe_width, ped):

ctapipe/image/tests/test_pixel_likelihood.py

@@ -1,10 +1,12 @@
 import numpy as np
+
 from ctapipe.image import (
-    neg_log_likelihood,
-    neg_log_likelihood_approx,
-    mean_poisson_likelihood_gaussian,
     chi_squared,
     mean_poisson_likelihood_full,
+    mean_poisson_likelihood_gaussian,
+    neg_log_likelihood,
+    neg_log_likelihood_approx,
+    neg_log_likelihood_numeric,
 )
 
 
@@ -22,21 +24,43 @@ def test_chi_squared():
 
 
 def test_mean_poisson_likelihoood_gaussian():
-    prediction = np.array([1, 1, 1], dtype="float")
+    prediction = np.array([50, 50, 50], dtype="float")
     spe = 0.5
 
     small_mean_likelihood = mean_poisson_likelihood_gaussian(prediction, spe, 0)
     large_mean_likelihood = mean_poisson_likelihood_gaussian(prediction, spe, 1)
 
     assert small_mean_likelihood < large_mean_likelihood
 
+    # Test that the mean likelihood of abunch of samples drawn from the gaussian
+    # behind the aprroximate log likelihood is indeed the precalculated mean
+
+    rng = np.random.default_rng(123456)
+
+    ped = 1
+
+    mean_likelihood = mean_poisson_likelihood_gaussian(prediction[0], spe, ped)
+
+    distribution_width = np.sqrt(ped**2 + prediction[0] * (1 + spe**2))
+
+    normal_samples = rng.normal(
+        loc=prediction[0], scale=distribution_width, size=100000
+    )
+
+    rel_diff = (
+        2 * neg_log_likelihood_approx(normal_samples, prediction[0], spe, ped) / 100000
+        - mean_likelihood
+    ) / mean_likelihood
+
+    assert np.abs(rel_diff) < 5e-4
+
 
 def test_mean_poisson_likelihood_full():
-    prediction = np.array([30.0, 30.0])
+    prediction = np.array([10.0, 10.0])
 
     spe = np.array([0.5])
 
-    small_mean_likelihood = mean_poisson_likelihood_full(prediction, spe, [0])
+    small_mean_likelihood = mean_poisson_likelihood_full(prediction, spe, [0.1])
     large_mean_likelihood = mean_poisson_likelihood_full(prediction, spe, [1])
 
     assert small_mean_likelihood < large_mean_likelihood
@@ -56,27 +80,33 @@ def test_full_likelihood():
 
     full_like_small = neg_log_likelihood(image_small, expectation_small, spe, pedestal)
     exp_diff = full_like_small - np.sum(
-        np.asarray([2.75630505, 2.62168656, 3.39248449])
+        np.asarray([1.37815294, 1.31084662, 1.69627197])
     )
 
     # Check against known values
-    assert exp_diff / np.sum(full_like_small) < 1e-4
+    assert np.abs(exp_diff / np.sum(full_like_small)) < 1e-4
 
     image_large = np.array([40, 50, 60])
     expectation_large = np.array([50, 50, 50])
 
     full_like_large = neg_log_likelihood(image_large, expectation_large, spe, pedestal)
     # Check against known values
     exp_diff = full_like_large - np.sum(
-        np.asarray([7.45489137, 5.99305388, 7.66226007])
+        np.asarray([3.72744569, 2.99652694, 3.83113004])
     )
 
-    assert exp_diff / np.sum(full_like_large) < 1e-4
+    assert np.abs(exp_diff / np.sum(full_like_large)) < 3e-4
 
     gaus_like_large = neg_log_likelihood_approx(
         image_large, expectation_large, spe, pedestal
     )
 
+    numeric_like_large = neg_log_likelihood_numeric(
+        image_large, expectation_large, spe, pedestal
+    )
+
     # Check thats in large signal case the full expectation is equal to the
     # gaussian approximation (to 5%)
-    assert np.all(np.abs((full_like_large - gaus_like_large) / full_like_large) < 0.05)
+    assert np.all(
+        np.abs((numeric_like_large - gaus_like_large) / numeric_like_large) < 0.05
+    )

docs/changes/2388.bug.rst

@@ -0,0 +1 @@
+Fixed a bug in the calculation of the full numeric pixel likelihood and the corresponding tests.