diff --git a/autolens/point/fit/positions/image/abstract.py b/autolens/point/fit/positions/image/abstract.py
index 9ce98701f..5109d5efc 100644
--- a/autolens/point/fit/positions/image/abstract.py
+++ b/autolens/point/fit/positions/image/abstract.py
@@ -15,6 +15,7 @@
   one observed position, useful for highly magnified systems where some images may be
   too close to separate.
 """
+
 from abc import ABC
 import numpy as np
 from typing import Optional
diff --git a/autolens/point/fit/positions/image/pair.py b/autolens/point/fit/positions/image/pair.py
index 9df5580d4..8c08e7a65 100644
--- a/autolens/point/fit/positions/image/pair.py
+++ b/autolens/point/fit/positions/image/pair.py
@@ -1,94 +1,108 @@
-import numpy as np
-from scipy.optimize import linear_sum_assignment
-
-import autoarray as aa
-
-from autolens.point.fit.positions.image.abstract import AbstractFitPositionsImagePair
-
-
-class FitPositionsImagePair(AbstractFitPositionsImagePair):
-    """
-    Fits the positions of a point source dataset using a `Tracer` object with an image-plane chi-squared where every
-    model position of the point-source is paired with its closest observed position, without allowing for repeated
-    pairings of the same observed position to model positions.
-
-    By not allowing for repeated pairings, this can produce behaviour such as a model position not being paired to
-    its closest observed position, but instead being paired to a further observed position, if doing so
-    means that the overall distances of pairings are reduced.
-
-    THIS FIT CURRENTLY GIVES UNRELIABLE RESULTS, BECAUSE IT GOES TO SOLUTIONS WHERE THE NUMBER OF MODEL POSITIONS
-    IS BELOW THE NUMBER OF DATA POSITIONS, REDUCING THE CHI-SQUARED TO LOW VALUES. PYAUTOLENS SHOULD BE UPDATED TO
-    PENALIZE THIS BEHAVIOUR BEFORE THIS FIT CAN BE USED. THIS REISDUAL MAP PROPERTY MAY ALSO NEED TO BE EXTENDED
-    TO ACCOUNT FOR NOISE.
-
-    The fit performs the following steps:
-
-    1) Determine the source-plane centre of the point source, which could be a free model parameter or computed
-       as the barycenter of ray-traced positions in the source-plane, using name pairing (see below).
-
-    2) Determine the image-plane model positions using the `PointSolver` and the source-plane centre of the point
-       source (e.g. ray tracing triangles to and from  the image and source planes), including accounting for
-       multi-plane ray-tracing.
-
-    3) Pair each model position with the observed position, not allowing for repeated pairings of the same
-       observed position to model positions, to compute the `residual_map`. This may result in some observed
-       positions not being paired to their closest model position, if doing so reduces the overall distances of
-       pairings.
-
-    5) Compute the chi-squared of each position as the square of the residual divided by the RMS noise-map value.
-
-    6) Sum the chi-squared values to compute the overall log likelihood of the fit.
-
-    Point source fitting uses name pairing, whereby the `name` of the `Point` object is paired to the name of the
-    point source dataset to ensure that point source datasets are fitted to the correct point source.
-
-    This fit object is used in the `FitPointDataset` to perform position based fitting of a `PointDataset`,
-    which may also fit other components of the point dataset like fluxes or time delays.
-
-    When performing a `model-fit`via an `AnalysisPoint` object the `figure_of_merit` of this object
-    is called and returned in the `log_likelihood_function`.
-
-    Parameters
-    ----------
-    name
-        The name of the point source dataset which is paired to a `Point` profile.
-    data
-        The positions of the point source in the image-plane which are fitted.
-    noise_map
-        The noise-map of the positions which are used to compute the log likelihood of the positions.
-    tracer
-        The tracer of galaxies whose point source profile are used to fit the positions.
-    solver
-        Solves the lens equation in order to determine the image-plane positions of a point source by ray-tracing
-        triangles to and from the source-plane.
-    profile
-        Manually input the profile of the point source, which is used instead of the one extracted from the
-        tracer via name pairing if that profile is not found.
-    """
-
-    @property
-    def residual_map(self) -> aa.ArrayIrregular:
-        residual_map = []
-
-        cost_matrix = np.linalg.norm(
-            np.array(
-                self.data,
-            )[:, np.newaxis]
-            - np.array(
-                self.model_data.array,
-            ),
-            axis=2,
-        )
-
-        data_indexes, model_indexes = linear_sum_assignment(cost_matrix)
-
-        for data_index, model_index in zip(data_indexes, model_indexes):
-            distance = np.sqrt(
-                self.square_distance(
-                    self.data[data_index], self.model_data.array[model_index]
-                )
-            )
-
-            residual_map.append(distance)
-
-        return aa.ArrayIrregular(values=residual_map)
+import numpy as np
+from scipy.optimize import linear_sum_assignment
+
+import autoarray as aa
+
+from autolens.point.fit.positions.image.abstract import AbstractFitPositionsImagePair
+
+
+class FitPositionsImagePair(AbstractFitPositionsImagePair):
+    """
+    Fits the positions of a point source dataset using a `Tracer` object with an image-plane chi-squared where every
+    model position of the point-source is paired with its closest observed position, without allowing for repeated
+    pairings of the same observed position to model positions.
+
+    By not allowing for repeated pairings, this can produce behaviour such as a model position not being paired to
+    its closest observed position, but instead being paired to a further observed position, if doing so
+    means that the overall distances of pairings are reduced.
+
+    **Under-prediction penalty**: the Hungarian assignment pairs ``min(n_observed, n_model)`` positions, which
+    historically meant that a model predicting *fewer* images than observed silently dropped the unmatched
+    observed positions from the chi-squared — rewarding under-prediction (samplers drove n_model below
+    n_observed to shrink the chi-squared). Unmatched observed positions now contribute their distance to the
+    nearest model position as a residual, and if the model predicts no images at all every observed position
+    contributes the ``no_image_residual`` floor. ``FitPositionsImagePairRepeat`` remains the model-fit default;
+    it additionally offers over-prediction policies.
+
+    The fit performs the following steps:
+
+    1) Determine the source-plane centre of the point source, which could be a free model parameter or computed
+       as the barycenter of ray-traced positions in the source-plane, using name pairing (see below).
+
+    2) Determine the image-plane model positions using the `PointSolver` and the source-plane centre of the point
+       source (e.g. ray tracing triangles to and from  the image and source planes), including accounting for
+       multi-plane ray-tracing.
+
+    3) Pair each model position with the observed position, not allowing for repeated pairings of the same
+       observed position to model positions, to compute the `residual_map`. This may result in some observed
+       positions not being paired to their closest model position, if doing so reduces the overall distances of
+       pairings.
+
+    5) Compute the chi-squared of each position as the square of the residual divided by the RMS noise-map value.
+
+    6) Sum the chi-squared values to compute the overall log likelihood of the fit.
+
+    Point source fitting uses name pairing, whereby the `name` of the `Point` object is paired to the name of the
+    point source dataset to ensure that point source datasets are fitted to the correct point source.
+
+    This fit object is used in the `FitPointDataset` to perform position based fitting of a `PointDataset`,
+    which may also fit other components of the point dataset like fluxes or time delays.
+
+    When performing a `model-fit`via an `AnalysisPoint` object the `figure_of_merit` of this object
+    is called and returned in the `log_likelihood_function`.
+
+    Parameters
+    ----------
+    name
+        The name of the point source dataset which is paired to a `Point` profile.
+    data
+        The positions of the point source in the image-plane which are fitted.
+    noise_map
+        The noise-map of the positions which are used to compute the log likelihood of the positions.
+    tracer
+        The tracer of galaxies whose point source profile are used to fit the positions.
+    solver
+        Solves the lens equation in order to determine the image-plane positions of a point source by ray-tracing
+        triangles to and from the source-plane.
+    profile
+        Manually input the profile of the point source, which is used instead of the one extracted from the
+        tracer via name pairing if that profile is not found.
+    """
+
+    no_image_residual = 1.0e4
+
+    @property
+    def residual_map(self) -> aa.ArrayIrregular:
+
+        model = np.asarray(self.model_data.array)
+        model = model[np.isfinite(model).all(axis=1)]
+
+        if model.shape[0] == 0:
+            return aa.ArrayIrregular(
+                values=[float(self.no_image_residual)] * len(self.data)
+            )
+
+        cost_matrix = np.linalg.norm(
+            np.array(
+                self.data,
+            )[:, np.newaxis]
+            - model,
+            axis=2,
+        )
+
+        data_indexes, model_indexes = linear_sum_assignment(cost_matrix)
+
+        # Residuals are ordered by observed position (matching the noise-map ordering): assigned
+        # positions get their Hungarian pairing distance; positions the assignment could not pair
+        # (under-prediction, n_model < n_observed) get their distance to the nearest model
+        # position, so a model that cannot produce an observed image is penalized, not rewarded.
+        residuals = np.min(cost_matrix, axis=1)
+
+        for data_index, model_index in zip(data_indexes, model_indexes):
+            residuals[data_index] = np.sqrt(
+                self.square_distance(self.data[data_index], model[model_index])
+            )
+
+        residual_map = [float(r) for r in residuals]
+
+        return aa.ArrayIrregular(values=residual_map)
diff --git a/autolens/point/fit/positions/image/pair_repeat.py b/autolens/point/fit/positions/image/pair_repeat.py
index e8d247531..52a1c3787 100644
--- a/autolens/point/fit/positions/image/pair_repeat.py
+++ b/autolens/point/fit/positions/image/pair_repeat.py
@@ -1,66 +1,219 @@
-import autoarray as aa
-
-from autolens.point.fit.positions.image.abstract import AbstractFitPositionsImagePair
-
-
-class FitPositionsImagePairRepeat(AbstractFitPositionsImagePair):
-    """
-    Fits the positions of a a point source dataset using a `Tracer` object with an image-plane chi-squared where every
-    model position of the point-source is paired with its closest observed position, allowing for repeated pairings of
-    the same observed position to model positions.
-
-    The fit performs the following steps:
-
-    1) Determine the source-plane centre of the point source, which could be a free model parameter or computed
-       as the barycenter of ray-traced positions in the source-plane, using name pairing (see below).
-
-    2) Determine the image-plane model positions using the `PointSolver` and the source-plane centre of the point
-       source (e.g. ray tracing triangles to and from  the image and source planes), including accounting for
-       multi-plane ray-tracing.
-
-    3) Pair each model position with the closest observed position, allowing for repeated pairings of the same
-       observed position to model positions, to compute the `residual_map`.
-
-    5) Compute the chi-squared of each position as the square of the residual divided by the RMS noise-map value.
-
-    6) Sum the chi-squared values to compute the overall log likelihood of the fit.
-
-    Point source fitting uses name pairing, whereby the `name` of the `Point` object is paired to the name of the
-    point source dataset to ensure that point source datasets are fitted to the correct point source.
-
-    This fit object is used in the `FitPointDataset` to perform position based fitting of a `PointDataset`,
-    which may also fit other components of the point dataset like fluxes or time delays.
-
-    When performing a `model-fit`via an `AnalysisPoint` object the `figure_of_merit` of this object
-    is called and returned in the `log_likelihood_function`.
-
-    Parameters
-    ----------
-    name
-        The name of the point source dataset which is paired to a `Point` profile.
-    data
-        The positions of the point source in the image-plane which are fitted.
-    noise_map
-        The noise-map of the positions which are used to compute the log likelihood of the positions.
-    tracer
-        The tracer of galaxies whose point source profile are used to fit the positions.
-    solver
-        Solves the lens equation in order to determine the image-plane positions of a point source by ray-tracing
-        triangles to and from the source-plane.
-    profile
-        Manually input the profile of the point source, which is used instead of the one extracted from the
-        tracer via name pairing if that profile is not found.
-    """
-
-    @property
-    def residual_map(self) -> aa.ArrayIrregular:
-        residual_map = []
-
-        for position in self.data:
-            distances = [
-                self.square_distance(model_position, position)
-                for model_position in self.model_data.array
-            ]
-            residual_map.append(self._xp.sqrt(self._xp.min(self._xp.array(distances))))
-
-        return aa.ArrayIrregular(values=self._xp.array(residual_map))
+import numpy as np
+
+import autoarray as aa
+import autogalaxy as ag
+
+from autolens.point.fit.positions.image.abstract import AbstractFitPositionsImagePair
+
+
+class FitPositionsImagePairRepeat(AbstractFitPositionsImagePair):
+    """
+    Fits the positions of a a point source dataset using a `Tracer` object with an image-plane chi-squared where every
+    model position of the point-source is paired with its closest observed position, allowing for repeated pairings of
+    the same observed position to model positions.
+
+    The fit performs the following steps:
+
+    1) Determine the source-plane centre of the point source, which could be a free model parameter or computed
+       as the barycenter of ray-traced positions in the source-plane, using name pairing (see below).
+
+    2) Determine the image-plane model positions using the `PointSolver` and the source-plane centre of the point
+       source (e.g. ray tracing triangles to and from  the image and source planes), including accounting for
+       multi-plane ray-tracing.
+
+    3) Pair each observed position with the closest model position, allowing for repeated pairings of the same
+       model position to observed positions, to compute the `residual_map`.
+
+    4) Apply the **unmatched-model policy** (see below): model-predicted images which no observed position paired
+       to — the over-prediction case, e.g. spurious images from a wrong mass model — may contribute an additional
+       chi-squared penalty.
+
+    5) Compute the chi-squared of each position as the square of the residual divided by the RMS noise-map value,
+       plus any policy penalty, and sum to compute the overall log likelihood of the fit.
+
+    **Over- and under-prediction**
+
+    A model tracer may predict more images than observed (over-prediction) or fewer (under-prediction):
+
+    - **Under-prediction** is always penalized: every *observed* position contributes a residual to its nearest
+      surviving model image, so a model that cannot produce an observed image pays the full image-plane distance
+      to the images it does produce. If the solver returns no images at all, every observed position contributes
+      the ``no_image_residual`` floor (a large, finite residual — loudly bad without breaking the sampler with
+      NaNs / infs).
+
+    - **Over-prediction** is controlled by the ``unmatched_model_policy`` class attribute:
+
+      - ``"ignore"`` — extra model images contribute nothing (the historical behaviour, now explicit).
+      - ``"penalize"`` — every finite model image which is not the nearest neighbour of any observed position
+        contributes its distance to the nearest observed position as an additional residual.
+      - ``"magnification_filter"`` (default) — as ``"penalize"``, but model images whose absolute magnification
+        is below ``magnification_threshold`` are exempt first. This is the standard observational convention
+        (e.g. Lenstool practice): a strongly demagnified extra image — typically the central image — is assumed
+        to be below the detection limit and does not count against the model, while a bright unobserved image
+        does.
+
+    The policy attributes are **class attributes** (the ``Analysis.LATENT_BATCH_MODE`` pattern), so a model-fit
+    can switch policy without new constructor plumbing::
+
+        class FitQuiet(al.FitPositionsImagePairRepeat):
+            unmatched_model_policy = "ignore"
+
+        al.AnalysisPoint(dataset=dataset, solver=solver, fit_positions_cls=FitQuiet)
+
+    Penalty terms are normalized by the mean of the position noise-map, and all policy computations are
+    fixed-shape (NaN-padded model positions are masked, never dropped), so the fit remains JAX-compilable.
+
+    Point source fitting uses name pairing, whereby the `name` of the `Point` object is paired to the name of the
+    point source dataset to ensure that point source datasets are fitted to the correct point source.
+
+    This fit object is used in the `FitPointDataset` to perform position based fitting of a `PointDataset`,
+    which may also fit other components of the point dataset like fluxes or time delays.
+
+    When performing a `model-fit`via an `AnalysisPoint` object the `figure_of_merit` of this object
+    is called and returned in the `log_likelihood_function`.
+
+    Parameters
+    ----------
+    name
+        The name of the point source dataset which is paired to a `Point` profile.
+    data
+        The positions of the point source in the image-plane which are fitted.
+    noise_map
+        The noise-map of the positions which are used to compute the log likelihood of the positions.
+    tracer
+        The tracer of galaxies whose point source profile are used to fit the positions.
+    solver
+        Solves the lens equation in order to determine the image-plane positions of a point source by ray-tracing
+        triangles to and from the source-plane.
+    profile
+        Manually input the profile of the point source, which is used instead of the one extracted from the
+        tracer via name pairing if that profile is not found.
+    """
+
+    unmatched_model_policy = "magnification_filter"
+    magnification_threshold = 0.1
+    no_image_residual = 1.0e4
+
+    @property
+    def _distance_matrix(self):
+        """
+        The (n_observed, n_model) matrix of distances between every observed and model position, with
+        non-finite (NaN-padded) model positions set to the ``no_image_residual`` sentinel so they are
+        never selected as a nearest neighbour.
+        """
+        data = self._xp.asarray(np.asarray(self.data))
+        model = self.model_data.array
+
+        distances = self._xp.sqrt(
+            self._xp.sum(
+                (data[:, None, :] - model[None, :, :]) ** 2,
+                axis=2,
+            )
+        )
+        finite = self._xp.isfinite(distances)
+        return self._xp.where(finite, distances, self.no_image_residual)
+
+    @property
+    def residual_map(self) -> aa.ArrayIrregular:
+        """
+        The distance of every observed position to its nearest model position (repeats allowed).
+
+        If the solver returns no images (or only NaN-padded rows), every observed position contributes the
+        ``no_image_residual`` floor.
+        """
+        model = self.model_data.array
+
+        if model.shape[0] == 0:
+            return aa.ArrayIrregular(
+                values=self._xp.full(len(self.data), float(self.no_image_residual))
+            )
+
+        return aa.ArrayIrregular(values=self._xp.min(self._distance_matrix, axis=1))
+
+    @property
+    def unmatched_model_mask(self):
+        """
+        Boolean mask over the model positions marking the "extra" images the policy penalizes: finite model
+        positions which are not the nearest neighbour of any observed position and (under
+        ``"magnification_filter"``) whose absolute magnification is above ``magnification_threshold``.
+        """
+        model = self.model_data.array
+
+        finite = self._xp.isfinite(model).all(axis=1)
+
+        nearest_indexes = self._xp.argmin(self._distance_matrix, axis=1)
+        matched = (
+            self._xp.arange(model.shape[0])[None, :] == nearest_indexes[:, None]
+        ).any(axis=0)
+
+        unmatched = self._xp.logical_and(finite, self._xp.logical_not(matched))
+
+        if self.unmatched_model_policy == "magnification_filter":
+            use_multi_plane = len(self.tracer.planes) > 2
+            plane_j = (
+                self.tracer.extract_plane_index_of_profile(profile_name=self.name)
+                if use_multi_plane
+                else -1
+            )
+            lens_calc = ag.LensCalc.from_tracer(
+                tracer=self.tracer,
+                use_multi_plane=use_multi_plane,
+                plane_j=plane_j,
+            )
+            safe_model = self._xp.where(
+                finite[:, None], model, self._xp.zeros_like(model)
+            )
+            magnifications = lens_calc.magnification_2d_via_hessian_from(
+                grid=safe_model, xp=self._xp
+            )
+            magnifications = (
+                magnifications.array
+                if hasattr(magnifications, "array")
+                else magnifications
+            )
+            detectable = self._xp.abs(magnifications) >= self.magnification_threshold
+            unmatched = self._xp.logical_and(unmatched, detectable)
+
+        return unmatched
+
+    @property
+    def unmatched_model_penalty_map(self):
+        """
+        The penalty residual of every model position: its distance to the nearest observed position where the
+        ``unmatched_model_mask`` is set, zero elsewhere.
+        """
+        distances_to_data = self._xp.min(self._distance_matrix, axis=0)
+        return self._xp.where(
+            self.unmatched_model_mask,
+            distances_to_data,
+            self._xp.zeros_like(distances_to_data),
+        )
+
+    @property
+    def n_unmatched_model_positions(self) -> int:
+        """
+        The number of model positions the policy counts as unexplained extras — a useful fit-quality
+        diagnostic alongside the chi-squared.
+        """
+        return self._xp.sum(self.unmatched_model_mask)
+
+    @property
+    def chi_squared(self) -> float:
+        """
+        The chi-squared of the fit: the observed-position chi-squared (as in all `AbstractFit` objects) plus,
+        for policies other than ``"ignore"``, the unmatched-model penalty terms normalized by the mean position
+        noise.
+        """
+        chi_squared = self._xp.sum(self.chi_squared_map.array)
+
+        if (
+            self.unmatched_model_policy == "ignore"
+            or self.model_data.array.shape[0] == 0
+        ):
+            return chi_squared
+
+        noise_mean = self._xp.mean(self._xp.asarray(np.asarray(self.noise_map)))
+
+        return chi_squared + self._xp.sum(
+            (self.unmatched_model_penalty_map / noise_mean) ** 2.0
+        )
diff --git a/test_autolens/point/fit/positions/image/test_pair.py b/test_autolens/point/fit/positions/image/test_pair.py
index cd92e5531..076336da9 100644
--- a/test_autolens/point/fit/positions/image/test_pair.py
+++ b/test_autolens/point/fit/positions/image/test_pair.py
@@ -25,3 +25,52 @@ def test__fit_positions_image_pair__three_observed_positions__model_not_repeated
 
     assert fit.model_data.in_list == [(4.0, 0.0), (3.0, 0.0), (0.0, 0.0)]
     assert fit.residual_map.in_list == [1.0, 3.0, 0.0]
+
+
+def test__under_prediction__unmatched_observed_positions_are_penalized():
+    # 3 observed, 1 model: the Hungarian assignment pairs one observed position; the other
+    # two historically dropped out of the chi-squared entirely (rewarding under-prediction).
+    # They now contribute their distance to the nearest model position, in data order.
+    point = al.ps.Point(centre=(0.1, 0.1))
+    galaxy = al.Galaxy(redshift=1.0, point_0=point)
+    tracer = al.Tracer(galaxies=[al.Galaxy(redshift=0.5), galaxy])
+
+    data = al.Grid2DIrregular([(0.0, 0.0), (3.0, 4.0), (6.0, 8.0)])
+    noise_map = al.ArrayIrregular([1.0, 1.0, 1.0])
+    model_data = al.Grid2DIrregular([(0.0, 0.0)])
+
+    solver = al.m.MockPointSolver(model_positions=model_data)
+
+    fit = al.FitPositionsImagePair(
+        name="point_0",
+        data=data,
+        noise_map=noise_map,
+        tracer=tracer,
+        solver=solver,
+    )
+
+    assert fit.residual_map.in_list == [0.0, 5.0, 10.0]
+    assert fit.chi_squared == pytest.approx(125.0, 1.0e-4)
+
+
+def test__under_prediction__no_model_images_hits_finite_floor():
+    point = al.ps.Point(centre=(0.1, 0.1))
+    galaxy = al.Galaxy(redshift=1.0, point_0=point)
+    tracer = al.Tracer(galaxies=[al.Galaxy(redshift=0.5), galaxy])
+
+    data = al.Grid2DIrregular([(0.0, 0.0), (3.0, 4.0)])
+    noise_map = al.ArrayIrregular([0.5, 1.0])
+    model_data = al.Grid2DIrregular(np.zeros(shape=(0, 2)))
+
+    solver = al.m.MockPointSolver(model_positions=model_data)
+
+    fit = al.FitPositionsImagePair(
+        name="point_0",
+        data=data,
+        noise_map=noise_map,
+        tracer=tracer,
+        solver=solver,
+    )
+
+    assert fit.residual_map.in_list == [1.0e4, 1.0e4]
+    assert np.isfinite(float(fit.log_likelihood))
diff --git a/test_autolens/point/fit/positions/image/test_pair_repeat.py b/test_autolens/point/fit/positions/image/test_pair_repeat.py
index 1d652b52c..91409c505 100644
--- a/test_autolens/point/fit/positions/image/test_pair_repeat.py
+++ b/test_autolens/point/fit/positions/image/test_pair_repeat.py
@@ -60,3 +60,111 @@ def test__fit_positions_image_pair_repeat__three_observed_positions__model_alloc
 
     assert fit.model_data.in_list == [(3.0, 1.0), (3.0, 4.0)]
     assert fit.residual_map.in_list == [np.sqrt(10.0), 0.0, 0.0]
+
+
+def test__over_prediction__unmatched_bright_model_image_is_penalized():
+    # 2 observed, 3 model: the third model image is far from every observed position and
+    # (with a no-mass tracer, magnification 1 everywhere) counts as a detectable extra.
+    point = al.ps.Point(centre=(0.1, 0.1))
+    galaxy = al.Galaxy(redshift=1.0, point_0=point)
+    tracer = al.Tracer(galaxies=[al.Galaxy(redshift=0.5), galaxy])
+
+    data = al.Grid2DIrregular([(0.0, 0.0), (3.0, 4.0)])
+    noise_map = al.ArrayIrregular([0.5, 1.0])
+    model_data = al.Grid2DIrregular([(0.0, 0.0), (3.0, 4.0), (10.0, 10.0)])
+
+    solver = al.m.MockPointSolver(model_positions=model_data)
+
+    fit = al.FitPositionsImagePairRepeat(
+        name="point_0",
+        data=data,
+        noise_map=noise_map,
+        tracer=tracer,
+        solver=solver,
+    )
+
+    # Matched residuals are zero; the whole chi-squared is the extra-image penalty:
+    # distance from (10, 10) to nearest observed (3, 4) over the mean noise (0.75).
+    assert int(fit.n_unmatched_model_positions) == 1
+    penalty_distance = np.sqrt(7.0**2 + 6.0**2)
+    assert fit.chi_squared == pytest.approx((penalty_distance / 0.75) ** 2, 1.0e-4)
+
+    class FitIgnore(al.FitPositionsImagePairRepeat):
+        unmatched_model_policy = "ignore"
+
+    fit_ignore = FitIgnore(
+        name="point_0",
+        data=data,
+        noise_map=noise_map,
+        tracer=tracer,
+        solver=solver,
+    )
+
+    assert fit_ignore.chi_squared == pytest.approx(0.0, abs=1.0e-8)
+
+
+def test__over_prediction__demagnified_central_image_exempt_under_magnification_filter():
+    # An isothermal lens demagnifies positions near its centre far below the 0.1 threshold:
+    # under the default magnification_filter policy the extra central image is exempt, while
+    # the explicit "penalize" policy charges for it.
+    lens = al.Galaxy(
+        redshift=0.5,
+        mass=al.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=1.0),
+    )
+    point = al.ps.Point(centre=(0.0, 0.0))
+    source = al.Galaxy(redshift=1.0, point_0=point)
+    tracer = al.Tracer(galaxies=[lens, source])
+
+    data = al.Grid2DIrregular([(0.0, 1.05), (0.0, -0.95)])
+    noise_map = al.ArrayIrregular([0.5, 0.5])
+    model_data = al.Grid2DIrregular([(0.0, 1.05), (0.0, -0.95), (0.0, 0.01)])
+
+    solver = al.m.MockPointSolver(model_positions=model_data)
+
+    fit = al.FitPositionsImagePairRepeat(
+        name="point_0",
+        data=data,
+        noise_map=noise_map,
+        tracer=tracer,
+        solver=solver,
+    )
+
+    assert int(fit.n_unmatched_model_positions) == 0
+    assert fit.chi_squared == pytest.approx(0.0, abs=1.0e-8)
+
+    class FitPenalize(al.FitPositionsImagePairRepeat):
+        unmatched_model_policy = "penalize"
+
+    fit_penalize = FitPenalize(
+        name="point_0",
+        data=data,
+        noise_map=noise_map,
+        tracer=tracer,
+        solver=solver,
+    )
+
+    assert int(fit_penalize.n_unmatched_model_positions) == 1
+    assert fit_penalize.chi_squared > 1.0
+
+
+def test__under_prediction__no_model_images_hits_finite_floor():
+    point = al.ps.Point(centre=(0.1, 0.1))
+    galaxy = al.Galaxy(redshift=1.0, point_0=point)
+    tracer = al.Tracer(galaxies=[al.Galaxy(redshift=0.5), galaxy])
+
+    data = al.Grid2DIrregular([(0.0, 0.0), (3.0, 4.0)])
+    noise_map = al.ArrayIrregular([0.5, 1.0])
+    model_data = al.Grid2DIrregular(np.zeros(shape=(0, 2)))
+
+    solver = al.m.MockPointSolver(model_positions=model_data)
+
+    fit = al.FitPositionsImagePairRepeat(
+        name="point_0",
+        data=data,
+        noise_map=noise_map,
+        tracer=tracer,
+        solver=solver,
+    )
+
+    assert fit.residual_map.in_list == [1.0e4, 1.0e4]
+    assert np.isfinite(float(fit.log_likelihood))