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（您可以使用 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 中定义缺失值来解决这个问题（正如您从代码中看到的那样），但事实并非如此。

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

我需要这个函数的高精度，因为我必须在 0 左右反转它。我将使用来自 nleqslv 的 nleqslv 函数来完成它图书馆。当然，如果函数有数值问题，反转将 return 错误的解决方案。

Answer 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

Answer 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 上添加一个线性项return。我还没有准确计算出“足够小”应该是多少，但我认为它在 1e-16 左右。

Answer 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 位数字

R中分数函数的数值放大问题

Numerical blowup problem in a fractional function in R

methods

r

numeric

numerical-methods