检查差异是否小于机器精度的正确/标准方法是什么？

Question

我经常遇到需要检查获得的差异是否高于机器精度的情况。似乎为此目的，R 有一个方便的变量：.Machine$double.eps. 但是，当我转向 R 源代码以获取有关使用此值的指南时，我看到了多种不同的模式。

例子

以下是stats库中的一些示例：

t.test.R

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

原理图

if (any(ev[neg] < - 9 * .Machine$double.eps * ev[1L]))

等等。

问题

1. 一个如何理解这些不同背后的理由10 *，100 *，50 *和sqrt()修饰？
2. 是否有关于.Machine$double.eps用于调整因精度问题而导致的差异的指南？

Answer 1

machine.eps的定义：它是最低值\xc2\xa0 eps\xc2\xa0，其中\xc2\xa0 1+eps\xc2\xa0不是\xc2\xa01

\n\n

根据经验（假设以 2 为基数的浮点表示）：
\n这eps会导致范围 1 .. 2 的差异，
\n对于范围 2 .. 4 的精度为2*eps
\n，依此类推。

\n\n

不幸的是，这里没有好的经验法则。这完全取决于您的程序的需求。

\n\n

在 R 中，我们有 all.equal 作为测试近似相等的内置方法。所以你可以使用类似的东西(x<y) | all.equal(x,y）

\n\n

i <- 0.1\n i <- i + 0.05\n i\nif(isTRUE(all.equal(i, .15))) { #code was getting sloppy &went to multiple lines\n    cat("i equals 0.15\\n") \n} else {\n    cat("i does not equal 0.15\\n")\n}\n#i equals 0.15\n

\n\n

Google mock 有许多用于双精度比较的浮点匹配器，包括DoubleEq和DoubleNear。您可以在数组匹配器中使用它们，如下所示：

\n\n

ASSERT_THAT(vec, ElementsAre(DoubleEq(0.1), DoubleEq(0.2)));\n

\n\n

更新：

\n\n

数值配方提供了一种推导来证明使用单边差商，sqrt是导数有限差分近似的步长的良好选择。

\n\n

维基百科文章网站 Numerical Recipes，第 3 版，第 5.7 节，第 229-230 页（有限数量的页面浏览量可在http://www.nrbook.com/empanel/上获得）。

\n\n

all.equal(target, current,\n           tolerance = .Machine$double.eps ^ 0.5, scale = NULL,\n           ..., check.attributes = TRUE)\n

\n\n

这些IEEE 浮点运算是计算机运算的众所周知的限制，并在多个地方进行了讨论：

\n\n

\n
• R 中的常见问题解答有一个专门针对它的完整问题：R FAQ 7.31 \n\n
  \n
  • Patrick Burns 的 R Inferno将第一个“Circle”专门用于解决这个问题（第 9 页起）
  \n
  • Math Meta中提出的算术和证明问题
  \n
  • David Goldberg，“每个计算机科学家应该了解浮点算术知识”，ACM 计算调查 23 , 1 (1991-03), 5-48 doi>10.1145/103162.103163（也可进行修订）
  \n
  • 浮点指南 - 每个程序员都应该了解的浮点运算知识
  \n
  • 0.30000000000000004.com比较不同编程语言的浮点运算
  \n
\n
• “浮点不准确”的规范重复（关于此问题的规范答案的元讨论）
\n
• 几个 Stack Overflow 问题，包括\n\n
  \n
  • 为什么十进制数不能用二进制精确表示？
  \n
  • 为什么浮点数不准确？
  \n
  • 浮点数学有问题吗？
  \n
\n
• 亚瑟·T·本杰明解释的数学技巧
\n

\n\n

.\n dplyr::near()是测试两个浮点数向量是否相等的另一个选项。

\n\n

该函数有一个内置的容差参数：tol = .Machine$double.eps^0.5可以调整。默认参数与 的默认参数相同all.equal()。

\n

Answer 2

机器精度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