Skip to content

Commit e95295b

Browse files
Jammy2211claude
authored andcommitted
fix: log_prior_from_value sign convention — density form across Prior subclasses
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) <noreply@anthropic.com>
1 parent 377ca1b commit e95295b

5 files changed

Lines changed: 95 additions & 47 deletions

File tree

autofit/mapper/prior/log_gaussian.py

Lines changed: 13 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -142,6 +142,17 @@ def parameter_string(self) -> str:
142142
return f"mean = {self.mean}, sigma = {self.sigma}"
143143

144144
def log_prior_from_value(self, value, xp=np):
145+
"""
146+
Compute the log prior density of a given physical value under this log-Gaussian prior.
147+
148+
The change-of-variables Jacobian for the log transform contributes
149+
``-log(value)``; the underlying Gaussian-in-log-space contributes the
150+
density-form quadratic via ``NormalMessage.log_prior_from_value``.
151+
Out-of-support (``value <= 0``) returns ``-inf``.
152+
153+
See ``NormalMessage.log_prior_from_value`` for the constant-dropping
154+
convention.
155+
"""
145156
if xp is np:
146157
if value <= 0:
147158
return float("-inf")
@@ -150,5 +161,5 @@ def log_prior_from_value(self, value, xp=np):
150161
) - np.log(value)
151162

152163
log_value = xp.log(value)
153-
base_log_prior = (log_value - self.mean) ** 2 / (2 * self.sigma ** 2)
154-
return xp.where(value > 0, base_log_prior - log_value, xp.inf)
164+
base_log_prior = -((log_value - self.mean) ** 2) / (2 * self.sigma ** 2)
165+
return xp.where(value > 0, base_log_prior - log_value, -xp.inf)

autofit/mapper/prior/log_uniform.py

Lines changed: 23 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -110,21 +110,36 @@ def with_limits(cls, lower_limit: float, upper_limit: float) -> "LogUniformPrior
110110

111111
__identifier_fields__ = ("lower_limit", "upper_limit")
112112

113-
def log_prior_from_value(self, value, xp=np) -> float:
113+
def log_prior_from_value(self, value, xp=np):
114114
"""
115-
Returns the log prior of a physical value, so the log likelihood of a model evaluation can be converted to a
116-
posterior as log_prior + log_likelihood.
117-
118-
This is used by certain non-linear searches (e.g. Emcee) in the log likelihood function evaluation.
115+
Returns the log prior density at a physical value, used by Emcee / Zeus
116+
/ MLE searches to form a log-posterior via ``log_likelihood + sum(log_priors)``.
117+
118+
For a log-uniform prior on ``[lower_limit, upper_limit]`` the density is
119+
``p(x) = 1 / (x * log(upper_limit / lower_limit))``, giving
120+
``log p(x) = -log(x) - log(log(upper_limit / lower_limit))``. The
121+
normalisation constant ``-log(log(upper_limit / lower_limit))`` is
122+
dropped (it is irrelevant to posterior shape), matching the convention
123+
used by ``UniformPrior.log_prior_from_value`` which drops ``-log(b - a)``
124+
to return ``0.0``.
125+
126+
Out-of-support (``value`` outside ``[lower_limit, upper_limit]``) returns
127+
``-inf`` on the JAX path; the NumPy path returns ``-log(value)`` without
128+
bounds-checking to mirror the existing ``UniformPrior`` NumPy semantics
129+
which trust the search to stay inside bounds.
119130
120131
Parameters
121132
----------
122-
value : float
123-
The physical value of this prior's corresponding parameter in a `NonLinearSearch` sample.
133+
value
134+
The physical value of this prior's corresponding parameter in a
135+
``NonLinearSearch`` sample.
124136
xp
125137
Array-module to dispatch on (``numpy`` or ``jax.numpy``). Default ``numpy``.
126138
"""
127-
return 1.0 / value
139+
if xp is np:
140+
return -np.log(value)
141+
in_bounds = (value >= self.lower_limit) & (value <= self.upper_limit)
142+
return xp.where(in_bounds, -xp.log(value), -xp.inf)
128143

