Improve numerical stability of conditional_ppf

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

src/bayestar_distance.c

@@ -78,13 +78,24 @@ double bayestar_distance_conditional_cdf(
 }
 
 
+typedef struct {
+    double p, mu, norm;
+} conditional_ppf_params;
+
+
 static void conditional_ppf_fdf(double r, void *params, double *f, double *df)
 {
-    const double p = ((double*) params)[0];
-    const double mu = ((double *) params)[1];
-    const double norm = ((double *) params)[2];
-    *f = bayestar_distance_conditional_cdf(r, mu, 1, norm) - p;
-    *df = bayestar_distance_conditional_pdf(r, mu, 1, norm);
+    const conditional_ppf_params *p = (conditional_ppf_params *)params;
+    const double _f = bayestar_distance_conditional_cdf(r, p->mu, 1, p->norm);
+    const double _df = bayestar_distance_conditional_pdf(r, p->mu, 1, p->norm);
+    if (p->p > 0.5)
+    {
+        *f = log(1 - _f) - log(1 - p->p);
+        *df = -_df / (1 - _f);
+    } else {
+        *f = log(_f) - log(p->p);
+        *df = _df / _f;
+    }
 }
 
 
@@ -104,6 +115,22 @@ static double conditional_ppf_df(double r, void *params)
 }
 
 
+static double conditional_ppf_initial_guess(double p, double mu)
+{
+    /* Initial guess: ignore r^2 term;
+     * distribution becomes truncated Gaussian */
+    const double lo = gsl_cdf_ugaussian_P(-mu);
+    const double z = gsl_cdf_ugaussian_Pinv(p + (1 - p) * gsl_cdf_ugaussian_P(-mu)) + mu;
+
+    if (z > 0)
+        return z;
+    else if (mu > 0)
+        return mu;  /* Fallback 1: mean */
+    else
+        return 0.5;  /* Fallback 2: constant value */
+}
+
+
 double bayestar_distance_conditional_ppf(
     double p, double mu, double sigma, double norm)
 {
@@ -121,9 +148,9 @@ double bayestar_distance_conditional_ppf(
 
     /* Set up variables for tracking progress toward the solution. */
     static const int max_iter = 50;
-    double params[] = {p, mu, norm};
+    conditional_ppf_params params = {p, mu, norm};
     int iter = 0;
-    double z = mu > 0 ? mu : 0.5; /* FIXME: better initial guess? */
+    double z = conditional_ppf_initial_guess(p, mu);
     int status;
 
     /* Set up solver (on stack). */