这些数字为什么不相等?

310

以下代码显然是错误的,问题出在哪里?

i <- 0.1
i <- i + 0.05
i
## [1] 0.15
if(i==0.15) cat("i equals 0.15") else cat("i does not equal 0.15")
## i does not equal 0.15

7
请参阅https://dev59.com/wlnUa4cB1Zd3GeqPdb1m和https://dev59.com/4HE85IYBdhLWcg3wdDIM。 [R Inferno](http://www.burns-stat.com/pages/Tutor/R_inferno.pdf)也是另一个很棒的阅读材料。 - Aaron left Stack Overflow
1
一个站点范围内的与语言无关的问答:浮点数运算是否存在问题? - Gregor Thomas
dplanet,我在下面添加了一个双精度算术的解决方案,适用于所有比较情况(“<=”,“> =”,“=”)。希望能有所帮助。 - Erdogan CEVHER
6个回答

411

通用(与语言无关)的原因

由于不是所有的数字都可以在IEEE浮点运算(几乎所有计算机用来表示十进制数并进行数学运算的标准)中精确表示,所以你不会总是得到你期望的结果。这一点尤其真实,因为一些简单的有限小数(如0.1和0.05)在计算机中并不能精确表示,所以对它们进行算术运算的结果可能与“已知”答案的直接表示不完全相同。

这是计算机算术的一个众所周知的限制,并在多个地方进行了讨论:

R FAQ有一个专门讨论这个问题的问题:R FAQ 7.31 Patrick Burns的《R地狱》将第一个“圆圈”专门用于解决这个问题(从第9页开始) David Goldberg的“计算机科学家应该了解的浮点算术知识”,《ACM计算调查》23,1(1991-03),5-48 doi>10.1145/103162.103163也有修订版本《浮点数指南-每个程序员都应该了解的浮点算术知识》 0.30000000000000004.com比较了不同编程语言中的浮点算术 包括几个Stack Overflow问题 为什么浮点数不准确? 为什么十进制数不能在二进制中精确表示? 浮点数运算是否有问题? “浮点数不准确”的规范副本(关于此问题的规范答案的元讨论)

比较标量

R中,解决这个问题的标准方法不是使用==,而是使用all.equal函数。或者更准确地说,由于all.equal函数在存在差异时会提供很多细节,所以使用isTRUE(all.equal(...))

if(isTRUE(all.equal(i,0.15))) cat("i equals 0.15") else cat("i does not equal 0.15")

产量
i equals 0.15

一些使用all.equal而不是==的更多示例(最后一个示例旨在显示这将正确显示差异)。
0.1+0.05==0.15
#[1] FALSE
isTRUE(all.equal(0.1+0.05, 0.15))
#[1] TRUE
1-0.1-0.1-0.1==0.7
#[1] FALSE
isTRUE(all.equal(1-0.1-0.1-0.1, 0.7))
#[1] TRUE
0.3/0.1 == 3
#[1] FALSE
isTRUE(all.equal(0.3/0.1, 3))
#[1] TRUE
0.1+0.1==0.15
#[1] FALSE
isTRUE(all.equal(0.1+0.1, 0.15))
#[1] FALSE

一些更详细的内容,直接从一个类似问题的回答中复制过来:
你遇到的问题是浮点数在大多数情况下无法精确表示小数,这意味着你经常会发现精确匹配失败。
当你说 R 稍微有点不准确时:
1.1-0.2
#[1] 0.9
0.9
#[1] 0.9

你可以用十进制来了解它的真实想法。
sprintf("%.54f",1.1-0.2)
#[1] "0.900000000000000133226762955018784850835800170898437500"
sprintf("%.54f",0.9)
#[1] "0.900000000000000022204460492503130808472633361816406250"

你可以看到这些数字是不同的,但表示方式有点不方便。如果我们用二进制(或者十六进制,它们是等价的)来看,我们会得到一个更清晰的图像:
sprintf("%a",0.9)
#[1] "0x1.ccccccccccccdp-1"
sprintf("%a",1.1-0.2)
#[1] "0x1.ccccccccccccep-1"
sprintf("%a",1.1-0.2-0.9)
#[1] "0x1p-53"

你可以看到它们之间的差异为2^-53,这很重要,因为这个数字是两个接近1的数值之间最小可表示的差异。
我们可以通过查看R的machine字段来确定任何给定计算机的最小可表示数字。
 ?.Machine
 #....
 #double.eps     the smallest positive floating-point number x 
 #such that 1 + x != 1. It equals base^ulp.digits if either 
 #base is 2 or rounding is 0; otherwise, it is 
 #(base^ulp.digits) / 2. Normally 2.220446e-16.
 #....
 .Machine$double.eps
 #[1] 2.220446e-16
 sprintf("%a",.Machine$double.eps)
 #[1] "0x1p-52"

