[Rd] qt() returns Inf with certain negative ncp values
GILLIBERT, Andre
Andre@G||||bert @end|ng |rom chu-rouen@|r
Tue Jun 14 15:39:41 CEST 2022
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.
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).
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
Consequently, both problems are related.
-----Message d'origine-----
De : R-devel <r-devel-bounces using r-project.org> De la part de Stephen Berman
Envoyé : mardi 14 juin 2022 11:18
À : r-devel using r-project.org
Objet : [Rd] qt() returns Inf with certain negative ncp values
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.
The Inf results below seem surprising:
> sapply(-1:-10, \(ncp) qt(1-1*(10^(-4+ncp)), 35, ncp))
[1] 3.6527153 3.0627759 2.4158355 1.7380812 1.0506904 0.3700821
[7] Inf -0.9279783 -1.5341759 -2.1085213
> sapply(seq(-6.9, -7.9, -0.1), \(ncp) qt(1-1*(10^(-4+ncp)), 35, ncp))
[1] -0.2268386 Inf Inf Inf -0.4857400 -0.5497784
[7] -0.6135402 -0.6770143 -0.7401974 -0.8030853 -0.8656810
These inputs also yield many repetitions of the following warning
message:
In qt(1 - 1 * (10^(-4 + ncp)), 35, ncp) :
full precision may not have been achieved in 'pnt{final}'
In particular, in the range -1:-10 I don't get this warning with ncp =
-1 through -4, but I do get it once with each of -5 and -8 through -10,
32 times with -6, and 50 times with -7. In the range -6.9:-7.9 I get the warning twice with each of -6.9 and -7.3 through -7.7, once with
-7.8 and -7.9, and 50 times with each of -7.0 through -7.2.
In case it matters:
> sessionInfo()
R Under development (unstable) (2022-06-05 r82452)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Linux From Scratch r11.0-165
Matrix products: default
BLAS: /usr/lib/R/lib/libRblas.so
LAPACK: /usr/lib/R/lib/libRlapack.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] compiler_4.3.0 tools_4.3.0
Thanks.
Steve Berman
______________________________________________
R-devel using r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
More information about the R-devel
mailing list