[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
## potential bug in betai; any hints at a fix or workaround?

**From**: |
Jonathan King |

**Subject**: |
potential bug in betai; any hints at a fix or workaround? |

**Date**: |
Thu, 29 Jan 1998 15:10:48 -0800 |

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?
Thanks.
jking

**potential bug in betai; any hints at a fix or workaround?**,
*Jonathan King* **<=**