fix(graphical): EP statistics completion — F6 truncated KL + F7(b) evidence recording

autofit/graphical/README.md

@@ -185,15 +185,30 @@ The EP approximation to the model evidence factorises as:
     log ∫ ∏ₖ qₖ = Σₖ (log hₖ − A(ηₖ)) − ( log h − A(Σₖ ηₖ) )          (16)
 
 and the per-factor `Ẑₐ` is carried on `MeanField.log_norm` by the
-projection. (Audit note: as of 2026-07 the `log_norm` bookkeeping is
-broken in three places — issue #1332 finding F7 — so Eq. (15) should
-not be used for model comparison until those fixes land; the
-decomposition itself is correct.)
+projection: a search-driven factor update records the sampler's
+log-evidence of the tilted fit there (`AbstractSearch.optimise` →
+`MeanField.from_priors(..., log_norm=...)`).
+
+All three legs of the 2026-07 audit's finding F7 (#1332) are now fixed:
+(a) `MeanField.__truediv__`/`__pow__` propagate `log_norm` (#1351),
+(b) the search-driven projection records `Ẑₐ` (this section),
+(c) the truncated-normal log-partition is complete (#1345). Evidence-
+correct model comparison additionally requires **nested-sampling factor
+searches** — MCMC/MLE searches carry no evidence estimate and contribute
+`log_norm = 0`.
 
 ## 6. Deterministic variables
 
-Three composition mechanisms exist; they are **not** interchangeable
-and reconciling them is an open design item (EP review Phase 5):
+Three composition mechanisms exist; they are **not** interchangeable.
+**Decision (EP review Phase 5, #1336, 2026-07-10): keep all three —
+no unification, no deprecation.** The trade-off users should choose by:
+compound priors / shared variables are statistically *tighter* (the
+relation holds exactly inside each factor), while `factor_out` trades
+that exactness for modularity (the deterministic variable receives its
+own messages and `q(z)` factorises from its parents). The declarative
+surface for deterministic quantities is the explicit compound pattern
+(e.g. `model.sigma * 2.355`), pinned by the seam tests (§8); the
+`model.<property>` sugar from #1153 was deliberately retired.
 
 1. **Graph-level deterministic variables**: `Factor(..., factor_out=v)`
    declares outputs computed by the factor

autofit/graphical/mean_field.py

@@ -123,7 +123,9 @@ def fixed_values(self):
         )
 
     @classmethod
-    def from_priors(cls, priors: Iterable[Prior]) -> "MeanField":
+    def from_priors(
+        cls, priors: Iterable[Prior], log_norm: float = 0.0
+    ) -> "MeanField":
         """
         Create a MeanField from a list of priors.
 
@@ -134,12 +136,22 @@ def from_priors(cls, priors: Iterable[Prior]) -> "MeanField":
         ----------
         priors
             A list of priors
+        log_norm
+            The log-normalisation carried by this mean field. For a
+            search-driven factor update this is the sampler's log-evidence of
+            the tilted-distribution fit — the per-factor ``Ẑₐ`` of README §5
+            Eq. (14), which the projection is documented to carry on
+            ``MeanField.log_norm`` (#1332 F7(b)). Defaults to 0.0 for callers
+            with no evidence to record.
 
         Returns
         -------
         A mean field
         """
-        return MeanField({prior: prior.message for prior in priors})
+        return MeanField(
+            {prior: prior.message for prior in priors},
+            log_norm=log_norm,
+        )
 
     pop = dict.pop
     values = dict.values

autofit/messages/truncated_normal.py

@@ -304,11 +304,21 @@ def sample(self, n_samples: Optional[int] = None) -> np.ndarray:
 
     def kl(self, dist : "TruncatedNormalMessage") -> float:
         """
-        Compute the Kullback-Leibler (KL) divergence between two truncated Gaussian distributions.
+        The exact Kullback-Leibler divergence KL(self || dist) between two
+        truncated Gaussians sharing the same support (#1332 F6).
 
-        This is an approximate KL divergence that assumes both distributions are truncated Gaussians
-        with the same support (i.e. the same lower and upper limits). If the supports differ, this
-        expression is invalid and should raise an error or be corrected for normalization.
+        With p, q truncated to the same ``[a, b]``, standardised bounds
+        ``α = (a − μ)/σ``, ``β = (b − μ)/σ``, truncation mass
+        ``Z = Φ(β) − Φ(α)`` and the *truncated* moments ``m_p``, ``V_p`` of p:
+
+            KL(p‖q) = log(σ_q/σ_p) + log(Z_q/Z_p)
+                      + ½·[ (V_p + (m_p − μ_q)²)/σ_q² − (V_p + (m_p − μ_p)²)/σ_p² ]
+
+        As the bounds recede (Z → 1, m_p → μ_p, V_p → σ_p²) this reduces to the
+        untruncated Gaussian KL. Previously this method *used* the untruncated
+        formula, which degrades as posterior mass approaches the bounds —
+        directly distorting the ``EPHistory`` convergence metric for models
+        built on ``af.TruncatedGaussianPrior`` (e.g. HowToFit chapter 3).
 
         Parameters
         ----------
@@ -323,10 +333,29 @@ def kl(self, dist : "TruncatedNormalMessage") -> float:
         if (self.lower_limit != dist.lower_limit) or (self.upper_limit != dist.upper_limit):
             raise ValueError("KL divergence between truncated Gaussians with different support is not implemented.")
 
+        from scipy.stats import norm, truncnorm
+
+        a_p = (self.lower_limit - self.mean) / self.sigma
+        b_p = (self.upper_limit - self.mean) / self.sigma
+        a_q = (self.lower_limit - dist.mean) / dist.sigma
+        b_q = (self.upper_limit - dist.mean) / dist.sigma
+
+        log_Z_p = np.log(norm.cdf(b_p) - norm.cdf(a_p))
+        log_Z_q = np.log(norm.cdf(b_q) - norm.cdf(a_q))
+
+        # Truncated mean and variance of p (scipy returns them for the
+        # standardised bounds with loc/scale applied).
+        m_p, V_p = truncnorm.stats(
+            a_p, b_p, loc=self.mean, scale=self.sigma, moments="mv"
+        )
+
+        e_zp2 = (V_p + (m_p - self.mean) ** 2) / self.sigma**2
+        e_zq2 = (V_p + (m_p - dist.mean) ** 2) / dist.sigma**2
+
         return (
             np.log(dist.sigma / self.sigma)
-            + (self.sigma**2 + (self.mean - dist.mean) ** 2) / 2 / dist.sigma**2
-            - 1 / 2
+            + (log_Z_q - log_Z_p)
+            + 0.5 * (e_zq2 - e_zp2)
         )
 
     @classmethod

autofit/non_linear/search/abstract_search.py

@@ -367,7 +367,22 @@ def optimise(
 
         result = self.fit(model=model, analysis=analysis)
 
-        new_model_dist = MeanField.from_priors(result.projected_model.priors)
+        # Record the sampler's log-evidence of this tilted-distribution fit on
+        # the projected mean field — the per-factor Ẑₐ that README §5 documents
+        # `MeanField.log_norm` as carrying (#1332 F7(b)). Previously always 0,
+        # so `EPMeanField.log_evidence` could not be trusted for model
+        # comparison in sampler-driven EP fits. Searches with no evidence
+        # estimate (MCMC / MLE) yield None and keep the 0.0 default — evidence-
+        # correct model comparison requires nested-sampling factor searches.
+        # (Both levels guarded: e.g. StaticResult carries no samples at all.)
+        log_evidence = getattr(
+            getattr(result, "samples", None), "log_evidence", None
+        )
+
+        new_model_dist = MeanField.from_priors(
+            result.projected_model.priors,
+            log_norm=log_evidence if log_evidence is not None else 0.0,
+        )
 
         status.result = result
 

test_autofit/graphical/functionality/test_ep_statistics_completion.py

@@ -0,0 +1,99 @@
+"""Regression tests for the EP statistics completion (PyAutoFit #1353;
+findings F6 and F7(b) from the #1332 audit). Numpy-only.
+"""
+import numpy as np
+import pytest
+from scipy.stats import truncnorm
+
+import autofit as af
+from autofit.graphical.mean_field import MeanField
+from autofit.messages.truncated_normal import TruncatedNormalMessage
+
+
+# --- F6: exact same-support truncated KL via truncated moments ---
+
+def _mc_kl_truncnorm(p, q, n=400_000, seed=3):
+    a_p = (p.lower_limit - p.mean) / p.sigma
+    b_p = (p.upper_limit - p.mean) / p.sigma
+    a_q = (q.lower_limit - q.mean) / q.sigma
+    b_q = (q.upper_limit - q.mean) / q.sigma
+    rng = np.random.default_rng(seed)
+    x = truncnorm.rvs(
+        a_p, b_p, loc=p.mean, scale=p.sigma, size=n, random_state=rng
+    )
+    return float(
+        np.mean(
+            truncnorm.logpdf(x, a_p, b_p, loc=p.mean, scale=p.sigma)
+            - truncnorm.logpdf(x, a_q, b_q, loc=q.mean, scale=q.sigma)
+        )
+    )
+
+
+def _old_untruncated_kl(p, q):
+    return (
+        np.log(q.sigma / p.sigma)
+        + (p.sigma**2 + (p.mean - q.mean) ** 2) / 2 / q.sigma**2
+        - 1 / 2
+    )
+
+
+@pytest.mark.parametrize(
+    "p, q",
+    [
+        # comfortable interior — old approximation was nearly right here
+        (
+            TruncatedNormalMessage(0.0, 1.0, -10.0, 10.0),
+            TruncatedNormalMessage(0.5, 2.0, -10.0, 10.0),
+        ),
+        # mass pressed against the bounds — old approximation degrades badly
+        (
+            TruncatedNormalMessage(0.9, 1.5, -1.0, 1.0),
+            TruncatedNormalMessage(-0.5, 0.7, -1.0, 1.0),
+        ),
+        # half-bounded, mean below the bound (prior-passing shape)
+        (
+            TruncatedNormalMessage(-0.5, 1.0, 0.0, np.inf),
+            TruncatedNormalMessage(1.0, 2.0, 0.0, np.inf),
+        ),
+    ],
+    ids=["interior", "near-bounds", "half-bounded"],
+)
+def test__truncated_kl_matches_monte_carlo(p, q):
+    analytic = float(p.kl(q))
+    mc = _mc_kl_truncnorm(p, q)
+    assert analytic == pytest.approx(mc, abs=0.02)
+
+
+def test__truncated_kl_reduces_to_gaussian_for_wide_bounds():
+    p = TruncatedNormalMessage(0.0, 1.0, -50.0, 50.0)
+    q = TruncatedNormalMessage(1.0, 2.0, -50.0, 50.0)
+    assert float(p.kl(q)) == pytest.approx(float(_old_untruncated_kl(p, q)), rel=1e-9)
+
+
+def test__truncated_kl_near_bounds_differs_from_old_formula():
+    # The case the audit quantified (errors 1.5% -> 140% as mass reaches the
+    # bounds): the exact value must both match MC and differ measurably from
+    # the old untruncated formula.
+    p = TruncatedNormalMessage(0.9, 1.5, -1.0, 1.0)
+    q = TruncatedNormalMessage(-0.5, 0.7, -1.0, 1.0)
+    exact = float(p.kl(q))
+    old = float(_old_untruncated_kl(p, q))
+    assert abs(exact - old) / abs(exact) > 0.05
+
+
+def test__truncated_kl_different_support_still_raises():
+    p = TruncatedNormalMessage(0.0, 1.0, -1.0, 1.0)
+    q = TruncatedNormalMessage(0.0, 1.0, -2.0, 2.0)
+    with pytest.raises(ValueError):
+        p.kl(q)
+
+
+# --- F7(b): from_priors records the per-factor evidence on log_norm ---
+
+def test__from_priors_log_norm_default_and_explicit():
+    priors = [af.GaussianPrior(0.0, 1.0), af.GaussianPrior(1.0, 2.0)]
+    assert MeanField.from_priors(priors).log_norm == 0.0
+    mf = MeanField.from_priors(priors, log_norm=-42.5)
+    assert mf.log_norm == -42.5
+    # and it survives the (now-fixed, #1351) operator algebra
+    assert (mf / MeanField.from_priors(priors, log_norm=-2.5)).log_norm == -40.0