为什么这些数字不相等?

7 浏览
0 Comments

为什么这些数字不相等?

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

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

admin 更改状态以发布 2023年5月24日
0
0 Comments

在 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

0
0 Comments

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

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

这是计算机算术的一个众所周知的局限性,在几个地方都有讨论:

比较标量

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
 #....
 #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 采用的就是这种方法,该方法自我说明为

这是一种安全的方法,用于比较两个浮点数向量是否(成对地)相等。这比使用 == 更安全,因为它内置了容差

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

0