Prevent passing denormals in calculation.

[jzmaddock]

Refs scipy/scipy#20693

[codecov]

Codecov Report

Attention: Patch coverage is 99.32886% with 1 lines in your changes are missing coverage. Please review.

Project coverage is 93.71%. Comparing base (cf0d343) to head (20f44d1).
Report is 4 commits behind head on develop.

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #1141      +/-   ##
===========================================
+ Coverage    93.69%   93.71%   +0.02%     
===========================================
  Files          772      774       +2     
  Lines        61168    61271     +103     
===========================================
+ Hits         57311    57422     +111     
+ Misses        3857     3849       -8     

Files	Coverage Δ
include/boost/math/concepts/std_real_concept.hpp	100.00% <ø> (ø)
...t/math/distributions/detail/hypergeometric_pdf.hpp	96.62% <100.00%> (ø)
...h/distributions/detail/hypergeometric_quantile.hpp	96.84% <100.00%> (ø)
...lude/boost/math/distributions/non_central_beta.hpp	91.96% <100.00%> (+0.16%)	⬆️
include/boost/math/distributions/non_central_t.hpp	97.62% <100.00%> (-0.03%)	⬇️
...e/boost/math/quadrature/detail/exp_sinh_detail.hpp	96.71% <100.00%> (+1.47%)	⬆️
...clude/boost/math/special_functions/jacobi_zeta.hpp	100.00% <ø> (ø)
test/exp_sinh_quadrature_test.cpp	99.13% <100.00%> (+0.05%)	⬆️
test/nc_t_pdf_data.ipp	100.00% <ø> (ø)
test/tanh_sinh_quadrature_test.cpp	100.00% <100.00%> (ø)
... and 4 more

... and 3 files with indirect coverage changes

---

Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update cf0d343...20f44d1. Read the comment docs.

[dschmitz89]

@jzmaddock : this looks impressive! One question for the addition to scipy. Do you have an external source for reference values for the noncentral T PDF? I did not check the complete (large) testing diff, might have missed it.

If not, I would implement the formula currently used in scipy with the high precision mpmath python library.

[dschmitz89]

Ah, this might be arbitrary precision implementation boost style that generates the test values after all?

[jzmaddock]

Nod. I used the hypergeometric formula and multiprecision arithmetic for the new test values.

