From bug-octave-request at bevo dot che dot wisc dot edu  Thu Jan 29 17:28:00 1998
Subject: potential bug in betai; any hints at a fix or workaround?
From: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
To: king at cogsci dot ucsd dot edu (Jonathan King)
Cc: help-octave at bevo dot che dot wisc dot edu, bug-octave@bevo.che.wisc.edu
Date: Thu, 29 Jan 1998 17:27:41 -0600

On 29-Jan-1998, Jonathan King <king at cogsci dot ucsd dot edu> 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 <Kurt dot Hornik at ci dot tuwien dot ac dot at>
| ## 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