bugGNU Octave - Bugs: bug #40362, chi2cdf/gamcdf/gammainc bug

 
 

bug #40362: chi2cdf/gamcdf/gammainc bug

Submitter:  None
Submitted:  Thu 24 Oct 2013 01:12:11 PM UTC
   
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Michael Leitner 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
   

Jump to the original submission

Wed 28 Mar 2018 03:11:21 AM UTC, comment #8: 

This issue looks fixed to me now with the new gammainc implementation for the upcoming Octave 4.4 release.

Mike Miller <mtmiller>
Group Member
Wed 01 Jun 2016 12:54:36 PM UTC, comment #7: 

Please, see bug # 47800. I think that the file gammainc.m there included can solve this bug.

Marco Caliari <caliari>
Group Member
Tue 10 Feb 2015 07:25:50 PM UTC, comment #6: 

This is listed as a crash but should be listed as an error.

T.J.Gaffney, If you see this message, could you show us what changes you made to d6lgit.f


Doug Stewart <dastew>
Wed 27 Nov 2013 02:40:48 PM UTC, comment #5: 

And why not steal the code from GSL? Both are GNU projects and have GPL v3 license...

José Luis García Pallero <jgpallero>
Wed 27 Nov 2013 03:45:16 AM UTC, comment #4: 

To compute the incomplete gamma function, GNU Scientific Library uses sometimes the continued fraction (which Octave uses) and sometimes a Taylor series.  (The two methods are discussed in section 5 of https://www.cs.purdue.edu/homes/wxg/selected_works/section_02/068.pdf)

I'm not able to follow how it makes the decision between the two, because it also uses a number of special cases.  However, I have checked that in these cases where Octave is failing, the series solution does converge.  By modifying d6lgit.f, I told it to try the series solution every time that the continued fraction fails.  This resolves the problem without adding any extra computation time for the cases that are working already.  Although, it wastes time in these special new cases.  So maybe a more elegant solution is needed.

Either way it comes down to:  We need to use the series computation in these cases, but Octave currently only uses the continued fraction computation.

T.J. Gaffney <gaffney2010>
Mon 25 Nov 2013 09:59:22 PM UTC, comment #3: 

In xgmainc.f, line 62, I changed ".LT." to ".LE." and this resolved the error, when the inputs are the same.

On this boundary case, the complementary incomplete Gamma function converges faster.

The problem persists when the two inputs are large and the second is a little bit larger.  In this case, neither function (incomplete Gamma or complementary incomplete Gamma) converges fast enough.

T.J. Gaffney <gaffney2010>
Fri 08 Nov 2013 12:58:37 PM UTC, comment #2: 

One option could be to use the Chi-square (and gamma related functions) from the GNU Scientific Library. I've chequed id and gets correct results

José Luis García Pallero <jgpallero>
Thu 07 Nov 2013 03:43:04 PM UTC, comment #1: 

I've also found this error, in my case using chi2inv:

>> chi2inv(0.005,1200)


 ***MESSAGE FROM ROUTINE D9LGIT IN LIBRARY SLATEC.
 ***FATAL ERROR, PROG ABORTED, TRACEBACK REQUESTED

  •  NO CONVERGENCE IN 200 TERMS OF CONTINUED FRACTION
  •  ERROR NUMBER = 3
  •  

 ***END OF MESSAGE
 
 ***JOB ABORT DUE TO FATAL ERROR.

  1.          ERROR MESSAGE SUMMARY

 LIBRARY    SUBROUTINE MESSAGE START             NERR     LEVEL     COUNT
 SLATEC     D9LGIT     NO CONVERGENCE IN 20         3         2         4

error: gammainc: exception encountered in Fortran subroutine xgammainc_
error: called from:
error:   /opt/octave/share/octave/3.6.4/m/statistics/distributions/gamcdf.m at line 58, column 12
error:   /opt/octave/share/octave/3.6.4/m/statistics/distributions/gaminv.m at line 84, column 13
error:   /opt/octave/share/octave/3.6.4/m/statistics/distributions/chi2inv.m at line 47, column 7

I've also found the value of 1200 as the starting bug point. This value is the same for my two different systems (Pentium M 32 bits and Intel Core i5).
In Matlab, chi2inv works fine regardless the value of the degrees of freedom

José Luis García Pallero <jgpallero>
Thu 24 Oct 2013 01:12:11 PM UTC, original submission:  

The octave function chi2cdf crashes if the argument is large and approximately equal to the number of degrees of freedom:


chi2cdf(1200,1200)
 ***MESSAGE FROM ROUTINE D9LGIT IN LIBRARY SLATEC.
 ***FATAL ERROR, PROG ABORTED, TRACEBACK REQUESTED
 *  NO CONVERGENCE IN 200 TERMS OF CONTINUED FRACTION
 *  ERROR NUMBER = 3
 *
 ***END OF MESSAGE

 ***JOB ABORT DUE TO FATAL ERROR.
0          ERROR MESSAGE SUMMARY
 LIBRARY    SUBROUTINE MESSAGE START             NERR     LEVEL     COUNT
 SLATEC     D9LGIT     NO CONVERGENCE IN 20         3         2        16
 SLATEC     D9LGIC     NO CONVERGENCE IN 30         1         2         8

error: gammainc: exception encountered in Fortran subroutine xgammainc_
error: called from:
error:   /usr/share/octave/3.6.2/m/statistics/distributions/gamcdf.m at line 58, column 12
error:   /usr/share/octave/3.6.2/m/statistics/distributions/chi2cdf.m at line 47, column 7


1200 is on my system the smallest crashing argument, and incidentally

chi2cdf(1200+eps(1200),1200)

works. For larger degrees of freedom, it also crashes for larger deviations from the center. Note that this is not an academic example, on the order of thousand degrees of freedom are not excessive for applications (in fact, at such high degrees of freedom asymptotic expansions should already converge rapidly enough to bypass that issue with a bit of work).

Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #29538:  chi_squared.cc added by None (772B - text/x-c++src - boost chi_squared does not seem to have the same numerical issue)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by mtmiller (Updated the item)
  • -email is unavailable- added by dastew (Posted a comment)
  • -email is unavailable- added by gaffney2010 (Posted a comment)
  • -email is unavailable- added by jgpallero (Posted a comment)
  •  

    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 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-03-28 mtmiller StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2015-02-22 mtmiller CategoryOctave Function Libraries
        Item GroupSegfault, Bus Error, etc. Incorrect Result
        StatusNone Confirmed
        Release3.6.2 dev
        Operating SystemGNU/Linux Any
    2013-11-04 None Attached File- Added chi_squared.cc, #29538

    Back to the top

    Powered by Savane 3.13-cf05.
    Corresponding source code