R:计算 log(exp(...)) 的最大浮点误差
R: Computing maximum floating-point error for log(exp(...))
我正在处理一些编程问题,我必须在标准 space 和 log-space 之间转换概率。为此,我试图找出 最大绝对误差 中的浮点误差 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
根据这个分析,我估计的最大绝对误差是 .Machine$double.eps
(即 2.220446e-16
)的一半。我不确定这是否有理论上的原因,或者我是否得到了错误的答案。
问题:有什么方法可以确定这是否真的是该计算的最大浮点误差?有没有理论上的方法可以计算最大值,或者这种网格搜索方法是否足够?
我认为您的答案是正确的。这里我把step细化到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
}
另请参阅:。
一般来说,我建议使用随机方法生成更多 x
s,例如:
x <- runif(10000000, 0, 2)
您的固定间隔值可能会偶然发现“有效”的模式。
也看你是关心绝对误差还是relative error。绝对误差应接近 .Machine$double.xmax
,而相对误差随着 x
接近零而增加。例如。 log(exp(1e-16))
被截断为零。
我正在处理一些编程问题,我必须在标准 space 和 log-space 之间转换概率。为此,我试图找出 最大绝对误差 中的浮点误差 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
根据这个分析,我估计的最大绝对误差是 .Machine$double.eps
(即 2.220446e-16
)的一半。我不确定这是否有理论上的原因,或者我是否得到了错误的答案。
问题:有什么方法可以确定这是否真的是该计算的最大浮点误差?有没有理论上的方法可以计算最大值,或者这种网格搜索方法是否足够?
我认为您的答案是正确的。这里我把step细化到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
}
另请参阅:
一般来说,我建议使用随机方法生成更多 x
s,例如:
x <- runif(10000000, 0, 2)
您的固定间隔值可能会偶然发现“有效”的模式。
也看你是关心绝对误差还是relative error。绝对误差应接近 .Machine$double.xmax
,而相对误差随着 x
接近零而增加。例如。 log(exp(1e-16))
被截断为零。