问题
我正在寻找尾部(1e-10 and 1 - 1e-10)
中正态分布的高精度值,因为我正在使用的R包将超出此范围的任何数字设置为这些值,然后调用qnorm
和qt
函数。
我注意到的是,在看尾部时,R中的qnorm
实现不是对称的。这对我来说是非常令人惊讶的,因为众所周知这种分布是对称的,而且我已经看到其他对称语言的实现。我已经检查了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
实现的已知问题吗?我在文档中找不到任何内容,对于Algorithm AS 241文件中所述的10^-314
,p值的算法应该是准确的16位数字。来自R doc的引用:
如果R代码实现7位数字版本,为什么要声明16位数字呢?还是“准确”但原始算法不对称且错误?
如果R同时实现了Algorithm AS 241的两个版本,我可以打开16位数字的版本吗?
或者,R中是否有更准确的
qnorm
版本?或者,解决我的问题的另一种解决方案,我需要在尾部使用高精度以实现分位数功能。
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
最佳答案
事实证明(正如Spencer Graves在his response中针对R-devel列表服务器上的相同问题所指出的那样)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
关于r - 从尾部的qnorm获取高精度值,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/43362644/