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

[Martin Maechler]

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

    > Hello,
    >> I asked about the following observations on r-help and it
    >> 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
deficient... but I had not been entirely satisfied, nor had I
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!}}.

In DPQ, (download the *source* package if you are interested),
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

Additionally, in the source (man/pnt.Rd) there are comments about a not yet
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.