problem with module inc_beta for incomplete beta function

[salbruce]

I have been using Alan Miller's inc_beta module (for calculating the incomplete beta function) for a while with no problems, but recently encountered a situation where it seems to fail. The source code is online at http://jblevins.org/mirror/amiller/bratio.f90

I'm hoping that someone more familiar with this routine than I am can suggest a fix.

The application is calculating a p-value for a large F-ratio (e.g. F=1363), and the sample size is also large (10,000). There is no question in this case that p will be very close to zero, but if nothing else, it seems like p should return as epsilon or zero, rather than choke. What does happen is that the STOP in the code is encountered:

131 CALL bgrat (b0, a0, y0, x0, w1, eps, ierr1)
IF (ierr1 > 0) STOP "Error in BGRAT"

I've made that exit more gracious, but would like a better fix.

The input values to subr. bgrat are:
0.500000000000000 = b
4999.00000000000 = a
0.255622364748191 = y
0.744377635251809 = x
0.000000000000000E+000 = w

Then the following condition is encountered, the routine bails out.

IF (u == zero) GO TO 100   

0.000000000000000E+000 = u

Details:
Intel Visual Fortran Compiler Integration for Microsoft Visual Studio* 2008, 11.1.3468.2008
Windows Vista Home SP 2
Intel Core 2 Quad CPU, 32 bit
Problem happens with or without compiler optimizations
I am compiling a dll.
Typical compiler settings:
/nologo /I"Release/" /DNICHERR=1 /arch:IA32 /warn:unused /warn:interfaces /fpconstant /iface:cvf /module:"Release/" /object:"Release/" /check:uninit /libs:static /threads /winapp /c

Thanks in advance for any suggestions.

-Bruce

[TimP]

It certainly appears that a zero result may be produced in extreme cases such as yours. Perhaps you might try changing the code so as to accept the zero result. /fp:source would seem a likely addition to your options, as you apparently don't care much about performance.
With due respect to Alan Miller, the gratuitious mixture of f66 and f90 styles makes the code unnecessarily difficult to peruse.

[mecej4]

The comments at the beginning of subroutine BGRAT state:

ASYMPTOTIC EXPANSION WHEN A IS > B. IT IS ASSUMED THAT A <= 15 AND B <= 1

Masking the error return for arguments outside the range of validity of the algorithm is probably a bad idea, as it introduces a bug that will be hidden to posterity.

See if the TOMS 179 code suits your needs.

[salbruce]

Thanks for the suggestions. I had seen the comment in BGRAT, but had not had problems with A > 15 before. For example, the following does not give an error in BGRAT:

0.500000000000000 = b7

1999.00000000000 = a

0.254162371158600 = y

0.745837628841400 = x

I will check out the TOMS 179 routine.

Thanks,
Bruce

[mecej4]

> I had seen the comment in BGRAT, but had not had problems with A > 15 before.

One should not assume that receiving no error messages is evidence of having no problems. Verification is needed to check that special function values returned by an approximation have enough precision.

Even with more benign values, a=5, b=0.5, x=0.2, y=1-x, in your notation, the incomplete beta function values given by Miller's routine and the one from Numerical Recipes agree to only eight significant digits when double precision calculations are employed. The NR result is less accurate.

BRATIO.f90: 0.85507239 4595920
NR : 0.85507239 9728096

[salbruce]

Agreed. Thanks for providing this comparison. I have checked results for various input values against regression results in SPSS and have found good agreement, at least to the precision you indicated. I am fine with that level of precision for my applications.

Do you use the NR routine a lot and have you been happy with it? Does it handle the more extreme values ok?

Thanks again,
Bruce

[mecej4]

> Do you use the NR routine a lot and have you been happy with it? Does it >handle the more extreme values ok?

I used NR in the mid-90s; like fast food, it was cheap, had variety and was always available. Since then, I gave up on NR and fast food for similar reasons.

I do not use any software, not just NR, without running some sanity checks. For special functions the classic Jahnke-Emde and Abramowicz-Stegun books are very useful as sources of reference values against which one can check software. These days, you may also use Maple, Matlab, etc. for such purposes.