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

编程问答 2022-04-23

问题描述

我正在解决一些必须在标准空间和对数空间之间转换概率的编程问题。为此，我试图找出 R 中浮点误差的最大绝对误差，用于计算 log(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.eps（即 2.220446e-16）大小的一半。我不确定这是否存在理论上的原因，或者我是否得到了错误的答案。

问题：有什么方法可以确定这是否真的是此计算的最大浮点误差？有没有理论上的方法来计算最大值，或者这种网格搜索方法是否足够？

解决方法

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

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

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

但是，一旦您的 x 非常大，您就会出现 Inf 错误，例如，

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

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

log(exp(x)) 的误差取决于 x 的值。如果你使用 float ，那么 x 也有一个取决于它的值的精度。可以使用 nextafter 中的 C 计算精度：

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
}

也看看：What is the correct/standard way to check if difference is smaller than machine precision?。

,

总的来说，我建议使用随机方法来生成更多 x，例如：

x <- runif(10000000,2)

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

这也取决于你是关心绝对误差还是relative error。绝对误差应接近 .Machine$double.xmax，而相对误差随着 x 接近零而增加。例如。 log(exp(1e-16)) 被截断为零。

floating-accuracyfloating-pointprecisionrr