你可以利用这个事实来创建一个“几乎相等”的函数,它可以检查差异是否接近浮点数中最小可表示的数字。事实上,这个函数已经存在了:all.equal
?all.equal
#....
#all.equal(x,y) is a utility to compare R objects x and y testing ‘near equality’.
#....
#all.equal(target, current,
#      tolerance = .Machine$double.eps ^ 0.5,
#      scale = NULL, check.attributes = TRUE, ...)
#....

所以all.equal函数实际上是检查数字之间的差异是否是两个尾数之间最小差异的平方根。
这个算法在非常小的被称为denormals的数字附近有些奇怪,但你不需要担心这个。
比较向量
上面的讨论假设比较的是两个单个值。在R中,没有标量,只有向量,隐式向量化是该语言的一个优势。对于逐元素比较向量的值,前面的原则仍然适用,但实现略有不同。==是向量化的(进行逐元素比较),而all.equal将整个向量作为单个实体进行比较。
使用前面的例子
a <- c(0.1+0.05, 1-0.1-0.1-0.1, 0.3/0.1, 0.1+0.1)
b <- c(0.15,     0.7,           3,       0.15)
==不会给出“预期”的结果,而all.equal不会执行逐元素比较。
a==b
#[1] FALSE FALSE FALSE FALSE
all.equal(a,b)
#[1] "Mean relative difference: 0.01234568"
isTRUE(all.equal(a,b))
#[1] FALSE

相反,必须使用一个循环遍历两个向量的版本。
mapply(function(x, y) {isTRUE(all.equal(x, y))}, a, b)
#[1]  TRUE  TRUE  TRUE FALSE

如果需要一个功能版本,可以编写它。
elementwise.all.equal <- Vectorize(function(x, y) {isTRUE(all.equal(x, y))})

可以称之为“只是”。
elementwise.all.equal(a, b)
#[1]  TRUE  TRUE  TRUE FALSE

或者,不必在更多的函数调用中包装all.equal,你可以直接复制all.equal.numeric的相关内部,并使用隐式向量化:
tolerance = .Machine$double.eps^0.5
# this is the default tolerance used in all.equal,
# but you can pick a different tolerance to match your needs

abs(a - b) < tolerance
#[1]  TRUE  TRUE  TRUE FALSE

这是采取的方法,它自己的文档如下所示:
这是一种安全的比较两个浮点数向量是否(逐对)相等的方法。这比使用==更安全,因为它内置了容差。
dplyr::near(a, b)
#[1]  TRUE  TRUE  TRUE FALSE

测试向量中值的出现

如果应用于浮点数值,标准的R函数%in%也可能遇到相同的问题。例如:

x = seq(0.85, 0.95, 0.01)
# [1] 0.85 0.86 0.87 0.88 0.89 0.90 0.91 0.92 0.93 0.94 0.95
0.92 %in% x
# [1] FALSE

我们可以定义一个新的中缀运算符,以允许在比较中引入容差,具体如下:
`%.in%` = function(a, b, eps = sqrt(.Machine$double.eps)) {
  any(abs(b-a) <= eps)
}

0.92 %.in% x
# [1] TRUE

dplyr::near 包裹在 any 中也可以用于向量化检查

any(dplyr::near(0.92, x))
# [1] TRUE

46

加上Brian的评论(以下是原因),你可以通过使用all.equal来克服这个问题:

# i <- 0.1
# i <- i + 0.05
# i
#if(all.equal(i, .15)) cat("i equals 0.15\n") else cat("i does not equal 0.15\n")
#i equals 0.15

根据Joshua的警告,这里是更新后的代码(感谢Joshua):

 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

20
在使用 if 语句时,all.equal 函数如果存在差异,不会返回 FALSE,因此需要用 isTRUE 将其包装起来。 - Joshua Ulrich
为什么R没有更新==,以使dplanet代码能够正常工作?我也遇到了同样的问题。 - Yacine Hajji

14

这个方法有些hackish,但很快速:

if(round(i, 10)==0.15) cat("i equals 0.15") else cat("i does not equal 0.15")

3
你可以使用 all.equal(..., tolerance) 参数。all.equal(0.147, 0.15, tolerance=0.05) 为真。 - smci

14

dplyr::near()是测试两个浮点数向量是否相等的选项。这是来自文档的示例:

