help-octave
[Top][All Lists]

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

 From: John W. Eaton Subject: potential bug in betai; any hints at a fix or workaround? Date: Thu, 29 Jan 1998 17:27:41 -0600

```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:
|
| ## Created: 2 August 1994
|
| 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