r - R 中分数函数的数值爆炸问题

Question

再见，

我正在 R 中使用这个函数:

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

该函数对于每个 x 都是平滑且定义良好的(至少从理论角度来看)，并且在 0 中极限接近 0.5(您可以通过使用霍皮塔定理来说服自己这一点)。

我有以下问题:

即事实上，由于限制，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 来解决这个问题(正如您从代码中看到的那样)，但事实并非如此。

Do you know how can I solve this numerical blow up problem?

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

Answer 1

我认为你在 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) 都将是 double x，因此它们的差值为 0。为了保持最大精度，可能需要在 0.5 的基础上添加一个线性项返回。我还没有准确地计算出“足够小”应该是多少，但我认为大约在 1e-16 左右。