Problème d'explosion numérique dans une fonction fractionnaire de R

Votre problème est que vous prenez le quotient de deux nombres avec de très petites valeurs absolues. Ces nombres ne sont représentés qu'avec une précision en virgule flottante.

Vous ne spécifiez pas pourquoi vous avez besoin de ces valeurs de fonction pour des valeurs x proches de zéro. Une option simple serait de forcer des nombres de haute précision:

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

Je pense que vous perdez de la précision dans l'évaluation de exp (x) -1 pour x proche de 0. En C si j'évalue votre fonction comme

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

Le problème disparaît. Ici expm1 est une fonction de bibliothèque mathématique qui calcule exp (x) - 1, sans perdre en précision pour les petits x. J'ai peur de ne pas savoir si R a ceci, mais vous espérez que ce serait le cas.

Je pense cependant que vous feriez mieux de tester | x | était suffisamment petit, plutôt que 0,0. Le fait est que pour un x assez petit, x * exp (x) et expm1 (x) seront, en tant que doubles, x, donc leur différence sera de 0. Pour conserver une précision maximale, vous devrez peut-être ajouter un terme linéaire à 0,5 vous revenir. Je n'ai pas déterminé précisément ce que «suffisamment petit devrait être, mais c'est quelque part autour de 1e-16 je pense.

Comme vous le constatez, vous rencontrez le problème proche de zéro. Les racines du numérateur et du dénominateur sont nulles. Et comme l'OP l'a mentionné, en utilisant L'Hépitcal, vous remarquez qu'en cela f (x) = 1/2 .

D'un point de vue numérique, les choses vont légèrement différemment. Les virgules flottantes auront toujours une erreur car tous les nombres réels ne peuvent pas être représentés sous forme de nombres à virgule flottante. Par exemple:

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) ) )
}

Le problème de l'évaluation numérique de exp (1E-3) -1 commence déjà au début, c'est-à-dire 1E-3 p>

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

1. 1E-3 ne peut pas être représenté sous forme de virgule flottante et est précis jusqu'à 17 chiffres.
2. IEEE donnera la valeur en virgule flottante la plus proche possible de la valeur vraie de x, qui a déjà une erreur due à (1). Encore exp (x) n'est précis que jusqu'à 17 chiffres.
3. En soustrayant 1, nous obtenons un tas de zéros au début, et maintenant notre résultat n'est précis que jusqu'à 14 chiffres.

Alors maintenant que nous savons que nous ne pouvons pas tout représenter exactement comme une virgule flottante, vous devez vous rendre compte que près de zéro, cela devient un peu gênant et le numérateur et le dénominateur deviennent de moins en moins précis, surtout près de 1E- 13.

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

Généralement, ce que vous faites près d'un tel point est d'utiliser une expansion de Taylor autour de zéro, et la fonction normale partout ailleurs:

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