检查差异是否小于机器精度的 correct/standard 方法是什么？
What is the correct/standard way to check if difference is smaller than machine precision?
我经常遇到需要检查获得的差异是否高于机器精度的情况。似乎为此目的，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))
integrate.R
rel.tol < max(50*.Machine$double.eps, 0.5e-28)
lm.influence.R
e[abs(e) < 100 * .Machine$double.eps * median(abs(e))] <- 0
princomp.R
if (any(ev[neg] < - 9 * .Machine$double.eps * ev[1L]))
等等
问题
- 如何理解所有这些不同的
10 *、100 *、50 * 和 sqrt() 修饰符背后的原因？
- 是否有关于使用
.Machine$double.eps 调整因精度问题引起的差异的指南？
machine.eps 的定义：它是 eps 的最低值 1+eps 不是 1
根据经验（假设以 2 为底的浮点表示）：
eps 使范围 1 .. 2,
有所不同
对于范围 2 .. 4，精度为 2*eps
等等。
不幸的是，这里没有好的经验法则。这完全取决于你的程序的需要。
在 R 中，我们有 all.equal 作为测试近似相等性的内置方法。所以你可以使用类似 (x<y) | all.equal(x,y)
i <- 0.1
i <- i + 0.05
i
if(isTRUE(all.equal(i, .15))) { #code was getting sloppy &went to multiple lines
cat("i equals 0.15\n")
} else {
cat("i does not equal 0.15\n")
}
#i equals 0.15
Google mock有多个floating point matchers用于双精度比较，包括DoubleEq和DoubleNear。您可以像这样在数组匹配器中使用它们：
ASSERT_THAT(vec, ElementsAre(DoubleEq(0.1), DoubleEq(0.2)));
更新：
Numerical Recipes 提供了一个推导来证明使用单边差分商，sqrt 是导数的有限差分近似步长的一个很好的选择。
维基百科文章站点 Numerical Recipes，第 3 版，第 5.7 节，第 229-230 页（有限的页面浏览量可在 http://www.nrbook.com/empanel/ 获得）。
all.equal(target, current,
tolerance = .Machine$double.eps ^ 0.5, scale = NULL,
..., check.attributes = TRUE)
这些 IEEE floating point arithmetic 是众所周知的计算机算法限制，在几个地方进行了讨论：
- R 中的常见问题解答有一个专门针对它的问题：R FAQ 7.31
- The R Inferno by Patrick Burns 将第一个 "Circle" 用于解决此问题（第 9 页以后）
- Arithmetic Sum Proof Problem 在 Math Meta
中提问
- David Goldberg，"What Every Computer Scientist Should Know About Floating-point Arithmetic," ACM 计算调查 23, 1 (1991-03), 5-48 doi>10.1145/103162.103163 (revision also available)
- The Floating-Point Guide - What Every Programmer Should Know About Floating-Point Arithmetic
- 0.30000000000000004.com 比较各种编程语言的浮点运算
- Canonical duplicate for "floating point is inaccurate"（关于此问题的规范答案的元讨论）
- 几个 Stack Overflow 问题，包括
- Why can't decimal numbers be represented exactly in binary?
- Why are floating point numbers inaccurate?
- Is floating point math broken?
- Arthur T. Benjamin
解释的数学技巧
。
dplyr::near() 是测试两个浮点数向量是否相等的另一种选择。
该函数有一个内置的公差参数：tol = .Machine$double.eps^0.5，可以调整。默认参数与 all.equal().
的默认值相同
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 是 <= 机器精度的一半时，将值 a 添加到值 b 不会更改 b。 使用 < 检查差异是否小于机器精度。修饰符可能会考虑典型案例，添加项未显示更改的频率。
在 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 有一个方便的变量：.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))
integrate.R
rel.tol < max(50*.Machine$double.eps, 0.5e-28)
lm.influence.R
e[abs(e) < 100 * .Machine$double.eps * median(abs(e))] <- 0
princomp.R
if (any(ev[neg] < - 9 * .Machine$double.eps * ev[1L]))
等等
问题
- 如何理解所有这些不同的
10 *、100 *、50 *和sqrt()修饰符背后的原因？ - 是否有关于使用
.Machine$double.eps调整因精度问题引起的差异的指南？
machine.eps 的定义：它是 eps 的最低值 1+eps 不是 1
根据经验（假设以 2 为底的浮点表示）：
eps 使范围 1 .. 2,
有所不同
对于范围 2 .. 4，精度为 2*eps
等等。
不幸的是，这里没有好的经验法则。这完全取决于你的程序的需要。
在 R 中，我们有 all.equal 作为测试近似相等性的内置方法。所以你可以使用类似 (x<y) | all.equal(x,y)
i <- 0.1
i <- i + 0.05
i
if(isTRUE(all.equal(i, .15))) { #code was getting sloppy &went to multiple lines
cat("i equals 0.15\n")
} else {
cat("i does not equal 0.15\n")
}
#i equals 0.15
Google mock有多个floating point matchers用于双精度比较，包括DoubleEq和DoubleNear。您可以像这样在数组匹配器中使用它们：
ASSERT_THAT(vec, ElementsAre(DoubleEq(0.1), DoubleEq(0.2)));
更新：
Numerical Recipes 提供了一个推导来证明使用单边差分商，sqrt 是导数的有限差分近似步长的一个很好的选择。
维基百科文章站点 Numerical Recipes，第 3 版，第 5.7 节，第 229-230 页（有限的页面浏览量可在 http://www.nrbook.com/empanel/ 获得）。
all.equal(target, current,
tolerance = .Machine$double.eps ^ 0.5, scale = NULL,
..., check.attributes = TRUE)
这些 IEEE floating point arithmetic 是众所周知的计算机算法限制，在几个地方进行了讨论：
- R 中的常见问题解答有一个专门针对它的问题：R FAQ 7.31
- The R Inferno by Patrick Burns 将第一个 "Circle" 用于解决此问题（第 9 页以后）
- Arithmetic Sum Proof Problem 在 Math Meta 中提问
- David Goldberg，"What Every Computer Scientist Should Know About Floating-point Arithmetic," ACM 计算调查 23, 1 (1991-03), 5-48 doi>10.1145/103162.103163 (revision also available)
- The Floating-Point Guide - What Every Programmer Should Know About Floating-Point Arithmetic
- 0.30000000000000004.com 比较各种编程语言的浮点运算
- Canonical duplicate for "floating point is inaccurate"（关于此问题的规范答案的元讨论）
- 几个 Stack Overflow 问题，包括
- Why can't decimal numbers be represented exactly in binary?
- Why are floating point numbers inaccurate?
- Is floating point math broken?
- Arthur T. Benjamin 解释的数学技巧
。
dplyr::near() 是测试两个浮点数向量是否相等的另一种选择。
该函数有一个内置的公差参数：tol = .Machine$double.eps^0.5，可以调整。默认参数与 all.equal().
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 是 <= 机器精度的一半时，将值 a 添加到值 b 不会更改 b。 使用 < 检查差异是否小于机器精度。修饰符可能会考虑典型案例，添加项未显示更改的频率。
在 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