Quantile for skew_normal returns invalid value

[yurybura]

Code example:

#include <boost/math/distributions/skew_normal.hpp>
#include <cstdlib>

int main()
{
  boost::math::skew_normal skew(573.39724735636185, 77.0, 4.0);
  const double x = boost::math::quantile(skew, 0.00285612015554148);
  const double y = boost::math::quantile(skew, 0.00285612015554149);
  const double z = boost::math::quantile(skew, 0.00285612015554150);
  std::printf("x=%g\n", x);
  std::printf("y=%g\n", y);
  std::printf("z=%g\n", z);
}

Result on Windows 10 + MSVC 19.16.27026.1 + Boost 1.67.0:

x=539.858
y=-inf
z=539.858

Result on Windows 10 + MSVC 19.16.27026.1 + Boost 1.69.0:

x=539.858
y=0
z=539.858

[NAThompson]

I just ran your code on develop, with both g++-8 and Apple clang. I get:

x=539.858
y=539.858
z=539.858

Could you clone the develop and see if it's still there in your environment?

[jzmaddock]

I can reproduce, it's a numerical instability that triggers a bug in the Newton-Raphson code, I'm just not sure at present if this is a one off or part of a larger issue...

[yurybura]

Test project: https://github.com/yurybura/test1
Build log: https://ci.appveyor.com/project/yurybura/test1

[pabristow]

I can also confirm and that it goes wrong before or here
delta becomes inf while delta1 and delta2 is 1.344e-14

 if(fabs(delta * 2) > fabs(delta2))
  {
     // last two steps haven't converged.
     delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;  <<<<<
     if (fabs(delta) > fabs(result))
        delta = sign(delta) * result; // protect against huge jumps!
     // reset delta2 so we don't take this branch next time round:
     delta1 = delta;
     delta2 = 3 * delta;
  }

at iteration down to 192 (from a very large max_iterations estimate of 200 - is this right?

but I am not expert enough to see why and how to correct it.

[jzmaddock]

Fixed in develop in the commit referenced above.

[yurybura]

Issue is returned in Boost 1.84.0.
The example returns on MSVC 19.38.33134.0 + Boost 1.84.0:

x=539.858
y=-1.43195e+250
z=539.858

[jzmaddock]

Confirmed, though strangely the cause is different this time.

Things to check and improve:

1. The derivative is going to zero - why?
2. The zero-derivative handling code in newton_raphson_iterate needs to be updated along the same lines as the previous fix.
3. The quantile function has no error handling for exhausting the max root iterations.
4. We can probably narrow the search bounds anyway before iterating.

[jzmaddock]

OK, ignore (1) and (2), this got reverted as part of #1002 which was evidently a bad idea. (3) and (4) we should still look at.

[jzmaddock]

This should now be fixed again, with better testing this time!