Численная задача о разрушении дробной функции в R - PullRequest
Новая
       121

Численная задача о разрушении дробной функции в R

0 голосов
/ 08 февраля 2019

Ciao,

Я работаю с этой функцией в R: 

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

Функция гладкая и хорошо определена для каждого x (по крайней мере, с теоретической точки зрения) ив 0 предел приближается к 0,5 (вы можете убедиться в этом, используя теорему Хопиталя).

У меня есть следующая проблема: enter image description here

т.е.что из-за ограничения 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 .Конечно, инверсия вернет неправильные решения, если функция имеет числовые проблемы.

Ответы [ 3 ]

0 голосов
/ 08 февраля 2019

Я думаю, что вы теряете точность в оценке exp (x) -1 для x, близкого к 0. В 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вернуть.Я не совсем понял, что должно быть достаточно маленьким, но я думаю, что это где-то около 1е-16.

0 голосов
/ 17 февраля 2019

Как вы заметили, вы столкнулись с проблемой около нуля.Корни как числителя, так и знаменателя равны нулю.И, как упоминалось в ОП, используя 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
  1. 1E-3 не может быть представлена ​​какс плавающей запятой, с точностью до 17 цифр.
  2. IEEE даст наиболее близкое возможное значение с плавающей запятой к истинному значению x, которое уже имеет ошибку из-за (1).Тем не менее, exp(x) имеет точность только до 17 цифр.
  3. Вычитая 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) ) )
}

Вышеупомянутое расширение с точностью до 13 цифр для х вмаленький регион

0 голосов
/ 08 февраля 2019

Ваша проблема в том, что вы берете частное от двух чисел с очень маленькими абсолютными значениями.Такие числа представлены только с точностью до плавающей запятой.

Вы не указываете, почему вам нужны эти значения функций для значений х, близких к нулю.Одним из простых вариантов будет приведение к высокоточным числам: 

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
...