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 + )