R:计算 log 的最大浮点误差(exp(...))

R: Computing maximum floating-point error for log(exp(...))

提问人:Ben 提问时间:6/21/2021 最后编辑:Ben 更新时间:6/22/2021 访问量:280

问:

我正在研究一些编程问题,我必须在标准空间和对数空间之间转换概率。为此,我试图找出输入为对数概率(即非正数)的计算中浮点误差的最大绝对误差Rlog(exp(...))

目前,我已经使用网格搜索计算了答案(请参阅下面的代码和绘图),但我不确定我计算的值是否正确。(我检查了其他一些范围,但图中显示的范围似乎获得了最大绝对误差。

#Set function for computing floating-point error of log(exp(...))
fp.error <- function(x) { abs(log(exp(x)) - x) }

#Compute and plot floating-point error over a grid of non-negative values
xx <- -(0:20000/10000)
ff <- fp.error(xx)
plot(xx, ff, col = '#0000FF10',
     main = 'Error in computation of log(exp(...))', 
     xlab = 'x', ylab = 'Floating-Point Error')

#Compute maximum floating-point error
fp.error.max <- max(ff)
fp.error.max
[1] 1.110223e-16

enter image description here

根据这个分析,我估计的最大绝对误差是 (即 ) 的一半大小。我不确定这是否有理论上的原因,或者我是否得到了错误的答案。.Machine$double.eps2.220446e-16

问题:有没有办法确定这是否真的是此计算的最大浮点误差?有没有理论上的方法可以计算最大值,或者这种网格搜索方法就足够了吗?

r 浮点 精度 float-accuracy

评论

0赞 chux - Reinstate Monica 6/21/2021
你只关心 [0...2.0] 吗?较大的值呢?
0赞 Ben 6/22/2021
@chux:编辑---我检查了其他一些范围,但显示的范围似乎得到了最大的绝对误差。
0赞 Sam Mason 6/22/2021
如果你正在使用对数概率/似然,那么如果可能的话,最好不要使用幂。有一些函数允许你尽可能地留在对数空间中,另请参阅 和 有关更稳定的函数logsumexplog1pexpm1

答:

2赞 ThomasIsCoding 6/21/2021 #1

我想你得到了正确的答案。在这里,我细化了小到,你会看到sqrt(.Machine$double.eps)

> x <- seq(0, 2, by = sqrt(.Machine$double.eps))

> max(abs(log(exp(x)) - x))
[1] 1.110725e-16

但是,一旦您的非常大,您就会遇到错误,例如,xInf

> (x <- .Machine$double.xmax)
[1] 1.797693e+308

> max(abs(log(exp(x)) - x))
[1] Inf

评论

0赞 Sam Mason 6/21/2021
请注意,您将获得任何 ,即大约 710Infx > log(.Machine$double.xmax)
0赞 ThomasIsCoding 6/21/2021
@SamMason是的。实际上源于 ,然后依次给出。我只想向 OP 显示可能的错误值。Infexplog(exp(x))Inf
0赞 Sam Mason 6/21/2021
注意,我的意思是建议使用来恢复会导致溢出的值logexp
0赞 ThomasIsCoding 6/21/2021
@SamMason OP 似乎不是在寻找解决方法,而是在寻找可能的错误。如果 OP 想要将错误风险降至最低,则应改用。log(exp(...))exp(log(x))
0赞 Ben 6/22/2021
@ThomasIsCoding:很抱歉,---我刚刚更新了我的问题。我实际上正在处理对数概率的输入,因此它们是非正数(这可能会对您的答案产生影响)。
2赞 GKi 6/21/2021 #2

的误差取决于 的值。如果使用浮点数,则 x 的精度也取决于其值。精度可以通过以下公式计算:log(exp(x))xnextafterC

library(Rcpp)
cppFunction("double getPrec(double x) {
  return nextafter(x, std::numeric_limits<double>::infinity()) - x;}")

getPrec(2)
#[1] 4.440892e-16

getPrec(exp(2))
#[1] 8.881784e-16

或不使用 :Rcpp

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
}

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

评论

0赞 Ben 6/22/2021
很抱歉,在我刚刚更新了我的问题---很痛苦。我实际上正在处理对数概率的输入,因此它们是非正的。
1赞 GKi 6/22/2021
对于这些值,最大误差将在以下范围内: 给出 so 值 从 to 将被存储为 并给出 给出 to 的范围。所以你计算出的误差说 a 是 a,但 a 存储为浮点数的范围可以从 所以我会说最大可能的误差在仅考虑 和 的精度时的范围内。getPrec(0)4.940656e-3240 - 4.940656e-324/20 + 4.940656e-324/20getPrec(exp(0))2.220446e-161 - 2.220446e-16/21 + 2.220446e-16/211.110223e-161111 +- 1.110223e-164.940656e-324 + 2.220446e-1601
2赞 Sam Mason 6/21/2021 #3

一般来说,我建议使用随机方法来生成更多的 s,例如:x

x <- runif(10000000, 0, 2)

您的规则间隔值可能会偶然发现一个“有效”的模式。

此外,这取决于您是关心绝对误差还是相对误差。绝对误差应接近,而相对误差随着接近零而增加。例如 被截断为零。.Machine$double.xmaxxlog(exp(1e-16))