Getting high precision values ​​from qnorm in tail

Question

Problem

I am looking for high-precision values ​​for a normal distribution in the tail (1e-10 and 1 - 1e-10), since the R package I use sets any number that is outside this range to these values, and then calls qnormand qt.

I noticed that the implementation qnormin R is not symmetrical when looking at the tails. This is quite surprising for me, as it is well known that this distribution is symmetric, and I have seen implementations in other languages ​​that are symmetric. I checked the function qt, and is also not symmetrical in the tails.

The following are the results of the qnorm function:

x       qnorm(x)                qnorm(1-x)              qnorm(1-x) + qnorm(x)
1e-2    -2.3263478740408408     2.3263478740408408      0.0 (i.e < machine epsilon)
1e-3    -3.0902323061678132     3.0902323061678132      0.0 (i.e < machine epsilon)
1e-4    -3.71901648545568       3.7190164854557084      2.8421709430404007e-14
1e-5    -4.2648907939228256     4.2648907939238399      1.014299755297543e-12
1e-10   -6.3613409024040557     6.3613408896974208      -1.2706634855419452e-08

It is obvious that with a value xclose to 0 or 1, this function breaks. Yes, in "normal" use this is not a problem, but I look at random cases and multiplying small probabilities by very large values, in which case the error (1e-08)becomes a large value.

Note. I tried this with 1-xand by entering the actual number 0.00001and 0.99999, and the problem with accuracy still exists.

Questions

First, is this problem known with implementations qnormand qt? I could not find anything in the documentation, the algorithm should be exact 16 digits for p values ​​from 10^-314, as described in Algorithm AS 241 paper.

Quote from R doc:

Wichura, M. J. (1988) AS 241: . , 37, 477-484.

16 .

R 7- , 16 ? "", ?

R Algorithm AS 241, 16- ?

, qnorm R? , .

R

>version 
platform       x86_64-w64-mingw32          
arch           x86_64                      
os             mingw32                     
system         x86_64, mingw32             
status                                     
major          3                           
minor          3.2                         
year           2016                        
month          10                          
day            31                          
svn rev        71607                       
language       R                           
version.string R version 3.3.2 (2016-10-31)
nickname       Sincere Pumpkin Patch   

Answer 1

( R-devel-service), qnorm() , . lower.tail.

:

options(digits=22)

## For values of p in [0, 0.5], specify lower tail probabilities 
qnorm(p = 1e-10)                      ## x: P(X <= x) == 1e-10
# [1] -6.3613409024040557

## For values of p in (0.5, 1], specify upper tail probabilities
qnorm(p = 1e-10, lower.tail=FALSE)    ## x: P(X > x)  == 1e-10     (correct approach)
# [1] 6.3613409024040557
qnorm(p = 1 - 1e-10)                  ## x: P(X <= x) == 1-(1e-1)  (incorrect approach)
# [1] 6.3613408896974208

, 1-1e-10 () , 1 ( ), 1e-10 0 ( ). ( R-FAQ 7.31!) , :

1 - (1 - 1e-10) == 1e-10
## [1] FALSE

, , qnorm() (, , ) , :

qnorm(1e-314)
## [1] -37.906647423565666
qnorm(1e-314, lower.tail=FALSE)
## [1] 37.906647423565666

## With this failing in just the way (and for just the reason) you'd now expect
qnorm(1-1e-314)
# [1] Inf
1 == (1-1e-314)
# [1] TRUE