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

问题描述：

我经常遇到必须检查所获得的差异是否超出机器精度的情况.为此，似乎R具有一个方便的变量:.Machine$double.eps.但是，当我转向R源代码获取有关使用此值的准则时，我会看到多种不同的模式.

I often end up in situations where it is necessary to check if the obtained difference is above machine precision. Seems like for this purpose R has a handy variable: .Machine$double.eps. However when I turn to R source code for guidelines about using this value I see multiple different patterns.

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

Here are a few examples from stats library:

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调整由于精度问题引起的差异的准则?

How can one understand the reasoning behind all those different 10 *, 100 *, 50 * and sqrt() modifiers?
Are there guidelines about using .Machine$double.eps for adjusting differences due to precision issues?

答

double的机器精度取决于其当前值. .Machine$double.eps给出值为1时的精度.可以使用C函数nextAfter获取其他值的机器精度.

The machine precision for double depends on its current value. .Machine$double.eps gives the precision when the values is 1. You can use the C function nextAfter to get the machine precision for other values.

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.检查差异是否小于机器精度.用<完成.修改者可能会考虑典型情况，添加项多久未显示更改.

Adding value a to value b will not change b when a is <= half of it's machine precision. Checking if the difference is smaler than machine precision is done with <. The modifiers might consider typical cases how often an addition did not show a change.

在 R 中，可以通过以下方式估算机器精度:

In R the machine precision can be estimated with:

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值表示一个范围.对于简单的加法，结果的范围取决于每个被加数的范围以及它们的和的范围.

Each double value is representing a range. For a simple addition, the range of the result depends on the reange of each summand and also the range of their sum.

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.

For higher precission Rmpfr could be used.

library(Rmpfr)
mpfr("2", 1024L)^54 - 2.9 + 4.1 - (mpfr("2", 1024L)^54 + 5.9)
#[1] -4.700000000000000621724893790087662637233734130859375

如果可以将其转换为整数，则可以使用gmp(Rmpfr中的内容).

In case it could be converted to integer gmp could be used (what is in Rmpfr).

library(gmp)
as.bigz("2")^54 * 10 - 29 + 41 - (as.bigz("2")^54 * 10 + 59)
#[1] -47

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

相关推荐