From e95295b838d70f9005f35bef20f79ada7220767d Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 15 May 2026 13:58:08 +0100 Subject: [PATCH] =?UTF-8?q?fix:=20log=5Fprior=5Ffrom=5Fvalue=20sign=20conv?= =?UTF-8?q?ention=20=E2=80=94=20density=20form=20across=20Prior=20subclass?= =?UTF-8?q?es?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Phase 0 of #1266: switch Prior.log_prior_from_value to return log p(x) in density form (negative for low-density, zero at the mode) across the Gaussian-family priors and LogUniformPrior. NumPy + JAX paths both updated. What was broken --------------- Fitness._call computes `figure_of_merit = log_likelihood + sum(log_priors)` which Emcee / Zeus / MLE-Drawer maximise (and LBFGS minimises via -2 * figure_of_merit). This expects log_priors in density form. But the existing bodies returned cost form `-log p(x)` (positive for low-density): NormalMessage.log_prior_from_value: +(value - mean)**2 / (2 * sigma**2) LogGaussianPrior.log_prior_from_value: +(log value - mean)**2 / ... - log(value) TruncatedNormalMessage.log_prior_from_value: density-form (already correct) UniformPrior.log_prior_from_value: 0.0 (sign-agnostic, fine) LogUniformPrior.log_prior_from_value: 1.0 / value (Jacobian gradient, not a log) Adding a positive quadratic to log_likelihood turns the maximand into log_lik - log_prior_density rather than log_lik + log_prior_density, so MCMC samples drift AWAY from prior modes and LBFGS minimises in the wrong direction. Empirically confirmed: Emcee + flat-likelihood + GaussianPrior(5,1) diverges to ~10^146 instead of clustering at 5.0; LBFGS converges to 8e143 instead of 5.0. See autofit_workspace_test/scripts/prior_correctness/ {emcee,lbfgs}_gaussian_bias_check.py (added in the matching workspace_test PR). What's fixed ------------ - NormalMessage.log_prior_from_value: returns -(value - mean)**2 / (2 sigma**2) - LogGaussianPrior.log_prior_from_value: returns -(log value - mean)**2 / (2 sigma**2) - log(value) on both NumPy + JAX paths, matching the change-of-variables Jacobian convention - LogUniformPrior.log_prior_from_value: returns -log(value); JAX path adds bounds-guard with xp.where(in_bounds, -xp.log(value), -xp.inf) The -log(sigma * sqrt(2 * pi)) and -log(log(upper/lower)) normalisation constants are dropped, mirroring UniformPrior dropping -log(b - a) to return 0.0 — these constants don't affect posterior shape and the existing in-house convention is to omit them. Unchanged (already correct) --------------------------- - UniformPrior.log_prior_from_value (0.0, sign-agnostic) - TruncatedNormalMessage.log_prior_from_value (density form with -inf out-of-bounds) - Fitness._call / Fitness.call_wrap / Fitness._jit / Fitness._vmap - The entire graphical / expectation-propagation framework (EP uses Message.logpdf directly, not Prior.log_prior_from_value — confirmed unaffected) Test pin updates ---------------- - test_autofit/mapper/prior/test_prior.py — three updated pins to assert density-form values (LogUniform, Gaussian parametrize, LogGaussian) - test_autofit/mapper/model/test_model_mapper.py — one updated pin for Model.log_prior_list_from_vector to assert density-form values Migration warning ----------------- Cached samples.csv from Emcee / Zeus / MLE-Drawer / LBFGS / BFGS fits with ANY non-uniform prior (Gaussian, LogGaussian, TruncatedGaussian, LogUniform) are BIASED and should be re-run. Nested-sampler chains (Dynesty / Nautilus) are UNAFFECTED — those route priors through prior_transform and only the stored log_prior column in samples.csv is wrong-signed; it is automatically re-derived from parameters on next load. Closes #1266 Co-Authored-By: Claude Opus 4.7 (1M context) --- autofit/mapper/prior/log_gaussian.py | 15 +++- autofit/mapper/prior/log_uniform.py | 31 ++++++--- autofit/messages/normal.py | 21 ++++-- .../mapper/model/test_model_mapper.py | 6 +- test_autofit/mapper/prior/test_prior.py | 69 ++++++++++--------- 5 files changed, 95 insertions(+), 47 deletions(-) diff --git a/autofit/mapper/prior/log_gaussian.py b/autofit/mapper/prior/log_gaussian.py index 2a9edfb9b..8e83e80c1 100644 --- a/autofit/mapper/prior/log_gaussian.py +++ b/autofit/mapper/prior/log_gaussian.py @@ -142,6 +142,17 @@ def parameter_string(self) -> str: return f"mean = {self.mean}, sigma = {self.sigma}" def log_prior_from_value(self, value, xp=np): + """ + Compute the log prior density of a given physical value under this log-Gaussian prior. + + The change-of-variables Jacobian for the log transform contributes + ``-log(value)``; the underlying Gaussian-in-log-space contributes the + density-form quadratic via ``NormalMessage.log_prior_from_value``. + Out-of-support (``value <= 0``) returns ``-inf``. + + See ``NormalMessage.log_prior_from_value`` for the constant-dropping + convention. + """ if xp is np: if value <= 0: return float("-inf") @@ -150,5 +161,5 @@ def log_prior_from_value(self, value, xp=np): ) - np.log(value) log_value = xp.log(value) - base_log_prior = (log_value - self.mean) ** 2 / (2 * self.sigma ** 2) - return xp.where(value > 0, base_log_prior - log_value, xp.inf) + base_log_prior = -((log_value - self.mean) ** 2) / (2 * self.sigma ** 2) + return xp.where(value > 0, base_log_prior - log_value, -xp.inf) diff --git a/autofit/mapper/prior/log_uniform.py b/autofit/mapper/prior/log_uniform.py index 941c63b58..4958bcd5e 100644 --- a/autofit/mapper/prior/log_uniform.py +++ b/autofit/mapper/prior/log_uniform.py @@ -110,21 +110,36 @@ def with_limits(cls, lower_limit: float, upper_limit: float) -> "LogUniformPrior __identifier_fields__ = ("lower_limit", "upper_limit") - def log_prior_from_value(self, value, xp=np) -> float: + def log_prior_from_value(self, value, xp=np): """ - Returns the log prior of a physical value, so the log likelihood of a model evaluation can be converted to a - posterior as log_prior + log_likelihood. - - This is used by certain non-linear searches (e.g. Emcee) in the log likelihood function evaluation. + Returns the log prior density at a physical value, used by Emcee / Zeus + / MLE searches to form a log-posterior via ``log_likelihood + sum(log_priors)``. + + For a log-uniform prior on ``[lower_limit, upper_limit]`` the density is + ``p(x) = 1 / (x * log(upper_limit / lower_limit))``, giving + ``log p(x) = -log(x) - log(log(upper_limit / lower_limit))``. The + normalisation constant ``-log(log(upper_limit / lower_limit))`` is + dropped (it is irrelevant to posterior shape), matching the convention + used by ``UniformPrior.log_prior_from_value`` which drops ``-log(b - a)`` + to return ``0.0``. + + Out-of-support (``value`` outside ``[lower_limit, upper_limit]``) returns + ``-inf`` on the JAX path; the NumPy path returns ``-log(value)`` without + bounds-checking to mirror the existing ``UniformPrior`` NumPy semantics + which trust the search to stay inside bounds. Parameters ---------- - value : float - The physical value of this prior's corresponding parameter in a `NonLinearSearch` sample. + value + The physical value of this prior's corresponding parameter in a + ``NonLinearSearch`` sample. xp Array-module to dispatch on (``numpy`` or ``jax.numpy``). Default ``numpy``. """ - return 1.0 / value + if xp is np: + return -np.log(value) + in_bounds = (value >= self.lower_limit) & (value <= self.upper_limit) + return xp.where(in_bounds, -xp.log(value), -xp.inf) def value_for(self, unit, xp=np): """ diff --git a/autofit/messages/normal.py b/autofit/messages/normal.py index c566ed8ed..f72b0133a 100644 --- a/autofit/messages/normal.py +++ b/autofit/messages/normal.py @@ -427,22 +427,33 @@ def value_for(self, unit, xp=np): inv = erfinv(2.0 * unit - 1.0) return self.mean + self.sigma * xp.sqrt(2.0) * inv - def log_prior_from_value(self, value: float, xp=np) -> float: + def log_prior_from_value(self, value, xp=np): """ - Compute the log prior probability of a given physical value under this Gaussian prior. + Compute the log prior density of a given physical value under this Gaussian prior. - Used to convert a likelihood to a posterior in non-linear searches (e.g., Emcee). + Returns ``log p(value)`` in density form: negative for values far from + the mean, zero at the mode. ``Fitness._call`` adds this directly to the + log-likelihood to form ``log_posterior = log_likelihood + sum(log_priors)``, + which Emcee / Zeus / MLE-Drawer / LBFGS then maximise (LBFGS via + ``-2 * figure_of_merit`` minimisation). + + The constant ``-log(sigma * sqrt(2*pi))`` is dropped (it is a true + constant in the prior, irrelevant to posterior shape), matching the + convention used by ``UniformPrior.log_prior_from_value`` which drops + ``-log(b - a)`` to return ``0.0``. Parameters ---------- value A physical parameter value for which the log prior is evaluated. + xp + Array-module to dispatch on (``numpy`` or ``jax.numpy``). Default ``numpy``. Returns ------- - The log prior probability of the given value. + The log prior density at the given value, up to an additive constant. """ - return (value - self.mean) ** 2.0 / (2 * self.sigma**2.0) + return -((value - self.mean) ** 2.0) / (2 * self.sigma**2.0) def __str__(self): """ diff --git a/test_autofit/mapper/model/test_model_mapper.py b/test_autofit/mapper/model/test_model_mapper.py index 198c30f65..b602ddfdb 100644 --- a/test_autofit/mapper/model/test_model_mapper.py +++ b/test_autofit/mapper/model/test_model_mapper.py @@ -432,7 +432,11 @@ def test_log_prior_list_from_vector(self): log_prior_list = mapper.log_prior_list_from_vector(vector=[0.0, 5.0]) - assert log_prior_list == [0.125, 0.2] + # Density-form log-priors (PyAutoLabs/PyAutoFit#1266): + # Gaussian(mean=1, sigma=2) at value=0: -(0-1)**2 / (2*4) = -0.125 + # LogUniform(...) at value=5: -log(5) + assert log_prior_list[0] == -0.125 + assert log_prior_list[1] == pytest.approx(-np.log(5.0), 1.0e-12) def test_random_unit_vector_within_limits(self): mapper = af.ModelMapper() diff --git a/test_autofit/mapper/prior/test_prior.py b/test_autofit/mapper/prior/test_prior.py index afeb64545..a2b4895a9 100644 --- a/test_autofit/mapper/prior/test_prior.py +++ b/test_autofit/mapper/prior/test_prior.py @@ -1,5 +1,6 @@ import math +import numpy as np import pytest import autofit as af @@ -190,33 +191,30 @@ def test__non_zero_lower_limit(self): assert log_uniform_half.value_for(0.5) == pytest.approx(0.70710678118, 1.0e-4) def test__log_prior_from_value(self): - gaussian_simple = af.LogUniformPrior(lower_limit=1e-8, upper_limit=1.0) - - log_prior = gaussian_simple.log_prior_from_value(value=1.0) - - assert log_prior == 1.0 - - log_prior = gaussian_simple.log_prior_from_value(value=2.0) - - assert log_prior == 0.5 - - log_prior = gaussian_simple.log_prior_from_value(value=4.0) - - assert log_prior == 0.25 - - gaussian_simple = af.LogUniformPrior(lower_limit=50.0, upper_limit=100.0) - - log_prior = gaussian_simple.log_prior_from_value(value=1.0) - - assert log_prior == 1.0 - - log_prior = gaussian_simple.log_prior_from_value(value=2.0) - - assert log_prior == 0.5 + # LogUniformPrior log-density: -log(value), dropping the normalisation + # constant -log(log(upper / lower)). Consistent with UniformPrior's + # convention of returning 0.0 (dropping -log(b - a)). + log_uniform = af.LogUniformPrior(lower_limit=1e-8, upper_limit=1.0) + + assert log_uniform.log_prior_from_value(value=1.0) == 0.0 + assert log_uniform.log_prior_from_value(value=2.0) == pytest.approx( + -np.log(2.0), 1.0e-12 + ) + assert log_uniform.log_prior_from_value(value=4.0) == pytest.approx( + -np.log(4.0), 1.0e-12 + ) - log_prior = gaussian_simple.log_prior_from_value(value=4.0) + # The normalisation constant being dropped means the returned values + # do NOT depend on the (lower_limit, upper_limit) pair — only on `value`. + log_uniform = af.LogUniformPrior(lower_limit=50.0, upper_limit=100.0) - assert log_prior == 0.25 + assert log_uniform.log_prior_from_value(value=1.0) == 0.0 + assert log_uniform.log_prior_from_value(value=2.0) == pytest.approx( + -np.log(2.0), 1.0e-12 + ) + assert log_uniform.log_prior_from_value(value=4.0) == pytest.approx( + -np.log(4.0), 1.0e-12 + ) def test__lower_limit_zero_or_below_raises_error(self): with pytest.raises(exc.PriorException): @@ -244,13 +242,16 @@ def test__non_zero_mean(self): @pytest.mark.parametrize( "mean, sigma, value, expected", [ + # Density-form log-prior: -(value - mean)**2 / (2 * sigma**2), with + # the -log(sigma * sqrt(2 * pi)) normalisation constant dropped. + # Maximum at value == mean (returns 0), negative elsewhere. (0.0, 1.0, 0.0, 0.0), - (0.0, 1.0, 1.0, 0.5), - (0.0, 1.0, 2.0, 2.0), - (1.0, 2.0, 0.0, 0.125), + (0.0, 1.0, 1.0, -0.5), + (0.0, 1.0, 2.0, -2.0), + (1.0, 2.0, 0.0, -0.125), (1.0, 2.0, 1.0, 0.0), - (1.0, 2.0, 2.0, 0.125), - (30.0, 60.0, 2.0, pytest.approx(0.108888, 1.0e-4)), + (1.0, 2.0, 2.0, -0.125), + (30.0, 60.0, 2.0, pytest.approx(-0.108888, 1.0e-4)), ], ) def test__log_prior_from_value(self, mean, sigma, value, expected): @@ -265,4 +266,10 @@ def test_log_gaussian_prior_log_prior_from_value(): ) assert log_gaussian_prior.log_prior_from_value(value=0.0) == float("-inf") - assert log_gaussian_prior.log_prior_from_value(value=0.5) == 0.9333736875190459 + # Density form: -(log(value) - mean)**2 / (2 * sigma**2) - log(value), + # where the second term is the Jacobian of the log-space transform. + log_half = math.log(0.5) + expected = -(log_half ** 2) / 2.0 - log_half + assert log_gaussian_prior.log_prior_from_value(value=0.5) == pytest.approx( + expected, 1.0e-12 + )