[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
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 29Jan1998, Jonathan King <address@hidden> wrote:
 I just convinced somebody around here to start using Octave (2.0.9)
 on i586pclinuxgnu. 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 = 1e20;
 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^.5eps; 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
 handwritten 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
 ## AdaptedBy: jwe

 with a 1995, 1996 copyright by the author. Is there a fix for this?
 An obvious workaround?
For Octave 2.0.10, the Mfile 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.6902e13 1.9771e23
Is that about right?
Thanks,
jwe