> GILLIBERT, Andre
> on Tue, 14 Jun 2022 13:39:41 + 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=1, ncp=0)# equal to 1e-189, which
> looks fine too pt(-30, df=1, ncp=0.01) # equal to
> 1.044e-14, which looks bad. It should be closer to zero
> than the previous one pt(-300, df=1, ncp=0.01) # equal
> to 1.044e-14, while it should be even closer to zero !
> pt(-3000, df=1, 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=1, ncp=0) # equal to -7.359, which looks fine
> qt(1e-13, df=1, ncp=0.01) # equal to -7.364, which
> looks fine qt(1.044e-14, df=1, ncp=0.01) # equal to
> -8.28, which looks fine qt(1.043e-14, df=1, 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 /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.
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
