Skip to content
Snippets Groups Projects

Improve numerical stability of conditional_ppf

Merged Improve numerical stability of conditional_ppf
leo-singer/ligo.skymap:marginal-ppf-reparametrize into master
Merged Leo Pound Singer requested to merge leo-singer/ligo.skymap:marginal-ppf-reparametrize into master
  1. Aug 18, 2020
    • Leo Pound Singer's avatar
      Improve initial guess for conditional_cdf · b0476ebb
      Leo Pound Singer authored
      b0476ebb
    • Leo Pound Singer's avatar
      Improve numerical stability of conditional_ppf · d64bc06c
      Leo Pound Singer authored 
      Improve the numerical stability of the method
`ligo.skymap.distance.conditional_ppf` by reparametrizing the equation
that is being solved. This method, which calculates the inverse of the
distance CDF, works by solving the equation `f(x) - p = 0` for `x`,
where `f(x)` is the distance CDF, and `p` is the desired probability.

The reparametrized equation is `log(1 - f(x)) - log(1 - p) = 0` if `p >
1/2` and `log(f(x)) - log(p) = 0` otherwise. This reparametrization is
effective because it improves the dynamic range in the tails of the
distribution. This same reparametrization had already proven effective
in the related method `ligo.skymap.distance.marginal_ppf`.

This change also fixes some rare corner cases where `marginal_ppf`
returned silly values becauses it uses `conditional_ppf` internally to
create its own initial guess. One example was the median distance for
the binary neutron star candidate S191205ah. Before this patch, the
result was negative and invalid:

    >>> from ligo.skymap.distance import marginal_ppf
    >>> from ligo.skymap.moc import uniq2pixarea
    >>> from ligo.skymap.io import read_sky_map
    >>> url = 'https://gracedb.ligo.org/apiweb/superevents/S191205ah/files/bayestar.multiorder.fits'
    >>> s = read_sky_map(url, moc=True)
    >>> marginal_ppf(0.5, s['PROBDENSITY'] * uniq2pixarea(s['UNIQ']),
    ...              s['DISTMU'], s['DISTSIGMA'], s['DISTNORM'])
    /Users/lpsinger/src/ligo.skymap/ligo/skymap/util/numpy.py:46:
    RuntimeWarning: invalid value encountered in marginal_ppf
      return func(*args, **kwargs)
    -223357.8508233767

After this patch, the result is positive and sensible:

    >>> marginal_ppf(0.5, s['PROBDENSITY'] * uniq2pixarea(s['UNIQ']),
    ...              s['DISTMU'], s['DISTSIGMA'], s['DISTNORM'])
    362.7485740018039
      d64bc06c