129144
def value_for(self, unit, xp=np):
130145
"""

autofit/messages/normal.py

Lines changed: 16 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -427,22 +427,33 @@ def value_for(self, unit, xp=np):
427427
inv = erfinv(2.0 * unit - 1.0)
428428
return self.mean + self.sigma * xp.sqrt(2.0) * inv
429429

430-
def log_prior_from_value(self, value: float, xp=np) -> float:
430+
def log_prior_from_value(self, value, xp=np):
431431
"""
432-
Compute the log prior probability of a given physical value under this Gaussian prior.
432+
Compute the log prior density of a given physical value under this Gaussian prior.
433433
434-
Used to convert a likelihood to a posterior in non-linear searches (e.g., Emcee).
434+
Returns ``log p(value)`` in density form: negative for values far from
435+
the mean, zero at the mode. ``Fitness._call`` adds this directly to the
436+
log-likelihood to form ``log_posterior = log_likelihood + sum(log_priors)``,
437+
which Emcee / Zeus / MLE-Drawer / LBFGS then maximise (LBFGS via
438+
``-2 * figure_of_merit`` minimisation).
439+
440+
The constant ``-log(sigma * sqrt(2*pi))`` is dropped (it is a true
441+
constant in the prior, irrelevant to posterior shape), matching the
442+
convention used by ``UniformPrior.log_prior_from_value`` which drops
443+
``-log(b - a)`` to return ``0.0``.
435444
436445
Parameters
437446
----------
438447
value
439448
A physical parameter value for which the log prior is evaluated.
449+
xp
450+
Array-module to dispatch on (``numpy`` or ``jax.numpy``). Default ``numpy``.
440451
441452
Returns
442453
-------
443-
The log prior probability of the given value.
454+
The log prior density at the given value, up to an additive constant.
444455
"""
445-
return (value - self.mean) ** 2.0 / (2 * self.sigma**2.0)
456+
return -((value - self.mean) ** 2.0) / (2 * self.sigma**2.0)
446457

447458
def __str__(self):
448459
"""

test_autofit/mapper/model/test_model_mapper.py

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -432,7 +432,11 @@ def test_log_prior_list_from_vector(self):
432432

433433
log_prior_list = mapper.log_prior_list_from_vector(vector=[0.0, 5.0])
434434

435-
assert log_prior_list == [0.125, 0.2]
435+
# Density-form log-priors (PyAutoLabs/PyAutoFit#1266):
436+
# Gaussian(mean=1, sigma=2) at value=0: -(0-1)**2 / (2*4) = -0.125
437+
# LogUniform(...) at value=5: -log(5)
438+
assert log_prior_list[0] == -0.125
439+
assert log_prior_list[1] == pytest.approx(-np.log(5.0), 1.0e-12)
436440

437441
def test_random_unit_vector_within_limits(self):
438442
mapper = af.ModelMapper()

test_autofit/mapper/prior/test_prior.py

Lines changed: 38 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
import math
22

3+
import numpy as np
34
import pytest
45

56
import autofit as af
@@ -190,33 +191,30 @@ def test__non_zero_lower_limit(self):
190191
assert log_uniform_half.value_for(0.5) == pytest.approx(0.70710678118, 1.0e-4)
191192

