R devel - Dec 2008 - function qt can fails if ndf < 1 (PR#13364)

Full_Name: Gerard Torrent
Version: R version 2.8.0 (2008-10-20)
OS: Linux 2.6.27.5-41.fc9.x86_64 #1 SMP
Submission from: (NULL) (85.52.227.233)


In some cases qt complains about NaNs and don't gives the correct result:
> qt(0.1,, 0.1)
[1] NaN
Warning message:
In qt(p, df, lower.tail, log.p) : NaNs produced

But the result can be found:
> pt(-1.60443e+06, 0.1)
[1] 0.09999997

If I replace the current method by bisection method in file file src/nmath/qt.c
the result is computed fine.

int cont=0;
do {
	cont++;
	nx = lx + (ux - lx)/2.0;
	if ((pp = pt(nx, ndf, TRUE, FALSE)) > p) {
		ux = nx; pu = pp;
	} else {
		lx = nx; pl = pp;
	}
}
while ((ux - lx) / fabs(nx) > accu && cont < 100);

You have missed the following item in NEWS for R-patched.

NEW FEATURES

     o   qt() now works for 0 < df < 1.


Please review what the FAQ has to say about checkiong your facts and not 
reporting things which are already fixed.

Fortunately this example had some use: it showed someone had replaced a 
working method in R-devel by one that infinite-looped, and I've reverted 
it.

On Wed, 3 Dec 2008, gerard at mail.generacio.com wrote:
> Full_Name: Gerard Torrent
> Version: R version 2.8.0 (2008-10-20)
> OS: Linux 2.6.27.5-41.fc9.x86_64 #1 SMP
> Submission from: (NULL) (85.52.227.233)
>
>
> In some cases qt complains about NaNs and don't gives the correct
result:
>
>> qt(0.1,, 0.1)
That gives a syntax error ....
> [1] NaN
> Warning message:
> In qt(p, df, lower.tail, log.p) : NaNs produced
>
> But the result can be found:
>> pt(-1.60443e+06, 0.1)
> [1] 0.09999997
It actually was documented so in 2.8.0

       df: degrees of freedom (> 0, maybe non-integer).  'df = Inf'
is
           allowed.  For 'qt' only values of at least one are currently
           supported unless 'ncp' is supplied.

!
> If I replace the current method by bisection method in file file
src/nmath/qt.c
> the result is computed fine.
>
> int cont=0;
> do {
> 	cont++;
> 	nx = lx + (ux - lx)/2.0;
> 	if ((pp = pt(nx, ndf, TRUE, FALSE)) > p) {
> 		ux = nx; pu = pp;
> 	} else {
> 		lx = nx; pl = pp;
> 	}
> }
> while ((ux - lx) / fabs(nx) > accu && cont < 100);
>
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595