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

[Nicholas Jankowski]

Update of bug #34362 (project octave):

                  Status:               Postponed => Patch Submitted        
                 Summary: [octave forge] (statistics) binopdf.m: Implement
high accuracy Holder algorithm for n >= 10 => [octave forge] (statistics)
binopdf.m: Implement high accuracy Loader algorithm for n >= 10

    _______________________________________________________

Follow-up Comment #8:

ok, that was a fun exercise. Attached is a patch and m-file for a modified
binopdf using the Loader algorithm.  It passes all tests, and I added a few
more checked against Matlab 2021b.  

1) new version does run about 2-3x slower than the current one, which is a
much simpler function just calling gammaln expansions. a bit of profiler work
found some things, but there's not enough detail in octave's profiler to find
specifics. (both versions spend 10% of their time on common_size, so that's
something that maybe could use a visit).  So that could use some work but I
didn't see obvious improvement areas. ( binopdf([1:50], 50, 0.5) is a decent
benchmark).

2) checking output accuracy of current version and this new octave version vs
matlab, R-implementation vs 'exact' VPA version for some of the help values
and high-n problem values identified in the Loader paper.  clipped off the
values after the first significant error for easier visual

x = 10, n = 20, p = 0.5

source       acc. digits  calculated value
Octave-vpa   --           0.176197052001953125
Octave-old   16           0.1761970520019532
Octave-new   exact        0.176197052001953125
Matlab2021b  exact        0.176197052001953125
R            exact        0.176197052001953125


x = 32, n = 65, p = 0.5

source       acc. digits  calculated value
Octave-vpa   --          
0.097841499903300731429394571092572618908889126032648 
Octave-old   13           0.09784149990329
Octave-new   16           0.09784149990330075
Matlab2021b  17           0.097841499903300738
R            16           0.09784149990330071


x = 0, n = 200, p = 0.02

source       acc. digits  calculated value
Octave-vpa   --          
0.017587946605721564931523103244847529749829653483380
Octave-old   16           0.01758794660572150
Octave-new   16           0.01758794660572150
Matlab2021b  16           0.01758794660572150
R            17           0.017587946605721566


x = 1e6, n = 2e6, p = 0.5

source       acc. digits  calculated value
Octave-vpa   --          
0.00056418951302406275121241988148617714907257387640426
Octave-old    9           0.000564189514
Octave-new   15           0.000564189513024063
Matlab2021b  17           0.00056418951302406278
R            16           0.0005641895130240628

x=0.3*n, n = 10^m, p = 0.3
 m Octave-old Octave-new Matlab2021b    R    Octave-vpa
 1         16         17          17   17   
2.6682793199999999999999999999999999999999999999999e-01
 2         14         15          16   17   
8.6783864753428087578960637912665039399024083554448e-02
 3         13         16          17   16   
2.7521003821268385527122960139188908419555755870387e-02
 4         11         17          17   17   
8.7053613650670565189701046818056804484495789345259e-03
 5         11         16          18   16   
2.7529546483974279546245692773645461883159351045387e-03
 6         10         15          16   16   
8.7056315463668078155717900870045623281585477425211e-04
 7          9         17          16   17   
2.7529631924020770334104364013945090603832813085629e-04
 8          7         15          16   15   
8.7056342482221605820748456181145242963308071359854e-05
 9          7         15          17   15   
2.7529632778422574377373933732840027323005242648868e-05
10          6         18          18   18   
8.7056342752407183443770621314073192664403859245026e-06
11          5         16          17   16   
2.7529632786966592551717359812182050924881279252233e-06
12          3         15          17   15   
8.7056342755109039224235471269393865017264362335144e-07




So the low value accuracy is comparable, and high order is much improved,
especially as seen in the last set. 

3) precision does lag a digit or two behind matlab and R. I don't see anything
obvious pointing to the issue, but I didn't extensively test after tackling
major numerical error spots to get it to this point. If an expert in such
things wants to peek and tweak, please do.

4) I did make one arbitrary change - old version errored on any complex input.
 matlab allows you to pass complex inputs and gives odd, inconsistent,
sometime complex probability outputs. rather than the error in the current
version, or try to mimic matlab's odd output, I changed it to output a NaN for
any iscomplex() input elements. easy enough to switch back to an error if
that's preferred.

Can push to OF Statistics default if the slowdown is an okay price for the
increased accuracy, unless someone finds anything that needs to be fixed or
maybe sees a way to address the significant slowdown. 

(file #52502, file #52503)
    _______________________________________________________

Additional Item Attachment:

File name: binopdf.m                      Size:11 KB
    <https://file.savannah.gnu.org/file/binopdf.m?file_id=52502>

File name: statistics_binopdf_loader_bug34362.diff Size:13 KB
   
<https://file.savannah.gnu.org/file/statistics_binopdf_loader_bug34362.diff?file_id=52503>



    _______________________________________________________

Reply to this item at:

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

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