Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 13 additions & 2 deletions autofit/mapper/prior/log_gaussian.py
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand All @@ -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)
31 changes: 23 additions & 8 deletions autofit/mapper/prior/log_uniform.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
"""
Expand Down
21 changes: 16 additions & 5 deletions autofit/messages/normal.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
"""
Expand Down
6 changes: 5 additions & 1 deletion test_autofit/mapper/model/test_model_mapper.py
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
69 changes: 38 additions & 31 deletions test_autofit/mapper/prior/test_prior.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
import math

import numpy as np
import pytest

import autofit as af
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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):
Expand All @@ -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
)
Loading