sqrt(2) ^ 2 == 2
#> [1] FALSE
library(dplyr)
near(sqrt(2) ^ 2, 2)
#> [1] TRUE

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


3

双精度浮点数中的泛化比较(“<=”,“>=”,“=”):

比较 a <= b:

IsSmallerOrEqual <- function(a,b) {   
# Control the existence of "Mean relative difference..." in all.equal; 
# if exists, it results in character, not logical:
if (   class(all.equal(a, b)) == "logical" && (a<b | all.equal(a, b))) { return(TRUE)
 } else if (a < b) { return(TRUE)
     } else { return(FALSE) }
}

IsSmallerOrEqual(abs(-2-(-2.2)), 0.2) # TRUE
IsSmallerOrEqual(abs(-2-(-2.2)), 0.3) # TRUE
IsSmallerOrEqual(abs(-2-(-2.2)), 0.1) # FALSE
IsSmallerOrEqual(3,3); IsSmallerOrEqual(3,4); IsSmallerOrEqual(4,3) 
# TRUE; TRUE; FALSE

比较a >= b:

IsBiggerOrEqual <- function(a,b) {
# Control the existence of "Mean relative difference..." in all.equal; 
# if exists, it results in character, not logical:
if (   class(all.equal(a, b)) == "logical" && (a>b | all.equal(a, b))) { return(TRUE)
 } else if (a > b) { return(TRUE)
     } else { return(FALSE) }
}
IsBiggerOrEqual(3,3); IsBiggerOrEqual(4,3); IsBiggerOrEqual(3,4) 
# TRUE; TRUE; FALSE

比较a = b:

IsEqual <- function(a,b) {
# Control the existence of "Mean relative difference..." in all.equal; 
# if exists, it results in character, not logical:
if (   class(all.equal(a, b)) == "logical" ) { return(TRUE)
 } else { return(FALSE) }
}

IsEqual(0.1+0.05,0.15) # TRUE

0

我曾经遇到过类似的问题,我使用了以下解决方案。

@我找到了一个关于不等间隔切割的解决方案。@我在R中使用了round函数。将选项设置为2位小数并不能解决问题。

options(digits = 2)
cbind(
  seq(      from = 1, to = 9, by = 1 ), 
  cut( seq( from = 1, to = 9, by = 1),          c( 0, 3, 6, 9 ) ),
  seq(      from = 0.1, to = 0.9, by = 0.1 ), 
  cut( seq( from = 0.1, to = 0.9, by = 0.1),    c( 0, 0.3, 0.6, 0.9 )),
  seq(      from = 0.01, to = 0.09, by = 0.01 ), 
  cut( seq( from = 0.01, to = 0.09, by = 0.01),    c( 0, 0.03, 0.06, 0.09 ))
)

根据选项(digits = 2),输出不等切割间隔的结果:
  [,1] [,2] [,3] [,4] [,5] [,6]
 [1,]    1    1  0.1    1 0.01    1
 [2,]    2    1  0.2    1 0.02    1
 [3,]    3    1  0.3    2 0.03    1
 [4,]    4    2  0.4    2 0.04    2
 [5,]    5    2  0.5    2 0.05    2
 [6,]    6    2  0.6    2 0.06    3
 [7,]    7    3  0.7    3 0.07    3
 [8,]    8    3  0.8    3 0.08    3
 [9,]    9    3  0.9    3 0.09    3


options(digits = 200)
cbind(
  seq(      from = 1, to = 9, by = 1 ), 
  cut( round(seq( from = 1, to = 9, by = 1), 2),          c( 0, 3, 6, 9 ) ),
  seq(      from = 0.1, to = 0.9, by = 0.1 ), 
  cut( round(seq( from = 0.1, to = 0.9, by = 0.1), 2),    c( 0, 0.3, 0.6, 0.9 )),
  seq(      from = 0.01, to = 0.09, by = 0.01 ), 
  cut( round(seq( from = 0.01, to = 0.09, by = 0.01), 2),    c( 0, 0.03, 0.06, 0.09 ))
)

根据round函数,返回等间隔切割的输出:
      [,1] [,2] [,3] [,4] [,5] [,6]
 [1,]    1    1  0.1    1 0.01    1
 [2,]    2    1  0.2    1 0.02    1
 [3,]    3    1  0.3    1 0.03    1
 [4,]    4    2  0.4    2 0.04    2
 [5,]    5    2  0.5    2 0.05    2
 [6,]    6    2  0.6    2 0.06    2
 [7,]    7    3  0.7    3 0.07    3
 [8,]    8    3  0.8    3 0.08    3
 [9,]    9    3  0.9    3 0.09    3

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接