diff --git a/autolens/analysis/analysis/lens.py b/autolens/analysis/analysis/lens.py
index 8a10aec33..70e54571a 100644
--- a/autolens/analysis/analysis/lens.py
+++ b/autolens/analysis/analysis/lens.py
@@ -157,8 +157,8 @@ def log_likelihood_penalty_from(
 
         Returns
         -------
-        The penalty value of the positions log likelihood, if the positions do not trace close in the source plane,
-        else a None is returned to indicate there is no penalty.
+        The summed penalty value of every positions log likelihood, if positions do not trace close in the
+        source plane, else an array-valued 0.0 to indicate there is no penalty.
         """
         log_likelihood_penalty = self._xp.array(0.0)
 
@@ -169,13 +169,12 @@ def log_likelihood_penalty_from(
                 if positions_likelihood is not None:
 
                     log_likelihood_penalty = (
-                        positions_likelihood.log_likelihood_penalty_from(
+                        log_likelihood_penalty
+                        + positions_likelihood.log_likelihood_penalty_from(
                             instance=instance, analysis=self, xp=self._xp
                         )
                     )
 
-                    log_likelihood_penalty += log_likelihood_penalty
-
             return log_likelihood_penalty
 
         return log_likelihood_penalty
diff --git a/test_autolens/imaging/model/test_analysis_imaging.py b/test_autolens/imaging/model/test_analysis_imaging.py
index 7549686c9..7a7ee230d 100644
--- a/test_autolens/imaging/model/test_analysis_imaging.py
+++ b/test_autolens/imaging/model/test_analysis_imaging.py
@@ -98,7 +98,13 @@ def test__positions__likelihood_overwrites__changes_likelihood(masked_imaging_7x
     )
     analysis_log_likelihood = analysis.log_likelihood_function(instance=instance)
 
-    assert analysis_log_likelihood == pytest.approx(-44097289521.734665, 1.0e-4)
+    analysis_log_likelihood_penalty = analysis.log_likelihood_penalty_from(instance=instance)
+    positions_log_likelihood_penalty = positions_likelihood.log_likelihood_penalty_from(
+        instance=instance, analysis=analysis
+    )
+
+    assert analysis_log_likelihood_penalty == pytest.approx(positions_log_likelihood_penalty, 1.0e-8)
+    assert analysis_log_likelihood == pytest.approx(-22048644768.175858, 1.0e-4)
 
 
 def test__positions__likelihood_overwrites__changes_likelihood__double_source_plane_example(masked_imaging_7x7):
@@ -123,7 +129,22 @@ def test__positions__likelihood_overwrites__changes_likelihood__double_source_pl
     )
     analysis_log_likelihood = analysis.log_likelihood_function(instance=instance)
 
-    assert analysis_log_likelihood == pytest.approx(-44097289521.734665, 1.0e-4)
+    analysis_log_likelihood_penalty = analysis.log_likelihood_penalty_from(instance=instance)
+    positions_log_likelihood_penalty_0 = positions_likelihood_0.log_likelihood_penalty_from(
+        instance=instance, analysis=analysis
+    )
+    positions_log_likelihood_penalty_1 = positions_likelihood_1.log_likelihood_penalty_from(
+        instance=instance, analysis=analysis
+    )
+
+    # The analysis penalty is the SUM of every entry in `positions_likelihood_list` -- a regression
+    # test against the historical bug where the loop returned 2x the last entry's penalty and
+    # discarded all earlier entries.
+    assert analysis_log_likelihood_penalty == pytest.approx(
+        positions_log_likelihood_penalty_0 + positions_log_likelihood_penalty_1, 1.0e-8
+    )
+    assert positions_log_likelihood_penalty_0 != pytest.approx(positions_log_likelihood_penalty_1, 1.0e-4)
+    assert analysis_log_likelihood == pytest.approx(-44140499627.74848, 1.0e-4)
 
 
 
diff --git a/test_autolens/interferometer/model/test_analysis_interferometer.py b/test_autolens/interferometer/model/test_analysis_interferometer.py
index c1bd2f4f4..8bfa06748 100644
--- a/test_autolens/interferometer/model/test_analysis_interferometer.py
+++ b/test_autolens/interferometer/model/test_analysis_interferometer.py
@@ -81,7 +81,7 @@ def test__positions__likelihood_overwrite__changes_likelihood(
     )
     analysis_log_likelihood = analysis.log_likelihood_function(instance=instance)
 
-    assert analysis_log_likelihood == pytest.approx(-44097289569.2342, 1.0e-4)
+    assert analysis_log_likelihood == pytest.approx(-22048644815.84869, 1.0e-4)
 
 
 def _pixelization_model():