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

 
 

bug #40362: chi2cdf/gamcdf/gammainc bug

Submitted by:  None
Submitted on:  Thu 24 Oct 2013 01:12:11 PM UTC  
 
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Crash
Status: NoneAssigned to: None
Originator Name: Michael LeitnerOriginator Email: -unavailable-
Open/Closed: OpenRelease: 3.6.2
Operating System: GNU/Linux

Add a New Comment (Rich MarkupRich Markup):
   

You are not logged in

Please log in, so followups can be emailed to you.

 

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:

1200 is on my system the smallest crashing argument, and incidentally
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 File(s):
   
   
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
  • -unavailable- added by gaffney2010 (Posted a comment)
  • -unavailable- added by jgpallero (Posted a comment)
  •  

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follows 1 latest change.

    Date Changed By Updated Field Previous Value => Replaced By
    Mon 04 Nov 2013 09:11:00 PM UTCNoneAttached File-=>Added chi_squared.cc, #29538

    Back to the top


    Powered by Savane 3.1-cleanup