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