bugGNU Octave - Bugs: bug #51157, betainc.cc: Inaccurate double...

 
 

bug #51157: betainc.cc: Inaccurate double precision results

Submitter:  None
Submitted:  Thu 01 Jun 2017 12:59:31 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Schinzilord Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Tue 27 Mar 2018 12:17:34 AM UTC, comment #5: 

Fixes from GSOC2018 special function project were incorporated during OctConf 2018.  Closing report.

Rik <rik5>
Group administrator
Thu 15 Jun 2017 04:18:43 PM UTC, comment #4: 

Octave has a Google Summer of Code programmer this year who's project is improving special functions.  He will be taking a look at replacing the Fortran code with something more accurate.

Rik <rik5>
Group administrator
Sat 10 Jun 2017 01:44:37 AM UTC, comment #3: 

I checked the dispatch sequence and Octave is passing the arguments to the dbetai.f function.  This function comes from slatec (http://www.netlib.org/slatec/fnlib/).  I looked a bit more and the version of the Fortran code used by Octave is up to date.  However, the last modification was in 1992 so this is some very old code.  It appears that the Fortran code itself is the source of the innaccuracies.

Rik <rik5>
Group administrator
Sat 10 Jun 2017 01:33:22 AM UTC, comment #2: 

Just to check, I made a quick C program which uses the GNU Scientific Library to calculate betainc and it arrives at the same values as Matlab and Wolfram Alpha.  So, it would appear that the Fortran code is simply too old and it should be replaced.

Rik <rik5>
Group administrator
Thu 01 Jun 2017 04:30:45 PM UTC, comment #1: 

It's more complicated than that.  The actual code that should be called for double precision values is dbetai.f.

Perhaps there is a dispatch problem when passing execution between the Octave interpreter written in C++ (betainc.cc), to liboctave also written in C++ (lo-specfun.cc), and finally to Fortran (dbetai.f).  That would be the first thing to check.

If it is dispatching to the correct Fortran code, then maybe the slatec Fortran routines are simply too old.

Rik <rik5>
Group administrator
Thu 01 Jun 2017 12:59:31 PM UTC, original submission:  

Double precision results of function betainc seem to be incorrect.

The following example should not point out the incompatibility with Matlab, but the inaccuracy of Octave implementation.
e.g. try the following:

betainc([0.00780;0.00782;0.00784],250.005,49750.995)

Octave's results are:

ans =
   1.000000000000000
   0.999999999998649
   1.000000000000000

while e.g. Matlab or Wolfram Alpha (CDF[BetaDistribution[250.005, 49750.995], 0.00780]) yields:

ans =
   0.999999999999989
   0.999999999999992
   0.999999999999995

To my opinion Octave function betainc.cc
http://octave.org/doxygen/4.0/dd/d3c/betainc_8cc_source.html
returns double precision, while the underlying
Fortran77 code betai.f:
http://octave.org/doxygen/4.0/d1/dc2/betai_8f_source.html
returns single precision only (according to the source header).
Betainc.cc distinguishes between double and single precision.

All implemented tests performed on betainc.cc go to single precision only (sqrt(eps) in asset function), even if declared as double precision tests.

The deviations have some real life consequences for e.g. Harrell-Davis Estimator, which would return negative weights if no error handling is performed.

Can you please comment on that?

Anonymous

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by None (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 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-03-27 rik5 StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2017-06-15 rik5 Release4.2.0 dev
        Operating SystemMicrosoft Windows Any
    2017-06-10 rik5 Operating SystemAny Microsoft Windows
    2017-06-10 rik5 Operating SystemMicrosoft Windows Any
    2017-06-01 rik5 StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-cf05.
    Corresponding source code