beta function: when b<epsilon return tgamma(b)?

[Eric Butter <egbutter@…>]

In the beta.hpp header for the beta function (i.e., *not* the beta distribution) around line 131, you have:

if((c == a) && (b < tools::epsilon<T>()))

return boost::math::tgamma(b, pol);

else if((c == b) && (a < tools::epsilon<T>()))

return boost::math::tgamma(a, pol);

Don't you actually want if b<epsilon, tgamma(a)? beta(a,b) reduces to (approximately) gamma(a)b^{-a for a>>b, which seems like the case for which you are trying to account?}

Cheers,

Eric

[John Maddock]

I believe the code to be correct:

if a >> b such that a + b == a (at the precision we're working to), then:

gamma(a) / gamma(a+b) == 1

so:

beta(a, b) = gamma(a) * gamma(b) / gamma(a+b) ~ gamma(b) ~ 1/b

[Eric Butter <egbutter@…>]

Ah yes that makes sense. I definitely misinterpreted that identity. Thanks for setting me straight!