# [Rd] qt() returns Inf with certain negative ncp values

Martin Maechler m@ech|er @end|ng |rom @t@t@m@th@ethz@ch
Tue Jun 14 18:00:24 CEST 2022

```>>>>> GILLIBERT, Andre
>>>>>     on Tue, 14 Jun 2022 13:39:41 +0000 writes:

> Hello,
>> was suggested that they may indicate an algorithmic
>> problem with qt(), so I thought I should report them
>> here.

Which is fine.
Usually you should *CAREFULLY* read the corresponding reference
documentation before posting.

In this case, we have on R's help page {on non-central pt():}

This computes the lower tail only, so the upper tail suffers
from cancellation and a warning will be given when this is
likely to be significant.

and (in ‘Note:’)

The code for non-zero ncp is principally intended to be used
for moderate values of ncp: it will not be highly accurate,
especially in the tails, for large values.

and further also that a simple inversion is used for computing
the non-central qt().

> I explored numerical accuracy issues of pt and qt with
> non-central parameters.  There seems to be problems when
> probabilities are small (less than 10^-12 or 10^-14).

Yes, the help (above) says  "especially in the tails",
i.e., this is also well known.

> A few examples: pnorm(-30)# equal to 4.9e-198, which looks
> fine pt(-30, df=10000, ncp=0)# equal to 1e-189, which
> looks fine too pt(-30, df=10000, ncp=0.01) # equal to
> 1.044e-14, which looks bad. It should be closer to zero
> than the previous one pt(-300, df=10000, ncp=0.01) # equal
> to 1.044e-14, while it should be even closer to zero !
> pt(-3000, df=10000, ncp=0.01) # still equal to 1.044e-14,
> while it should be even closer to zero !

> qnorm(1e-13) # equal to -7.349, which looks fine qt(1e-13,
> df=10000, ncp=0) # equal to -7.359, which looks fine
> qt(1e-13, df=10000, ncp=0.01) # equal to -7.364, which
> looks fine qt(1.044e-14, df=10000, ncp=0.01) # equal to
> -8.28, which looks fine qt(1.043e-14, df=10000, ncp=0.01)
> # equal to -Inf, which is far too negative...

> The source code shows that the non-central qt() works by
> inverting the non-central pt()
> https://github.com/wch/r-source/blob/trunk/src/nmath/qnt.c

exactly; as the help page also says ..

> Consequently, both problems are related.

Indeed, and known and documented for a long time..

Still, this lack of a better algorithm had bothered me (as R
Core member) in the past quite a bit, and I had implemented other
approximations for cases where the current algorithm is
finished exploring or finding solutions in all relevant cases.

In the mean time I had created CRAN package 'DPQ' (Density,
Probability, Quantile computations) which also contains
quite a few functions related to better/alternative computations
of pt(*, ncp=*)  which I call pnt(), not the least because R's
implementation of the algorithm is in   <Rsrc>/src/nmath/pnt.c
and the C function is called pnt().

Till now, I have not found a student or a collaborator to
finally get this project further  {{hint, hint!}}.

there's a help page listing the current approaches I have

https://search.r-project.org/CRAN/refmans/DPQ/html/pnt.html
or
https://rdrr.io/cran/DPQ/man/pnt.html

implemented one, and there are even two R scripts exhibiting
bogous (and already fixed) behavior of the non-central t CDF:

https://rdrr.io/rforge/DPQ/src/tests/t-nonc-tst.R   and
https://rdrr.io/rforge/DPQ/src/tests/pnt-prec.R

Indeed, this situation *can* be improved, but it needs dedicated work
of people somewhat knowledgable in applied math etc.

Would you (readers ..) be interested in helping?

Best,
Martin

Martin Maechler
ETH Zurich  and  R Core team

PS: I'm adding code to explore this specific issue (better
inversion for those cases where pnt() is not the problem)
to my DPQ package  just these hours, notably a simple function
qtU() which only uses pt() and uniroot() to compute
(non-central) t-quantiles.

```