我经常遇到必须检查所获得的差异是否超出机器精度的情况。为此,似乎R具有一个方便的变量:.Machine$double.eps
。但是,当我转向R源代码获取有关使用此值的准则时,会看到多种不同的模式。
例子
以下是stats
库中的一些示例:
检验
if(stderr < 10 *.Machine$double.eps * abs(mx))
chisq.test.R
if(abs(sum(p)-1) > sqrt(.Machine$double.eps))
整合
rel.tol < max(50*.Machine$double.eps, 0.5e-28)
影响力
e[abs(e) < 100 * .Machine$double.eps * median(abs(e))] <- 0
PRINCOMR
if (any(ev[neg] < - 9 * .Machine$double.eps * ev[1L]))
等等
问题
10 *
,100 *
,50 *
和sqrt()
修饰符背后的原因? .Machine$double.eps
调整由于精度问题引起的差异的准则? 最佳答案
double
的机器精度取决于其当前值。 .Machine$double.eps
在值为1时给出精度。您可以使用C函数nextAfter
获得其他值的机器精度。
library(Rcpp)
cppFunction("double getPrec(double x) {
return nextafter(x, std::numeric_limits<double>::infinity()) - x;}")
(pr <- getPrec(1))
#[1] 2.220446e-16
1 + pr == 1
#[1] FALSE
1 + pr/2 == 1
#[1] TRUE
1 + (pr/2 + getPrec(pr/2)) == 1
#[1] FALSE
1 + pr/2 + pr/2 == 1
#[1] TRUE
pr/2 + pr/2 + 1 == 1
#[1] FALSE
当
a
是b
的一半机器精度时,将b
值添加到a
值不会更改<=
。使用<
检查差异是否比机器精度小。修饰符可能会考虑典型情况,即添加项未显示更改的频率。在R中,可以通过以下方式估算机器精度:
getPrecR <- function(x) {
y <- log2(pmax(.Machine$double.xmin, abs(x)))
ifelse(x < 0 & floor(y) == y, 2^(y-1), 2^floor(y)) * .Machine$double.eps
}
getPrecR(1)
#[1] 2.220446e-16
每个
double
值代表一个范围。对于一个简单的加法,结果的范围取决于每个被加数的范围以及它们的和的范围。library(Rcpp)
cppFunction("std::vector<double> getRange(double x) {return std::vector<double>{
(nextafter(x, -std::numeric_limits<double>::infinity()) - x)/2.
, (nextafter(x, std::numeric_limits<double>::infinity()) - x)/2.};}")
x <- 2^54 - 2
getRange(x)
#[1] -1 1
y <- 4.1
getRange(y)
#[1] -4.440892e-16 4.440892e-16
z <- x + y
getRange(z)
#[1] -2 2
z - x - y #Should be 0
#[1] 1.9
2^54 - 2.9 + 4.1 - (2^54 + 5.9) #Should be -4.7
#[1] 0
2^54 - 2.9 == 2^54 - 2 #Gain 0.9
2^54 - 2 + 4.1 == 2^54 + 4 #Gain 1.9
2^54 + 5.9 == 2^54 + 4 #Gain 1.9
对于更高的精度,可以使用
Rmpfr
。library(Rmpfr)
mpfr("2", 1024L)^54 - 2.9 + 4.1 - (mpfr("2", 1024L)^54 + 5.9)
#[1] -4.700000000000000621724893790087662637233734130859375
如果可以将其转换为整数,则可以使用
gmp
(Rmpfr中的内容)。library(gmp)
as.bigz("2")^54 * 10 - 29 + 41 - (as.bigz("2")^54 * 10 + 59)
#[1] -47
关于r - 检查差异是否小于机器精度的正确/标准方法是什么?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/59229545/