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

[Rik]

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

I've never seen input values to nchoosek like 2^128.  I think that is well
outside the norm which is probably why this hasn't been noticed before.  I'm
not saying your application doesn't call for these numbers, just that
overwhelming number of use cases don't call for these values.

Generally nchoosek blows up very quickly such that it isn't conceivable to try
the direct approach of


log (nchoosek (n, k))


Something simple like


nchoosek (1e6,4)
warning: nchoosek: possible loss of precision
warning: called from
    nchoosek at line 151 column 7

ans = 4.1666e+22


is already complaining about lost precision.

Back to betaln:


lnb = gammaln (a) + gammaln (b) - gammaln (a + b)


The accuracy is limited not by the gammaln() function, but by the addition
operation using IEEE 8-byte doubles.  eps() is 2e-16 and if the difference in
scale between two arguments is greater than 1e16 then the smaller argument is
lost.


octave:19> format long
octave:20> 1 + 1e-16
ans = 1
octave:21> sprintf ('%.17d\n', 1 + 2e-16)
ans = 1.0000000000000002


I would go back to the base definition of nchoosek as


c = n! / k!(n-k)!


Take the logarithm of both sides and use the properties of logarithms to
write


log (c) = log (n!) - log (k!) - log ((n-k)!)


Then using the relationship


gamma (n + 1) = n!


I can write


log (c) = gammaln (n+1) - gammaln (k+1) - gammaln (n-k+1)


A quick anonymous function for this is


bignk = @(n,k) gammaln (n+1) - gammaln (k+1) - gammaln (n-k+1)


This still has a problem in that the expression "N-K+1" can't be calculated
without loss of precision with IEEE doubles.  So, I tried switching to
variable precision numbers which are part of the symbolic package.  You will
need to have python3-sympy installed on your computer.  After that, install
and load the symbolic package with


pkg install -forge symbolic
pkg load symbolic


Finally, create some numbers and try it out.


octave:47> n = vpa (2, 50)
n = (sym) 2.0000000000000000000000000000000000000000000000000
octave:48> n = n^128
n = (sym) 340282366920938463463374607431768211456.00000000000
octave:49> k = vpa (2, 50)
k = (sym) 2.0000000000000000000000000000000000000000000000000
octave:50> k = k^32
k = (sym) 4294967296.0000000000000000000000000000000000000000
octave:51> logc = bignk (n,k)
logc = (sym) 290091236578.66962544142734259366989135742187500000


Maybe this is what you are looking for?



    _______________________________________________________

Reply to this item at:

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

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