Add the absolute-increment rule for the Hessian

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

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)