192193
def test__log_prior_from_value(self):
193-
gaussian_simple = af.LogUniformPrior(lower_limit=1e-8, upper_limit=1.0)
194-
195-
log_prior = gaussian_simple.log_prior_from_value(value=1.0)
196-
197-
assert log_prior == 1.0
198-
199-
log_prior = gaussian_simple.log_prior_from_value(value=2.0)
200-
201-
assert log_prior == 0.5
202-
203-
log_prior = gaussian_simple.log_prior_from_value(value=4.0)
204-
205-
assert log_prior == 0.25
206-
207-
gaussian_simple = af.LogUniformPrior(lower_limit=50.0, upper_limit=100.0)
208-
209-
log_prior = gaussian_simple.log_prior_from_value(value=1.0)
210-
211-
assert log_prior == 1.0
212-
213-
log_prior = gaussian_simple.log_prior_from_value(value=2.0)
214-
215-
assert log_prior == 0.5
194+
# LogUniformPrior log-density: -log(value), dropping the normalisation
195+
# constant -log(log(upper / lower)). Consistent with UniformPrior's
196+
# convention of returning 0.0 (dropping -log(b - a)).
197+
log_uniform = af.LogUniformPrior(lower_limit=1e-8, upper_limit=1.0)
198+
199+
assert log_uniform.log_prior_from_value(value=1.0) == 0.0
200+
assert log_uniform.log_prior_from_value(value=2.0) == pytest.approx(
201+
-np.log(2.0), 1.0e-12
202+
)
203+
assert log_uniform.log_prior_from_value(value=4.0) == pytest.approx(
204+
-np.log(4.0), 1.0e-12
205+
)
216206

217-
log_prior = gaussian_simple.log_prior_from_value(value=4.0)
207+
# The normalisation constant being dropped means the returned values
208+
# do NOT depend on the (lower_limit, upper_limit) pair — only on `value`.
209+
log_uniform = af.LogUniformPrior(lower_limit=50.0, upper_limit=100.0)
218210

219-
assert log_prior == 0.25
211+
assert log_uniform.log_prior_from_value(value=1.0) == 0.0
212+
assert log_uniform.log_prior_from_value(value=2.0) == pytest.approx(
213+
-np.log(2.0), 1.0e-12
214+
)
215+
assert log_uniform.log_prior_from_value(value=4.0) == pytest.approx(
216+
-np.log(4.0), 1.0e-12
217+
)
220218

221219
def test__lower_limit_zero_or_below_raises_error(self):
222220
with pytest.raises(exc.PriorException):
@@ -244,13 +242,16 @@ def test__non_zero_mean(self):
244242
@pytest.mark.parametrize(
245243
"mean, sigma, value, expected",
246244
[
245+
# Density-form log-prior: -(value - mean)**2 / (2 * sigma**2), with
246+
# the -log(sigma * sqrt(2 * pi)) normalisation constant dropped.
247+
# Maximum at value == mean (returns 0), negative elsewhere.
247248
(0.0, 1.0, 0.0, 0.0),
248-
(0.0, 1.0, 1.0, 0.5),
249-
(0.0, 1.0, 2.0, 2.0),
250-
(1.0, 2.0, 0.0, 0.125),
249+
(0.0, 1.0, 1.0, -0.5),
250+
(0.0, 1.0, 2.0, -2.0),
251+
(1.0, 2.0, 0.0, -0.125),
251252
(1.0, 2.0, 1.0, 0.0),
252-
(1.0, 2.0, 2.0, 0.125),
253-
(30.0, 60.0, 2.0, pytest.approx(0.108888, 1.0e-4)),
253+
(1.0, 2.0, 2.0, -0.125),
254+
(30.0, 60.0, 2.0, pytest.approx(-0.108888, 1.0e-4)),
254255
],
255256
)
256257
def test__log_prior_from_value(self, mean, sigma, value, expected):
@@ -265,4 +266,10 @@ def test_log_gaussian_prior_log_prior_from_value():
265266
)
266267

267268
assert log_gaussian_prior.log_prior_from_value(value=0.0) == float("-inf")
268-
assert log_gaussian_prior.log_prior_from_value(value=0.5) == 0.9333736875190459
269+
# Density form: -(log(value) - mean)**2 / (2 * sigma**2) - log(value),
270+
# where the second term is the Jacobian of the log-space transform.
271+
log_half = math.log(0.5)
272+
expected = -(log_half ** 2) / 2.0 - log_half
273+
assert log_gaussian_prior.log_prior_from_value(value=0.5) == pytest.approx(
274+
expected, 1.0e-12
275+
)

0 commit comments

Comments
 (0)