[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Octave-bug-tracker] [bug #34362] [octave forge] (statistics) binopdf.m:
From: |
Nicholas Jankowski |
Subject: |
[Octave-bug-tracker] [bug #34362] [octave forge] (statistics) binopdf.m: Implement high accuracy Holder algorithm for n >= 10 |
Date: |
Fri, 3 Dec 2021 15:27:15 -0500 (EST) |
User-agent: |
Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:78.0) Gecko/20100101 Firefox/78.0 |
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/