Faster / more precise logarithm of complimentary error function?

20 次查看(过去 30 天)
If there is a way to calculate the logarithm of complimentary error function with better numerical precision (and preferably also faster) than doing it in the straightforward manner, e.g.
log(erfc(-5:5))
ans =
0.6931 0.6931 0.6931 0.6908 0.6112 0 -1.8496 -5.3649 -10.7204 -17.9878 -27.2009
As far as I know, many of numerical approximations involve product that includes exponential in some form e.g. see Wikipedia, so it feels that some potential gain could be achieved through it.

回答(5 个)

Walter Roberson
Walter Roberson 2017-9-23
编辑:Walter Roberson 2017-9-23
No, the algorithms involve the sum of exponentiation times another factor, not products that you could potentially turn into sums.
erfc() is a built-in in MATLAB. It is going to call into some built-in library function written in C or C++ . To have any chance of matching the speed, your replacement would also have to be in C or C++.
  1 个评论
Gleb Tikhonov
Gleb Tikhonov 2017-9-23
Thanks for the links. However your answer seems to be somewhat a bit contradictory - e.g. in the slide 31 the formulas (4) and (5) look like they could easily be log transformed into sums. It is also super clear that the build-in erfc should be faster than whatever mine MATLAB written approximation for it, but here the specifics is that I do not need the erfc itself, but instead the logarithm of it. Similarly I am concerned about the precision of log(erfc(x)), not the precision of the erfc(x).
What I am looking for is some alternative to following R code, which allows to compute the logarithm of Gaussian CDF, which then could be easily adjusted to get the erfc value:
log(2) + pnorm(-sqrt(2)*c(-5:5), log=TRUE)
[1] 0.6931472 0.6931472 0.6931361 0.6908056 0.6112323 0.0000000 -1.8496055 -5.3649413 -10.7203630 -17.9877783 -27.2008895
It is also important to mention that whatever algorithm is implemented there, it can calculate the logarithm of CDF even when the numerical precision is not enough for the CDF value itself:
log(pnorm(-100))
[1] -Inf
pnorm(-100, log=TRUE)
[1] -5005.524
Any ideas on whether I can get something similar in MATLAB without recoding the approximation algorithm myself?

请先登录,再进行评论。


John D'Errico
John D'Errico 2017-9-23
If you want speed, no, you are not going to beat log(erfc(x)), at least not without writing custom code, preferably in some compiled language.
If you want higher accuracy for large positive x, then yes you could get more accuracy in the upper tail. It won't be incredibly fast though. I'd probably start with the simple asymptotic expansion given on the wikipedia page.
Be careful though, as the classic asymptotic expansion actually diverges if you take too many terms. But for moderately small x, you would need perhaps more terms to get an adequate number of digits. It depends on how many digits of accuracy you need/want. The point is, you need to do some careful work to choose an expansion that will be fast AND highly accurate. The problem will be for intermediate values of x, perhaps on the order of 1 to 2.

ES
ES 2017-9-23
>> format long
>> log(erfc(-5:5))

Shep Bryan
Shep Bryan 2020-3-31
If you are mainly worried about underlow, this is a handy function that avoids numerical issues:
function out = logerfc(X)
out = log(erfcx(X)) - X.^2;
end

Luca Citi
Luca Citi 2023-5-14
编辑:Luca Citi 2023-5-14
Not much faster but certainly more accurate when erfc(x)⩬1:
log1p(erf(-x))
So switching between log1p(erf(-x)) and log(erfc(x)) depending on the argument shoud give good accuracy everywhere.

类别

Help CenterFile Exchange 中查找有关 Interpolation of 2-D Selections in 3-D Grids 的更多信息

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by