One thing I should have mentioned: the integration method though accurate, is about 10x slower than the alternatives :( Other than truncation to zero, there really is no alternative though... as I'm typing though I'm just wondering if wolframalpha would generate a series for the integral...

[ckormanyos]

just wondering if wolframalpha would generate a series

I'm not all that familiar with this function but I've recently been doing a lot of computer algebra for Series and Pade approximants in other projects. Does anyone have a link to representation of the integral mentioned? I coould try for an approximation.

[jzmaddock]

I'm not all that familiar with this function but I've recently been doing a lot of computer algebra for Series and Pade approximants in other projects. Does anyone have a link to representation of the integral mentioned? I coould try for an approximation.

Cool, it's:

integral y^v exp(-1/2 ((y - mu * x / sqrt(x * x + v)))^2) dy for x = 0 to infinity

Which is the integral part of this formula https://wikimedia.org/api/rest_v1/media/math/render/svg/0059c7e6e3afb0e4496fecdcb3b85296c162f065

Ideally we need this whenever y < 0 in which domain the result is expected to be tiny.

[ckormanyos]

try for an approximation

The first thing I tried is the exponent part alone without integrating over $dy$, wondering if series expansion and subsequent integration would work. And it is normalized by dividing by the large exponent term at the end of the expression. If the series is continued in $x$, then integration of $y^{\nu}$ should be straightforward enough. It's been a long day, but @jzmaddock your intuition is correct, it seems like a series could do the trick.

W-Alpha didn't do it. I needed a bit more power and I ended up doing this first try locally.

Normal[Series[Exp[-1/2 ((y - mu*x/Sqrt[x*x + v]))^2], {x, Infinity, 4}] / Exp[-(1/2) (y - mu)^2]]

Note that we use Normal[] to truncate the series.

Here is a picture of the input and its output.

[ckormanyos]

OK I included the scaling by the exponential term added the multiplication by $y^{\nu}$ in the integrand and did the definite integral. It gives results like the following to Order $4$. See below.

Is there A numerical example I could use to check the approximation to Order $4$

[ckormanyos]

Actually, @jzmaddock I don't know if expanding at $x{\rightarrow}{\infty}$ was correct. Are the presumptions in my approximation above even correct?

[jzmaddock]

Actually, @jzmaddock I don't know if expanding at was correct. Are the presumptions in my approximation above even correct?

Since we're really interested in x < 0, I would say expansion either at x == 0 or x = -infinity would be better?

I also note that the curve being integrated, is extremely "pointy" - all the area is located in a tiny area - so this might be getting out of hand, but a series around the maxima (for y) would be optimal I guess?

If we fix the SciPy test case, with x = -1, v = 8 and mu=16 then https://www.wolframalpha.com/input?i=x%5E8+exp%28%28%28x+-+16+*+-1+%2F+sqrt%281+%2B+8%29%29%29%5E2+%2F+-2%29%3B shows the curve nicely, most of the peak in this case is outside of [0, INF] but that's not always the case. Interestingly, with all the variables fixed like this, wolfram alpha gives a simple(-ish) closed form for the integral, and for the summit of the curve. Does this help at all? Possibly not, since evaluating the definite integral in this case is the difference between two large values?

The best way to sanity check the series would be against the integral itself:

            boost::math::quadrature::exp_sinh<T, Policy> integrator;
            // Remove this line to check just the integral part:
            T integral = pow(v, v / 2) * exp(-v * mu * mu / (2 * (x * x + v)));
            if (integral != 0)
            {
               integral *= integrator.integrate([&x, v, mu](T y) 
                  {
                     T p;
                     if (v * log(y) < tools::log_max_value<T>())
                        p = pow(y, v) * exp(boost::math::pow<2>((y - mu * x / sqrt(x * x + v))) / -2);
                     else
                        p = exp(log(y) * v + boost::math::pow<2>((y - mu * x / sqrt(x * x + v))) / -2);
                     return p; 
                  });
            }

Extracted from non_central_t.hpp.

[ckormanyos]

Hi John (@jzmaddock) OK I am getting closer.

I just tried:

Integrate[(y^nu) Exp[-((y - ((mu x)/Sqrt[x x + nu]))^2)/2], {y, 0, Infinity}]

and received a closed form answer.

[jzmaddock]

There seems to be something wrong there if the answer depends only on x^2 ?

[jzmaddock]

I think it should have been Integrate[ [x^v exp(-1/2 ((x - mu * t / sqrt(t^2 + v)))^2)], {x, 0, Infinity}] but wolframalpha chokes on that.

[ckormanyos]

think it should have been

Geez this is challenging. Like this then...

Integrate[x^v Exp[-1/2 ((x - mu t/Sqrt[t^2 + v]))^2], {x, 0, Infinity}]

Whith the following answer.

But I think there are some very large terms and this is numerically unstable. I haven't quite figured out what to do with that yet. I'll try off and onn. Does not seem to be critical. If I ever come up with something, I'll report.

[ckormanyos]

Then I evaluated the numerical point at the test case, but internally the computer algebra system enlarges the precision. So the cancellations do not influence the result. But I haven't figured out what to do at fixed precison.

[ckormanyos]

It's a nightmare. Here are the two added terms in the result. The cancellations quench the result for machine precision.

I wish I knew more about asymptotics and perturbation theory... Arrrggghhh

[jzmaddock]

Well the good news is that 2.29e-9 is the correct answer ;)

But yes, you've basically re-invented the hypergeometric approximation, and the cancellation is indeed terrible :(

We might well be out of ideas at this point, but I appreciate the effort!

[ckormanyos]

might well be out of ideas at this point

I tried expanding in parameters - the first trick of asymptotics. Expanding in ${\mu}$ worked but still suffered from severe cancellations. Expansion in ${\nu}$ did not converge.

You know, this would be a great spot for the infamous double-double? I don't know if you'd like to throw Multiprecision at it, but if you double the working precision, you get, as you know, the digits back. But Multiprecision doesn't quite work with standalone Math. Somewhere down the line, Math could benefit from its own kind of a double-double.

Other than that, I'm out of ideas also.

[jzmaddock]

The issue with the hypergeomentrics is that when v is small, the largest term in the hypergeometric series is of the order (mu^2)! which will quite rapidly overflow. And even if it doesn't you end up subtracting two arbitrarily large numbers :( So I don't think even a double-double would do it, better to suck it up and use numerical integration I would guess.

[ckormanyos]

The issue with the hypergeomentrics is that when v is small, the largest term in the hypergeometric series is of the order (mu^2)! which will quite rapidly overflow. And even if it doesn't you end up subtracting two arbitrarily large numbers :( So I don't think even a double-double would do it, ...

Indeed. I even tried expansion in ${\mu}$ and although the expansion was promising, the cancellation was also present.

better to suck it up and use numerical integration I would guess.

Yes. Thanks John. I'll hang up my hat on this one now.