bugGNU Octave - Bugs: bug #34362, [octave forge] (statistics)...

 
 

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

Submitter:  Rik <rik5>
Submitted:  Fri 23 Sep 2011 04:42:28 AM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  3 - Low Item Group:  Feature Request
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 27 Jul 2022 10:28:33 PM UTC, comment #12: 

Thanks for that, just pushed a fix!

Andreas Bertsatos <pr0m1th3as>
Wed 27 Jul 2022 07:11:44 PM UTC, comment #11: 

@Andreas, just noticed a tiny typo in the function. doesn't change output, but don't want the stray operator to cause confusion later:


diff --git a/inst/distributions/binopdf.m b/inst/distributions/binopdf.m
index 0024e59..80c2eff 100644
--- a/inst/distributions/binopdf.m
+++ b/inst/distributions/binopdf.m
@@ -174,7 +174,7 @@ function retval = loader_expansion (x, n, p, nx, q)
     delta_nx(!nx_precalc) = delta_fn (nx(!nx_precalc));

     ## calculate exp(log(pdf));
-    retval = exp ((delta_n - delta_x - delta_nx + -...
+    retval = exp ((delta_n - delta_x - delta_nx - ...
                        deviance (x, n .* p) - ...
                         deviance (nx, n .* q)) - ...
                          0.5 * (log(2*pi) + log (x) + log (1-x./n)));


Nicholas Jankowski <nrjank>
Group Member
Mon 25 Jul 2022 08:55:32 PM UTC, comment #10: 
Nicholas Jankowski <nrjank>
Group Member
Mon 25 Jul 2022 08:04:27 PM UTC, comment #9: 

Updated binopdf with Nicholas' new implementation.
Commit 2a6072d at https://github.com/gnu-octave/statistics/commit/59830039e91ad012e0bc17a8fe678a5272870d9c
Closed

Andreas Bertsatos <pr0m1th3as>
Tue 14 Dec 2021 05:46:58 AM UTC, 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)

Nicholas Jankowski <nrjank>
Group Member
Mon 06 Dec 2021 05:22:31 PM UTC, comment #7: 

curious - both the 2002 paper and the current R code store a number of precalculated values of formula and frequently used repeating or irrational fractions for precision and or speed. They store 28 digits:

eg:

d_n = [0.081061466795327258219670264,
  0.041340695955409294093822081, 0.0276779256849983391487892927,
  0.020790672103765093111522771, 0.0166446911898211921631948653,
  0.013876128823070747998745727, 0.0118967099458917700950557241,
  0.010411265261972096497478567, 0.0092554621827127329177286366,
  0.008330563433362871256469318, 0.0075736754879518407949720242,
  0.006942840107209529865664152, 0.0064089941880042070684396310,
  0.005951370112758847735624416, 0.0055547335519628013710386899];


it's been 20 years, so if we were to reimplement this function, is there any value to extending the length of stored numbers? future proof the code for 256bit logic :) ?

Nicholas Jankowski <nrjank>
Group Member
Fri 03 Dec 2021 08:27:15 PM UTC, comment #6: 

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

Nicholas Jankowski <nrjank>
Group Member
Fri 03 Dec 2021 06:04:22 PM UTC, comment #5: 

peeking through the Loader pdf, and updated version of which is hosted by R-project at
https://www.r-project.org/doc/reports/CLoader-dbinom-2002.pdf
as they seem to have implemented her routine ~2014 or so.

 Curiously enough, a thesis looking at the accuracy of the Loader algorithm cites the Octave 2011 mailing list for the only source of the 2000 version of her paper. Google agrees. thesis:
https://ubt.opus.hbz-nrw.de/opus45-ubtr/frontdoor/deliver/index/docId/758/file/DissertationDimitriadis.pdf which

Attaching the 2002 version here as well for archival purposes.

(file #52423)

Nicholas Jankowski <nrjank>
Group Member
Fri 07 Jun 2019 12:44:49 AM UTC, comment #4: 

The only item that I could find is the Holder algorithm
paper referred to in Rik's comment below. This is:
loader2000Fast.pdf. This looks OK but I will look a bit
more soon.

Michael Godfrey <godfrey>
Group Member
Thu 06 Jun 2019 11:58:41 PM UTC, comment #3: 

Nir: This was a mere 8 years ago! I looked for the code,
but have not found it. I will look some more...

Michael Godfrey <godfrey>
Group Member
Sat 01 Jun 2019 09:26:31 PM UTC, comment #2: 

Michael, would you mind posting your program? I would be interested in helping adapt it.

Nir Krakauer <nir_krakauer>
Fri 23 Sep 2011 06:48:28 AM UTC, comment #1: 

I just compiled a C version of the Loader algorithm
as used in R.  Simple tests get the right answers,
but it needs quite a bit of clean up.

Is it true that a C-coded version will be OK?
The code contains the following copyright:

  •  This program is free software; you can redistribute it and/or modify
  •  it under the terms of the GNU General Public License as published by
  •  the Free Software Foundation; either version 2 of the License, or
  •  (at your option) any later version.

==========================================
It is an updated version of the Catherine Loader code from
2000.

Does anyone want to take this over?

I could do a bit more clean up and then post it.

Michael

Michael Godfrey <godfrey>
Group Member
Fri 23 Sep 2011 04:42:28 AM UTC, original submission:  

The implementation of binopdf using gammaln is not particularly accurate when the number of tests, n, is large.  The Holder algorithm, see the attached paper, describes a more accurate algorithm which is already widely in use by other statistics packages such as R.

This report is being filed to make a record of the need for a better algorithm.  Anyone with Octave script file experience is welcome to try implementing the algorithm.

Rik <rik5>
Group administrator

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #52502:  binopdf.m added by nrjank (11KiB - text/plain - 1st attempt at a version of binopdf using the Loader algorithm.)
file #52503:  statistics_binopdf_loader_bug34362.diff added by nrjank (13KiB - application/octet-stream - 1st attempt at a version of binopdf using the Loader algorithm.)
file #52423:  CLoader-dbinom-2002.pdf added by nrjank (161KiB - application/pdf - updated document (2002) with the Loader algorithm)
file #24016:  loader2000Fast.pdf added by rik5 (141KiB - application/pdf)

 

Depends on the following items: None found

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -email is unavailable- added by pr0m1th3as (Posted a comment)
  • -email is unavailable- added by nrjank (Updated the item)
  • -email is unavailable- added by nir_krakauer (Posted a comment)
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by rik5 (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 14 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-07-25 nrjank StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2021-12-14 nrjank Attached File- Added binopdf.m, #52502
        Attached File- Added statistics_binopdf_loader_bug34362.diff, #52503
        StatusPostponed Patch Submitted
        Summary[octave forge] (statistics) binopdf.m: Implement high accuracy Holder algorithm for n &gt;= 10 [octave forge] (statistics) binopdf.m: Implement high accuracy Loader algorithm for n >= 10
    2021-12-03 nrjank Attached File- Added CLoader-dbinom-2002.pdf, #52423
    2019-06-04 mtmiller Summarybinopdf.m: Implement high accuracy Holder algorithm for n &gt;= 10 [octave forge] (statistics) binopdf.m: Implement high accuracy Holder algorithm for n >= 10
    2019-06-01 mtmiller Carbon-CopyRemoved 80942 -
    2019-06-01 rik5 CategoryLibraries Octave Package
    2016-11-18 mtmiller Priority5 - Normal 3 - Low
        StatusNone Postponed
    2012-03-13 rik5 Dependencies- bugs #35827 is dependent
    2011-09-23 rik5 Attached File- Added loader2000Fast.pdf, #24016

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code