Faster / more precise logarithm of complimentary error function?

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.

Answer 1

[Walter Roberson]

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

http://www.math.kobe-u.ac.jp/icms2006/icms2006-video/slides/115.pdf

https://hal.inria.fr/inria-00261360/document

[Gleb Tikhonov]

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?

Answer 2

[John D'Errico]

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.

Answer 3

[ES]

>> format long
>> log(erfc(-5:5))

Answer 4

[Shep Bryan]

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

Answer 5

[Luca Citi]

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.