potential bug in betai; any hints at a fix or workaround?

[John W. Eaton]

On 29-Jan-1998, Jonathan King <address@hidden> wrote:

| I just convinced somebody around here to start using Octave (2.0.9)
| on i586-pc-linux-gnu.  He was very happy with it, until he found
| what appears to be a bug while writing his very first octave
| function, based on some previous (working) C code:
| 
| **start code here**
| 
| function p = rtop2 (r,df)
| 
| % RTOP(r, df) computes p values from a vector of correlations coefficients
| %
| % Usage: p = rtop2(r, df)
| %
| % where r is a vector of correlation coefficients, and df is a scalar
| % that gives the required degrees of freedom for each r.
| 
| EPSILON = 1e-20;
| one = 1 + EPSILON;
| 
| tsq = r .* r .* (df ./ ((one + r) .* (one -r)));
| p = betai(0.5*df, 0.5, (df ./ (df + tsq)));
| 
| p;
| 
| **end code here
| 
| Everything looked fine until he called it with a vector including
| values of in the range of (+/-)[.56141,.70710] with a df of 144.
| (The upper limit of the range is actually 1/2^.5-eps; the lower
| limit appears to be just about e^(-1/3^.5).)  In this range, he
| started to get answers less than zero, which is, um, wrong, since
| these are probabilities, albeit tiny ones.
| 
| Betai seems to be the problem here; both Matlab 5.1 and a
| hand-written C version essentially yanked from Numerical Recipes
| provided what seem to be the correct answers.
| 
| Looking at the code for betai.m, it's not immediately clear to me
| where the real problem is; the version we have is stated to be:
| 
| ## Author: KH <address@hidden>
| ## Created: 2 August 1994
| ## Adapted-By: jwe
| 
| with a 1995, 1996 copyright by the author.  Is there a fix for this?
| An obvious work-around?

For Octave 2.0.10, the M-file implemntation for betai is replaced by a
function from the Slatec library.  Here's what I get with my current
sources:

  octave:4> rtop2 ([.56141,.70710], 144)
  ans =

     1.6902e-13   1.9771e-23

Is that about right?

Thanks,

jwe