从尾部的qnorm获得高精度值

Question

问题

我正在寻找尾部正态分布的高精度值(1e-10 and 1 - 1e-10),因为我正在使用的R包将任何超出此范围的数字设置为这些值,然后调用qnormand qt函数.

我注意到的是,qnorm在查看尾部时,R 中的实现是不对称的.这对我来说是非常令人惊讶的,因为众所周知这个分布是对称的,我已经看到其他语言中的实现是对称的.我已检查过该qt功能,它在尾部也不对称.

以下是qnorm函数的结果:

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

很明显,在x接近0或1 的值时,此功能会崩溃.是的,在"正常"使用中这不是问题,但我正在研究边缘情况并将小概率乘以非常大的值,在这种情况下,误差(1e-08)变为大值.

注意:我已尝试使用1-x并输入实际数字0.00001,0.99999并且准确性问题仍然存在.

问题

首先,这是一个已知的问题qnorm和qt实现？我在文档中找不到任何内容,算法应该是p值的准确16位数,10^-314如算法AS 241论文中所述.

引自R doc:

Wichura,MJ(1988)Algorithm AS 241:正态分布的百分点.Applied Statistics,37,477-484.

它提供高达约16位的精确结果.

如果R代码实现7位数版本,为什么它要求16位数字？或者它是"准确的",但原始算法不对称和错误？

如果R确实实现了两个版本的算法AS 241,我可以打开16位数版本吗？

或者,qnormR中是否有更准确的版本？或者,我的问题的另一个解决方案,我需要在分位数函数的尾部高精度.

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

事实证明（正如斯宾塞·格雷夫斯（Spencer Graves）在R-devel list-serve 上对同一问题的回答中所指出的那样）实际上确实如广告中所宣传的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