[R] Why does qt() return Inf with certain negative ncp values?
Martin Maechler
m@ech|er @end|ng |rom @t@t@m@th@ethz@ch
Tue Jun 14 17:49:47 CEST 2022
>>>>> on Mon, 13 Jun 2022 19:28:43 -0700 writes:
> Well, this will likely require close numerical analysis of
> the algorithm used for the calculation, which I can't help
> with, but I do note for the first case (I didn't bother
> checking the second) it may be useful to note that besides
> returning Inf, I also got with -7 as ncp:
> "There were 50 or more warnings" all of which were: "In
> qt(1 - 1 * (10^(-11)), 35, -7) : full precision may not
> have been achieved in 'pnt{final}' "
> I got this warning for all the other values of ncp also.
> So it looks like: Congratulations, you seem to have
> defeated the algorithm.
> Bert Gunter
Thank you, Bert, for the extra information and Stephen for
reporting.
[@Stephen: no reason to use sapply(); qt() etc are all vectorized!]
Part of this is "known in principle"
as the help page mentions in several places that the non-central t
algorithm for pt() has deficiencies and
that in the non-central case, qt() is defined by inversion,
i.e., a (primitive !) version of uniroot(), using pt(*, ncp=.) :
> 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 indeed, you use values very much "in the tails".
and that *does* give a warning, as Bert mentioned.
[ You really should heed such warnings.
{{careless people using suppressWarnings() all the time is
really very irresponsible programming/scientific practice!
}}]
>>>>> on Mon, 13 Jun 2022 19:28:43 -0700 writes:
> Well, this will likely require close numerical analysis of
> the algorithm used for the calculation, which I can't help
> with, but I do note for the first case (I didn't bother
> checking the second) it may be useful to note that besides
> returning Inf, I also got with -7 as ncp:
> "There were 50 or more warnings" all of which were: "In
> qt(1 - 1 * (10^(-11)), 35, -7) : full precision may not
> have been achieved in 'pnt{final}' "
> I got this warning for all the other values of ncp also.
> So it looks like: Congratulations, you seem to have
> defeated the algorithm.
> Bert Gunter
Thank you, Bert, for the extra information and Stephen for
reporting.
[@Stephen: no reason to use sapply(); qt() etc are all vectorized!]
Part of this is "known in principle"
as the help page mentions in several places that the non-central t
algorithm for pt() has deficiencies and
that in the non-central case, qt() is defined by inversion,
i.e., a (primitive !) version of uniroot(), using pt(*, ncp=.) :
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 indeed, you use values very much "in the tails".
and the use *does* give warning, as Bert mentioned.
You really should heed such warnings.
{{careless people using suppressWarnings() all the time is
really very irresponsible programming/scientific practice! }}
There's much more to be said here, but as I see Stephen has in
the mean time diverted the topic to the R-devel mailing list --
quite appropriately --
--> https://stat.ethz.ch/pipermail/r-devel/2022-June/081785.html
I will answer there.
Best,
Martin
--
Martin Maechler
ETH Zurich and R Core team
More information about the R-help
mailing list