From 7b15a908042f53c8f1fee9e511f0697736c58523 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Thu, 24 Aug 2023 15:53:36 +0200 Subject: [PATCH 01/12] Fixed pixel likelihood formulation --- ctapipe/image/pixel_likelihood.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py index c29db09ac94..8a8412d151d 100644 --- a/ctapipe/image/pixel_likelihood.py +++ b/ctapipe/image/pixel_likelihood.py @@ -30,6 +30,7 @@ import numpy as np from scipy.integrate import quad from scipy.stats import poisson +from scipy.special import factorial __all__ = [ "neg_log_likelihood_approx", @@ -95,13 +96,14 @@ 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 + # This is really 2 times the full log likelihood + neg_log_l = 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, @@ -131,7 +133,7 @@ def neg_log_likelihood_numeric( likelihood = epsilon - ns = np.arange(*poisson(np.max(prediction)).ppf(confidence)) + ns = np.arange(poisson(np.max(prediction)).ppf(confidence) + 1) ns = ns[ns >= 0] @@ -140,7 +142,8 @@ def neg_log_likelihood_numeric( _l = ( prediction**n * np.exp(-prediction) - / theta + / np.sqrt(2 * np.pi * theta) + / factorial(n) * np.exp(-((image - n) ** 2) / (2 * theta)) ) likelihood += _l From 6766c7946bb0f9a5d5cd5b73fd07bdb0dca81fb8 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Thu, 24 Aug 2023 16:07:13 +0200 Subject: [PATCH 02/12] Fixed tests, passing now --- ctapipe/image/pixel_likelihood.py | 5 +++-- ctapipe/image/tests/test_pixel_likelihood.py | 17 ++++++++++++----- 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py index 8a8412d151d..fec2247d876 100644 --- a/ctapipe/image/pixel_likelihood.py +++ b/ctapipe/image/pixel_likelihood.py @@ -96,8 +96,9 @@ def neg_log_likelihood_approx(image, prediction, spe_width, pedestal): """ theta = pedestal**2 + prediction * (1 + spe_width**2) - # This is really 2 times the full log likelihood - neg_log_l = np.log(2 * np.pi * 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) diff --git a/ctapipe/image/tests/test_pixel_likelihood.py b/ctapipe/image/tests/test_pixel_likelihood.py index 9a2b5bff49e..3ea62b8f855 100644 --- a/ctapipe/image/tests/test_pixel_likelihood.py +++ b/ctapipe/image/tests/test_pixel_likelihood.py @@ -2,6 +2,7 @@ from ctapipe.image import ( neg_log_likelihood, neg_log_likelihood_approx, + neg_log_likelihood_numeric, mean_poisson_likelihood_gaussian, chi_squared, mean_poisson_likelihood_full, @@ -56,11 +57,11 @@ 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]) @@ -68,15 +69,21 @@ def test_full_likelihood(): 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 + ) From db009423eda3a61f880614d84af799ec7f8f04b8 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Thu, 24 Aug 2023 16:20:37 +0200 Subject: [PATCH 03/12] Sort imports --- ctapipe/image/pixel_likelihood.py | 3 ++- ctapipe/image/tests/test_pixel_likelihood.py | 6 +++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py index fec2247d876..b2361e213a6 100644 --- a/ctapipe/image/pixel_likelihood.py +++ b/ctapipe/image/pixel_likelihood.py @@ -29,8 +29,9 @@ import numpy as np from scipy.integrate import quad -from scipy.stats import poisson from scipy.special import factorial +from scipy.stats import poisson + __all__ = [ "neg_log_likelihood_approx", diff --git a/ctapipe/image/tests/test_pixel_likelihood.py b/ctapipe/image/tests/test_pixel_likelihood.py index 3ea62b8f855..b5a8373557a 100644 --- a/ctapipe/image/tests/test_pixel_likelihood.py +++ b/ctapipe/image/tests/test_pixel_likelihood.py @@ -1,11 +1,11 @@ import numpy as np from ctapipe.image import ( + chi_squared, + mean_poisson_likelihood_full, + mean_poisson_likelihood_gaussian, neg_log_likelihood, neg_log_likelihood_approx, neg_log_likelihood_numeric, - mean_poisson_likelihood_gaussian, - chi_squared, - mean_poisson_likelihood_full, ) From 592d3146132b00363c65c78a5f9043dda24fe079 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Thu, 24 Aug 2023 16:23:19 +0200 Subject: [PATCH 04/12] Added changelog --- docs/changes/2388.bug.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/changes/2388.bug.rst diff --git a/docs/changes/2388.bug.rst b/docs/changes/2388.bug.rst new file mode 100644 index 00000000000..0d7d716c1aa --- /dev/null +++ b/docs/changes/2388.bug.rst @@ -0,0 +1 @@ +Fixed a bug in the calculation of the full numeric pixel likelihood and the corresponding tests. \ No newline at end of file From afd557add3dc04f49f69433291b884d8aa15db17 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Thu, 24 Aug 2023 16:56:01 +0200 Subject: [PATCH 05/12] Changed cofidence docstring --- ctapipe/image/pixel_likelihood.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py index b2361e213a6..0c36e58d7b4 100644 --- a/ctapipe/image/pixel_likelihood.py +++ b/ctapipe/image/pixel_likelihood.py @@ -121,8 +121,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: loat, 0 < x < 1 + Upper end of Poisson confidence interval of maximum prediction. + Determines upper end of poisson integration. Returns ------- From 7c0cf914b991e4736789fa7ff4b43aa3b5967a16 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Thu, 24 Aug 2023 17:12:38 +0200 Subject: [PATCH 06/12] black+isort --- ctapipe/image/pixel_likelihood.py | 1 - 1 file changed, 1 deletion(-) diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py index 0c36e58d7b4..cb47226d66b 100644 --- a/ctapipe/image/pixel_likelihood.py +++ b/ctapipe/image/pixel_likelihood.py @@ -32,7 +32,6 @@ from scipy.special import factorial from scipy.stats import poisson - __all__ = [ "neg_log_likelihood_approx", "neg_log_likelihood_numeric", From 77992f3ca096af8802fb11abcd74eaf42e06d21e Mon Sep 17 00:00:00 2001 From: gschwefer Date: Thu, 24 Aug 2023 17:12:38 +0200 Subject: [PATCH 07/12] black+isort --- ctapipe/image/pixel_likelihood.py | 1 - ctapipe/image/tests/test_pixel_likelihood.py | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py index 0c36e58d7b4..cb47226d66b 100644 --- a/ctapipe/image/pixel_likelihood.py +++ b/ctapipe/image/pixel_likelihood.py @@ -32,7 +32,6 @@ from scipy.special import factorial from scipy.stats import poisson - __all__ = [ "neg_log_likelihood_approx", "neg_log_likelihood_numeric", diff --git a/ctapipe/image/tests/test_pixel_likelihood.py b/ctapipe/image/tests/test_pixel_likelihood.py index b5a8373557a..cf2f1f2447c 100644 --- a/ctapipe/image/tests/test_pixel_likelihood.py +++ b/ctapipe/image/tests/test_pixel_likelihood.py @@ -1,4 +1,5 @@ import numpy as np + from ctapipe.image import ( chi_squared, mean_poisson_likelihood_full, From 76f881aa63224fd98b484c567b193caf536587b3 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Fri, 25 Aug 2023 08:12:00 +0200 Subject: [PATCH 08/12] Variable rename --- ctapipe/image/pixel_likelihood.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py index cb47226d66b..dff6cfa2270 100644 --- a/ctapipe/image/pixel_likelihood.py +++ b/ctapipe/image/pixel_likelihood.py @@ -135,11 +135,11 @@ def neg_log_likelihood_numeric( likelihood = epsilon - ns = np.arange(poisson(np.max(prediction)).ppf(confidence) + 1) + 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 From 327a49a57ef901a1d57f51aec3534dbc093a4385 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Tue, 29 Aug 2023 16:28:46 +0200 Subject: [PATCH 09/12] Fix typo and factor 2 in mean full mean likelihood --- ctapipe/image/pixel_likelihood.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py index dff6cfa2270..b02bbb30d4c 100644 --- a/ctapipe/image/pixel_likelihood.py +++ b/ctapipe/image/pixel_likelihood.py @@ -120,7 +120,7 @@ def neg_log_likelihood_numeric( Width of single p.e. peak (:math:`σ_γ`). pedestal: ndarray Width of pedestal (:math:`σ_p`). - confidence: loat, 0 < x < 1 + confidence: float, 0 < x < 1 Upper end of Poisson confidence interval of maximum prediction. Determines upper end of poisson integration. @@ -227,7 +227,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): From dedd1809cd0c966d84939b1c909086e563fff15f Mon Sep 17 00:00:00 2001 From: gschwefer Date: Tue, 29 Aug 2023 17:07:11 +0200 Subject: [PATCH 10/12] Add unit test for mean gaussian likelihood --- ctapipe/image/tests/test_pixel_likelihood.py | 28 +++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/ctapipe/image/tests/test_pixel_likelihood.py b/ctapipe/image/tests/test_pixel_likelihood.py index cf2f1f2447c..97561b8cc3c 100644 --- a/ctapipe/image/tests/test_pixel_likelihood.py +++ b/ctapipe/image/tests/test_pixel_likelihood.py @@ -24,7 +24,7 @@ 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) @@ -32,13 +32,35 @@ def test_mean_poisson_likelihoood_gaussian(): 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 From ed57223201edab96da357c398d5aef6905d2b70d Mon Sep 17 00:00:00 2001 From: gschwefer Date: Tue, 29 Aug 2023 17:11:06 +0200 Subject: [PATCH 11/12] Updated docstring --- ctapipe/image/pixel_likelihood.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py index b02bbb30d4c..1d124d63852 100644 --- a/ctapipe/image/pixel_likelihood.py +++ b/ctapipe/image/pixel_likelihood.py @@ -72,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 From 1d664a6a2999bce6cbad66ca11b19da5a176fffb Mon Sep 17 00:00:00 2001 From: gschwefer Date: Tue, 29 Aug 2023 17:11:06 +0200 Subject: [PATCH 12/12] Updated docstring --- ctapipe/image/pixel_likelihood.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py index b02bbb30d4c..2da4ee6c070 100644 --- a/ctapipe/image/pixel_likelihood.py +++ b/ctapipe/image/pixel_likelihood.py @@ -72,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