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)