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

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

## Post a Comment

Add a New Comment

Comment Type & Canned Response:

## Discussion

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.

Nicholas Jankowski <nrjank>
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>
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>
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.

Nicholas Jankowski <nrjank>
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>
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>
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>
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>

## Attached Files

(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)

## Dependencies

Depends on the following items: None found

Items that depend on this one

Digest:
bug dependencies.

## Mail Notification Carbon-Copy List

Carbon-Copy List
• -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)
•

## Votes

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 project members can vote.

## History

Follow 12 latest changes.

Date Changed by Updated Field Previous Value => Replaced by
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 >= 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 >= 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 Forge 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

Powered by Savane 3.9