Bug and Fix for Inverse Chi-Squared Method

[Janes, Robert (Columbus)]

Hello,

I ran into a pretty serious issue with the cdf_chisq_Pinv method.  Thank 
goodness for open-source because I was able to come up with a simple fix.

In case Outlook considers a makefile sort of an executable I've copied it here 
rather than attach it.  You may have to add tabs back into it.
sample: testing
./testing

testing: obj/testing.o obj/imsl_c.o
gfortran -o testing obj/imsl_c.o obj/testing.o \
  -L/home/rjanes/gsl/gsl-2.6/.libs -l:libgsl.a -lgslcblas

obj/%.o: %.f
gfortran -c -g -o obj/$*.o $*.f

obj/%.o: %.c
gcc -c -g -o obj/$*.o $*.c

The problem is that for large degrees of freedom the method will abort with an 
error, "inverse failed to converge".  In looking at the code and running test 
cases I found that the limit of 32 iterations is too small.  I found the method 
failed many times once degrees of freedom got above 1000 or so.  I increased 
the iteration limit to 50 and did not find a single failure for a huge range of 
degrees of freedom.  I found successful iterations in this range were taking 
28, 29 and 31 iterations.  So hitting 32 maybe isn't much of a surprise.  But 
it suggests that maybe this isn't too efficient an iteration, at least for 
large degrees of freedom.  But I'll leave that to you experts.

Sorry about the Fortran code.  I'm porting a large application from Intel 
Fortran / IMSL to gfortran / gsl.  Rather than change the base application 
code, I just implemented the IMSL calls with gsl.

Of the attached files, testing.f is the test Fortran program.  imsl.c is the 
implementation of the IMSL routine with gsl.  And gammainv.c is the updated gsl 
code.  The only change was 32 -> 50.

Let me know of you need anything else from me.  Thank you.


Rob Janes

SGS Transportation

Senior Software Engineer

SGS – CMX

2860 N. National Road Suite A

US – 47201 – Columbus, Indiana

Phone: +1 - 812 - 378 - 7966

Fax: +1 812 - 378 - 3393

Email: address@hidden



Information in this email and any attachments is confidential and intended 
solely for the use of the individual(s) to whom it is addressed or otherwise 
directed. Please note that any views or opinions presented in this email are 
solely those of the author and do not necessarily represent those of the 
Company. Finally, the recipient should check this email and any attachments for 
the presence of viruses. The Company accepts no liability for any damage caused 
by any virus transmitted by this email. All SGS services are rendered in 
accordance with the applicable SGS conditions of service available on request 
and accessible at https://www.sgs.com/en/terms-and-conditions

gammainv.c
Description: gammainv.c

imsl_c.c
Description: imsl_c.c

testing.f
Description: testing.f