R中分数函数中的数值爆炸问题答案

【问题标题】：Numerical blowup problem in a fractional function in RR中分数函数中的数值爆炸问题
【发布时间】：2019-07-02 16:18:37
【问题描述】：

巧，

我在 R 中使用这个函数：

betaFun = function(x){
  if(x == 0){
    return(0.5)
  }
  return( ( 1+exp(x)*(x-1) )/( x*(exp(x)-1) ) )
}

函数对于每个 x 都是平滑且定义明确的（至少从理论的角度来看），并且在 0 中极限接近 0.5（您可以通过使用 Hopital 定理来说服自己）。

我有以下问题：

即事实上，由于限制，R 错误地计算了值，我得到了 0 的爆炸。

这里我报告数值问题：

x = c(1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-13)  
sapply(x, betaFun)

[1] 5.000083e-01 5.000442e-01 2.220446e+00 0.000000e+00 0.000000e+00 1.111111e+10

如您所见，评估非常奇怪，尤其是最后一个。我以为我可以通过在 0 中定义缺失值来解决这个问题（从代码中可以看出），但事实并非如此。

你知道我该如何解决这个数值爆炸问题吗？

这个函数需要高精度，因为我必须将它反转到 0 左右。我将使用 nleqslv 库中的 nleqslv 函数来实现。当然，如果函数有数值问题，反演会返回错误的解。

【问题讨论】：

curve((1+exp(x)*(x-1) )/( x*(exp(x)-1)), -10, 10, n = 1e5) 看不到这一点。请提供一个完整的可重现示例。
@Roland 请注意，您选择的点最多只能达到 1E-4。 OP 计算 1E-13
能否请您使绘图与使用的公式兼容？对于所描述的域[-10,10]，给定的函数并没有那么接近渐近线，渐近线在y=0 和y=1，并且奇点在x=0 附近更加紧密，其中值为0.5，不是0。

标签： r methods numeric numerical-methods

【解决方案1】：

你的问题是你取两个绝对值非常小的数字的商。此类数字仅表示为浮点精度。

对于接近于零的 x 值，您没有说明为什么需要这些函数值。一种简单的选择是强制转换为高精度数字：

library(Rmpfr)  
betaFun = function(x){
  x <- mpfr(as.character(x), precBits = 256) 
  #if x is calculated, you should switch to high precision numbers for its calculation
  #this step could be removed then

  #do calculation with high precision, 
  #then coerce to normal precision (assuming that is necessary)
  ifelse(x == 0, 0.5, as((1 + exp(x) * (x - 1)) / (x * (exp(x) - 1)), "numeric"))
}  

x = c(1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-13, 0) 
betaFun(x)
#[1] 0.5000083 0.5000001 0.5000000 0.5000000 0.5000000 0.5000000 0.5000000

【讨论】：

Thx，除了 x = 0 外，一切正常。添加一个 if 语句在 x = 0 处手动赋值就足够了。
酷，它保持矢量结构！

【解决方案2】：

我认为当 x 接近 0 时，您在评估 exp(x)-1 时会失去准确性。在 C 中，如果我将您的函数评估为

double  f2( double x)
{   return (x==0)   ? 0.5
            : (x*exp(x) - expm1(x))/( x*expm1(x));
}

问题消失了。这里 expm1 是一个数学库函数，它计算 exp(x) - 1，而不会损失小 x 的准确性。恐怕我不知道 R 是否有这个，但你希望它会。

不过，我认为您最好测试一下 |x|。足够小，而不是 0.0。关键是，对于足够小的 x，x*exp(x) 和 expm1(x) 都将是双精度 x，因此它们的差将为 0。为了保持最大精度，可能需要在 0.5 中添加一个线性项返回。我还没有准确计算出“足够小”应该是多少，但我认为它在 1e-16 左右。

【讨论】：

是的，R中有一个expm1函数。

【解决方案3】：

如您所见，您遇到的问题接近为零。分子和分母的根都为零。正如 OP 所提到的，使用 L'Hôpitcal，您会注意到 f(x) = 1/2。

从数字的角度来看，情况略有不同。浮点数总是会出错，因为并非每个实数都可以表示为浮点数。例如：

exp(1E-3)  -1 = 0.0010005001667083845973138522822409868               # numeric
exp(1/1000)-1 = 0.001000500166708341668055753993058311563076200580... # true
                                  ^

数值计算exp(1E-3)-1的问题已经从开头开始，即1E-3

1E-3 = x   = 0.0010000000000000000208166817117216851
exp(x)     = 1.0010005001667083845973138522822409868
exp(x) - 1 = 0.0010005001667083845973138522822409868

1E-3 不能表示为浮点数，精确到 17 位。
IEEE 将给出可能与 x 的真实值最接近的浮点值，由于 (1) 已经存在错误。仍然exp(x) 最多只能精确到 17 位。
通过减 1，我们得到一串零，现在我们的结果只能精确到 14 位。

既然我们知道我们不能将所有内容都表示为浮点数，您应该意识到接近零会变得有点尴尬，分子和分母都变得越来越不准确，尤其是在 1E-13 附近。

numerator_numeric(1E-13) = 1.1102230246251565E-16
numerator_true(1E-13)    = 5.00000000000033333333333...E-27

一般来说，你在这个点附近做的是在零附近使用泰勒展开，在其他任何地方都使用正常函数：

betaFun = function(x){
  if(-1E-1 < x && x < 1E-1){
    return(0.5 + x/12. - x^3/720. + x^5/30240.)
  }
  return( ( 1+exp(x)*(x-1) )/( x*(exp(x)-1) ) )
}

上述展开对于小区域中的 x 精确到 13 位

【讨论】：