[Octave-bug-tracker] [bug #34362] [octave forge] (statistics) binopdf.m: Implement high accuracy Holder algorithm for n >= 10

[Nicholas Jankowski]

Follow-up Comment #6, bug #34362 (project octave):

looking at the paper octave's binopdf currently uses the form in (2), noted as
inaccurate for large n. mentions Matlab's binopdf function using it as well,
at least back in 2002. running a number of the examples in octave and matlab,
results do deviate 7+ sig figs out, so maybe they don't anymore.

anyway, the approach splits the calculation into 
p(x,n,p) = p(x,n,x/n)*exp(-D(x,n,p))

with the second part being called the Deviance term, more or less equal to 
D(x,n,p) = x ln (x/np) + (n-x) ln((n-x)/(n(1-p)))

and the first term rearranges to:
p(x,n,x/n) = sqrt(n/(2pi*x*(n-x))) * exp (d(n)-d(x)-d(n-x))
  where d(n) = ln(n!*exp(n)/ (sqrt(2pi*n)*n^n))

the R implementation of Loader's algo can be found at 
https://github.com/SurajGupta/r-source/blob/master/src/nmath/dbinom.c

with the deviance term sub function as:
https://github.com/SurajGupta/r-source/blob/master/src/nmath/bd0.c

will look and see how straight forward it is to bring into octave's binopdf

    _______________________________________________________

Reply to this item at:

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

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