From eb7dd580bfd8a1a4bed29ce98b6bbeab7b8f5613 Mon Sep 17 00:00:00 2001 From: Remis Date: Wed, 2 Sep 2026 09:49:19 +0100 Subject: [PATCH] Add the absolute-increment rule for the Hessian ApproximateSD has always sized a single relative increment so that moving every parameter together changes the function by 0.5% of its own value. Colquhoun, Hatton & Hawkes (2003) J Physiol 547:699-728, p. 702, describe a different rule, and it is the one their Tables 2 and 3 were made with: find the increment for each parameter separately that changes the log likelihood by a specified absolute amount, and where no such increment exists, omit that parameter's row and column from the Hessian before inverting it. delta_L selects that rule. The per-parameter part is what matters: sizing one step by moving every parameter together moves along the diagonal ridge that an anisotropic likelihood has, where the function changes slowly, so the step that results is far too large for the directions across the ridge. kept and dropped say which parameters survived; hessian and covariance cover those only, while sd and correlations stay the full size of theta with nan in the dropped places, so sd[i] always refers to theta[i] and a table of results cannot silently misalign. A function that raises where it cannot be evaluated is treated as too large a step rather than propagating the error. The default path is untouched: the existing pinned Hessian, covariance, SD and correlation values are unchanged. 12 new tests, on quadratics whose Hessian is known exactly, so that a failure means the estimator is wrong rather than that a reference value drifted. One of them records that the flat-parameter case is precisely what the default rule cannot do. The suite goes from 36 to 48 passed. Co-Authored-By: Claude Opus 5 --- ekdist/errors.py | 181 +++++++++++++++++++++++++-- tests/test_errors_absolute_deltas.py | 138 ++++++++++++++++++++ 2 files changed, 312 insertions(+), 7 deletions(-) create mode 100644 tests/test_errors_absolute_deltas.py diff --git a/ekdist/errors.py b/ekdist/errors.py index 2156f08..e6b8092 100644 --- a/ekdist/errors.py +++ b/ekdist/errors.py @@ -74,16 +74,109 @@ def __lik_contour(self, x, num, theta, func, data): return -result.fun class ApproximateSD: - def __init__(self, theta, func, arg, delta_step=0.0001): - self.hessian = self.hessian_matrix(theta, func, arg, delta_step) + """Approximate standard deviations and correlations from the Hessian. + + The Hessian of the function being minimised is estimated by finite + differences, inverted to give the covariance matrix, and from that the + standard deviation of each estimate and the correlation between each pair. + ``func`` is the quantity that was *minimised* -- minus the log likelihood -- + so it increases as the parameters move away from the fitted values. + + Everything turns on the increments used for the finite differences. Two + rules are available. + + **The default**, unchanged since this class was written, sizes a single + relative increment so that moving *every* parameter together by that + relative amount changes ``func`` by 0.5% of its own value. + + **The absolute rule**, selected by giving ``delta_L``, is the one described + by Colquhoun, Hatton & Hawkes (2003), p. 702: the increment is found + separately **for each parameter**, and is the one that changes ``func`` by + ``delta_L`` -- an absolute amount, not a fraction of the value. Where no + such increment can be found, because the fit is insensitive to that + parameter, its row and column are omitted from the Hessian before it is + inverted, exactly as the paper describes, and its index is listed in + :attr:`dropped`. + + Which rule to use matters when the likelihood surface is strongly + anisotropic. Moving all the parameters together, as the default does, moves + along the diagonal ridge that such a surface has, where the function changes + slowly; the increment that results can be far too large for the directions + across the ridge. + + Parameters + ---------- + theta : ndarray + The fitted parameters. + func : callable + ``func(theta, arg) -> float``, the quantity minimised. + arg : object + Passed through to ``func``. + delta_step : float + Relative size of the first increment tried, before tuning. + delta_L : float, optional + Selects the absolute rule, and is the change in ``func`` each increment + is tuned to produce. Note that the paper's 0.1 "log units" are natural + logarithms; a likelihood in log10 wants 0.1/ln(10) = 0.0434. + max_relative_delta : float + An increment larger than this multiple of its own parameter counts as + not found, and the parameter is dropped. Absolute rule only. + + Attributes + ---------- + hessian, covariance : ndarray + Square, over the parameters actually used -- smaller than ``theta`` if + any were dropped. + kept, dropped : tuple of int + Indices into ``theta``. + sd, correlations : ndarray + Full size, indexed by parameter, with ``nan`` for dropped parameters, + so that ``sd[i]`` always refers to ``theta[i]``. + deltas : ndarray + The increments used; ``nan`` for dropped parameters. + """ + + def __init__(self, theta, func, arg, delta_step=0.0001, + delta_L=None, max_relative_delta=100.0): + theta = np.asarray(theta, float) + if delta_L is None: + self.deltas = self.__optimal_deltas(theta, func, arg, delta_step) + self.dropped = () + else: + self.deltas = self.__absolute_deltas( + theta, func, arg, delta_step, delta_L, max_relative_delta) + self.dropped = tuple( + int(i) for i in np.flatnonzero(~np.isfinite(self.deltas))) + self.kept = tuple(i for i in range(theta.size) if i not in self.dropped) + + self.hessian = self.hessian_matrix(theta, func, arg, delta_step, + deltas=self.deltas, keep=self.kept) self.covariance = nplin.inv(self.hessian) - self.sd = np.sqrt(self.covariance.diagonal()) - self.correlations = self.correlation_matrix(self.covariance) - def hessian_matrix(self, theta, LLfunc, args, delta_step=0.0001): - """ """ + self.sd = np.full(theta.size, np.nan) + self.sd[list(self.kept)] = np.sqrt(self.covariance.diagonal()) + reduced = self.correlation_matrix(self.covariance) + self.correlations = np.full((theta.size, theta.size), np.nan) + for a, i in enumerate(self.kept): + for b, j in enumerate(self.kept): + self.correlations[i, j] = reduced[a, b] + + def hessian_matrix(self, theta, LLfunc, args, delta_step=0.0001, + deltas=None, keep=None): + """Second derivatives of ``LLfunc`` by central differences. + + ``deltas`` and ``keep`` are supplied by :meth:`__init__`; called + directly with neither, this behaves as it always has. + """ + theta = np.asarray(theta, float) + if deltas is None: + deltas = self.__optimal_deltas(theta, LLfunc, args, delta_step) + if keep is None: + keep = tuple(range(theta.size)) + if len(keep) != theta.size: + return self.__hessian_subset(theta, LLfunc, args, deltas, keep) + hess = np.zeros((theta.size, theta.size)) - deltas = self.__optimal_deltas(theta, LLfunc, args, delta_step) # Diagonal elements of Hessian coe11 = np.array([theta.copy(), ] * theta.size) + np.diag(deltas) coe33 = np.array([theta.copy(), ] * theta.size) - np.diag(deltas) @@ -112,6 +205,80 @@ def hessian_matrix(self, theta, LLfunc, args, delta_step=0.0001): (4 * deltas[i] * deltas[j])) return hess + def __hessian_subset(self, theta, LLfunc, args, deltas, keep): + """The Hessian over ``keep`` only, the rest omitted as in the paper.""" + keep = list(keep) + n = len(keep) + hess = np.zeros((n, n)) + L0 = LLfunc(theta, args) + for a, i in enumerate(keep): + up, down = theta.copy(), theta.copy() + up[i] += deltas[i] + down[i] -= deltas[i] + hess[a, a] = ((LLfunc(up, args) - 2.0 * L0 + LLfunc(down, args)) + / deltas[i] ** 2) + for a, i in enumerate(keep): + for b, j in enumerate(keep): + if a >= b: + continue + pp, pm, mp, mm = (theta.copy(), theta.copy(), + theta.copy(), theta.copy()) + pp[i] += deltas[i]; pp[j] += deltas[j] + pm[i] += deltas[i]; pm[j] -= deltas[j] + mp[i] -= deltas[i]; mp[j] += deltas[j] + mm[i] -= deltas[i]; mm[j] -= deltas[j] + value = ((LLfunc(pp, args) - LLfunc(pm, args) + - LLfunc(mp, args) + LLfunc(mm, args)) + / (4 * deltas[i] * deltas[j])) + hess[a, b] = hess[b, a] = value + return hess + + def __change_in(self, theta, func, args, index, delta): + """How much ``func`` rises when parameter ``index`` moves by ``delta``. + + ``inf`` if it cannot be evaluated there, which counts as too large a + step and makes the caller try a smaller one. + """ + trial = theta.copy() + trial[index] += delta + try: + value = func(trial, args) + except (ArithmeticError, ValueError, RuntimeError): + return np.inf + if not np.isfinite(value): + return np.inf + return value - self.__L0 + + def __absolute_deltas(self, theta, func, args, step_factor, delta_L, + max_relative): + """One increment per parameter, each changing ``func`` by ``delta_L``. + + The rule of Colquhoun, Hatton & Hawkes (2003), p. 702. Returns ``nan`` + for a parameter whose increment could not be found, which is how a + parameter the fit is insensitive to gets dropped. + """ + self.__L0 = func(theta, args) + deltas = np.full(theta.size, np.nan) + for i in range(theta.size): + scale = abs(theta[i]) if theta[i] != 0.0 else 1.0 + biggest = max_relative * scale + delta = step_factor * scale + change = self.__change_in(theta, func, args, i, delta) + if change < delta_L: # too small: grow it + while change < delta_L and delta < biggest: + delta *= 2.0 + change = self.__change_in(theta, func, args, i, delta) + else: # too big: shrink it + smallest = 1e-12 * scale + while change > delta_L and delta > smallest: + delta /= 2.0 + change = self.__change_in(theta, func, args, i, delta) + delta *= 2.0 # the last one that was + change = self.__change_in(theta, func, args, i, delta) + if np.isfinite(change) and change >= delta_L and delta <= biggest: + deltas[i] = delta + return deltas + def __tune_deltas(self, theta, func, args, Lcrit, deltas, increase=True): factor = [1, 2] if increase else [-1, 0.5] count = 0 diff --git a/tests/test_errors_absolute_deltas.py b/tests/test_errors_absolute_deltas.py new file mode 100644 index 0000000..368d5ce --- /dev/null +++ b/tests/test_errors_absolute_deltas.py @@ -0,0 +1,138 @@ +"""The absolute-increment rule for the Hessian. + +``ApproximateSD`` has always sized one relative increment so that moving every +parameter together changes the function by 0.5% of its value. Colquhoun, Hatton +& Hawkes (2003), p. 702, describe a different rule: find the increment for +**each** parameter separately that changes the log likelihood by a specified +absolute amount, and where no such increment exists, omit that parameter's row +and column from the Hessian. ``delta_L`` selects that rule. + +The tests use quadratics, whose Hessian is known exactly, so that a failure +means the estimator is wrong rather than that the reference values drifted. +""" + +import numpy as np +import numpy.testing as npt +import pytest + +from ekdist import errors + + +def quadratic(H, centre, offset=100.0): + """A function with Hessian ``H``, minimised at ``centre``.""" + H = np.asarray(H, float) + centre = np.asarray(centre, float) + + def f(theta, args=None): + d = np.asarray(theta, float) - centre + return offset + 0.5 * d @ H @ d + + return f + + +WELL_CONDITIONED = np.array([[4.0, 1.0], [1.0, 9.0]]) +CENTRE = np.array([2.0, 3.0]) + + +@pytest.mark.parametrize('kwargs', [{}, {'delta_L': 0.1}, {'delta_L': 2.0}]) +def test_both_rules_recover_a_known_hessian(kwargs): + asd = errors.ApproximateSD(CENTRE, quadratic(WELL_CONDITIONED, CENTRE), + None, **kwargs) + npt.assert_allclose(asd.hessian, WELL_CONDITIONED, rtol=1e-6) + npt.assert_allclose( + asd.sd, np.sqrt(np.diag(np.linalg.inv(WELL_CONDITIONED))), rtol=1e-6) + assert asd.dropped == () + + +def test_absolute_rule_gives_the_right_correlation(): + covariance = np.linalg.inv(WELL_CONDITIONED) + expected = covariance[0, 1] / np.sqrt(covariance[0, 0] * covariance[1, 1]) + asd = errors.ApproximateSD(CENTRE, quadratic(WELL_CONDITIONED, CENTRE), + None, delta_L=0.1) + npt.assert_allclose(asd.correlations[0, 1], expected, rtol=1e-6) + npt.assert_allclose(np.diag(asd.correlations), [1.0, 1.0], rtol=1e-9) + + +def test_each_increment_produces_the_requested_change(): + """That is what 'the increment needed to change L by delta_L' means.""" + delta_L = 0.1 + f = quadratic(WELL_CONDITIONED, CENTRE) + asd = errors.ApproximateSD(CENTRE, f, None, delta_L=delta_L) + L0 = f(CENTRE) + for i, delta in enumerate(asd.deltas): + moved = CENTRE.copy() + moved[i] += delta + change = f(moved) - L0 + # the search doubles or halves, so it lands within a factor of two + assert delta_L <= change < 4 * delta_L, (i, change) + + +def test_increments_differ_between_parameters(): + """The point of the rule: a stiff parameter gets a smaller increment.""" + H = np.array([[1000.0, 0.0], [0.0, 1.0]]) + asd = errors.ApproximateSD(CENTRE, quadratic(H, CENTRE), None, delta_L=0.1) + assert asd.deltas[0] < asd.deltas[1] / 10 + + +# ------------------------------------------------- insensitive parameters + +FLAT = np.array([[4.0, 1.0, 0.0], [1.0, 9.0, 0.0], [0.0, 0.0, 1e-14]]) +CENTRE3 = np.array([2.0, 3.0, 5.0]) + + +def test_insensitive_parameter_is_dropped(): + asd = errors.ApproximateSD(CENTRE3, quadratic(FLAT, CENTRE3), None, + delta_L=0.1) + assert asd.dropped == (2,) + assert asd.kept == (0, 1) + assert asd.hessian.shape == (2, 2) + assert all(isinstance(i, int) for i in asd.dropped) + + +def test_dropping_leaves_the_others_exact(): + """The point of omitting a row and column rather than giving up.""" + asd = errors.ApproximateSD(CENTRE3, quadratic(FLAT, CENTRE3), None, + delta_L=0.1) + npt.assert_allclose(asd.hessian, FLAT[:2, :2], rtol=1e-6) + covariance = np.linalg.inv(FLAT[:2, :2]) + npt.assert_allclose(asd.sd[:2], np.sqrt(np.diag(covariance)), rtol=1e-6) + + +def test_dropped_entries_are_nan_and_indexed_by_parameter(): + """sd[i] must still refer to theta[i], or a table of results misaligns.""" + asd = errors.ApproximateSD(CENTRE3, quadratic(FLAT, CENTRE3), None, + delta_L=0.1) + assert asd.sd.shape == (3,) + assert asd.correlations.shape == (3, 3) + assert np.isnan(asd.sd[2]) and np.isnan(asd.deltas[2]) + assert np.all(np.isnan(asd.correlations[2, :])) + assert np.all(np.isnan(asd.correlations[:, 2])) + assert np.all(np.isfinite(asd.correlations[:2, :2])) + + +def test_the_flat_case_is_what_the_default_rule_cannot_do(): + """Motivation for the option, kept as a test so it stays true.""" + with pytest.raises(np.linalg.LinAlgError): + errors.ApproximateSD(CENTRE3, quadratic(FLAT, CENTRE3), None) + + +def test_a_function_that_raises_is_treated_as_too_large_a_step(): + """Our HJC likelihood raises on a Q matrix it cannot handle.""" + f = quadratic(WELL_CONDITIONED, CENTRE) + + def fragile(theta, args=None): + if abs(theta[0] - CENTRE[0]) > 0.05: + raise ArithmeticError('too far') + return f(theta, args) + + asd = errors.ApproximateSD(CENTRE, fragile, None, delta_L=1e-4) + assert asd.deltas[0] <= 0.05 + assert np.isfinite(asd.sd[0]) + + +def test_zero_valued_parameter_does_not_give_a_zero_increment(): + centre = np.array([0.0, 3.0]) + asd = errors.ApproximateSD(centre, quadratic(WELL_CONDITIONED, centre), + None, delta_L=0.1) + assert asd.deltas[0] > 0 + npt.assert_allclose(asd.hessian, WELL_CONDITIONED, rtol=1e-6)