[Octave-bug-tracker] [bug #65230] Poor accuracy of betaln for high ratios of a/b

[anonymous]

Follow-up Comment #2, bug#65230 (group octave):

OP here.  Thank you for confirming, and for reporting the behavior of Matlab. 
I see that the Matlab documentation also doesn't explicitly mention accuracy,
but the "Algorithms" section does seem to indicate that the same
implementation based on gammaln is used.

I don't really know what the logarithm of the beta function is commonly used
for by others, but betaln has returned 0 for every real case I've tried.  For
example, I needed to compute log(nchoosek(n, k)) where (n, k) was something
like (2^128, 2^30) or (2^256, 2^40).  I tried

function [retval] = log_nchoosek (n, k)
  retval = -log1p(n) - betaln(k + 1, n - k + 1);
end

which doesn't work for the sorts of values I was interested in.  Replacing
betaln with gsl_sf_lnbeta seemed to do the job.

Regarding this being a corner case, those are mighty big corners!  Of the
input space where both arguments are greater than 1, betaln seems to return 0
for something like 88%, NaN for another 2%, and arguably has poor accuracy for
most of the rest.  Still, I suppose I know what you mean.  This issue appears
to have existed for over 17 years without much complaint.

If the current implementation is retained, I think it would be helpful for the
documentation to disclose the implementation and mention the accuracy
implications.



    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?65230>

_______________________________________________
Message sent via Savannah
https://savannah.gnu.org/