We need to compute this function:

$f(x) = x + \sqrt {\frac 2 \pi} \cdot \sigma_x \cdot \frac {e^{- ( \frac {x} {\sqrt 2 \cdot \sigma_x} )^2}} {1 + erf(\frac {x} {\sqrt 2 \cdot \sigma_x})}$

We noticed that, for negative values of $\frac {x} {\sqrt 2 \cdot \sigma_x}$ lower than a certain value (and different depending on which software we use to compute it, say about $-5$ for instance), there is numerical instability (oscillations).

I saw this post suggesting how to achieve better stability, although on a rather different function:

Numerical Instability with this cumulate

Do you think a similar method would work here? And if so, should I express the $erf$ as an integral?
Or would you suggest some other method? Taylor expansion? Anything else...?


EDIT after evaluation of the responses

Thank you all very much for your feedback!

I think there is a complication here that I had not explicitly foreseen.

First I will make the following substitution, so your alternative expressions of $erf$ can be used directly:

$x = y \cdot \sqrt 2 \cdot \sigma_x$

$f(y) = \sqrt 2 \cdot \sigma_x \cdot y + \sqrt {\frac 2 \pi} \cdot \sigma_x \cdot \frac {e^{- y^2}} {1 + erf(y)}$

Now, the theory we used to derive this equation requires that:

$f(y) > 0, \forall y$
$\lim_{y \to -\infty} f(y) = 0$

So, we expect $f(y)$ to be very slowly approaching $0$ as $y$ becomes more negative; but $f(y)$ should never be negative.

In fact, if I plot the second term of $f(y)$, it looks like, as $y$ becomes more and more negative:

$\sqrt {\frac 2 \pi} \cdot \sigma_x \cdot \frac {e^{- y^2}} {1 + erf(y)} \approx - \sqrt 2 \cdot \sigma_x \cdot y$

This is in substantial agreement with Claude's approximation for the $erf$ part.
I.e., using:

$\frac{e^{-y^2}}{1+\text{erf}(y)} \approx -\sqrt{\pi } \,y\Bigg[1+\frac{1}{2 y^2}-\frac{1}{2 y^4}+\frac{5}{4 y^6}+O\left(\frac{1}{y^8}\right)\Bigg]$

I get:

$\sqrt {\frac 2 \pi} \cdot \sigma_x \cdot \frac {e^{- y^2}} {1 + erf(y)} \approx \sqrt {\frac 2 \pi} \cdot \sigma_x \cdot (-\sqrt{\pi }) \,y\Bigg[1+\frac{1}{2 y^2}-\frac{1}{2 y^4}+\frac{5}{4 y^6}+O\left(\frac{1}{y^8}\right)\Bigg] =$
$= - \sqrt {2} \cdot \sigma_x \cdot \,y\Bigg[1+\frac{1}{2 y^2}-\frac{1}{2 y^4}+\frac{5}{4 y^6}+O\left(\frac{1}{y^8}\right)\Bigg]$

which indeed approaches $- \sqrt 2 \cdot \sigma_x \cdot y$ when $y$ becomes very negative.

With the above approximation:

$f(y) \approx \sqrt 2 \cdot \sigma_x \cdot y - \sqrt {2} \cdot \sigma_x \cdot \,y\Bigg[1+\frac{1}{2 y^2}-\frac{1}{2 y^4}+\frac{5}{4 y^6}+O\left(\frac{1}{y^8}\right)\Bigg] = $ $= - \sqrt {2} \cdot \sigma_x \cdot \,y\Bigg[\frac{1}{2 y^2}-\frac{1}{2 y^4}+\frac{5}{4 y^6}+O\left(\frac{1}{y^8}\right)\Bigg] $

Numerically, this 'does the trick', i.e. for $y \le -4$ it seems really very close to the function I need to compute, as indicated by Claude. And in fact with the rational function version even up to $y \le -1$.

This seems to imply that only polynomial approximations are viable, not the $erfc$-based ones, or am I wrong?

Any thoughts?


I would suggest a series expansion for $x \leq- 4$ $$f(x)=\frac{e^{-x^2}}{1+\text{erf}(x)}=-\sqrt{\pi } \,x\Bigg[1+\frac{1}{2 x^2}-\frac{1}{2 x^4}+\frac{5}{4 x^6}+O\left(\frac{1}{x^8}\right)\Bigg]$$ For $x=-5$, the "exact" value is $9.03318$ while this truncated expansion gives $9.03313$.

Much better would be the Padé approximants $$f(x)=-\sqrt{\pi } \,x \,\frac{8 x^6+84 x^4+210 x^2+105 } {8 x^6+80 x^4+174 x^2+48 }$$ which, for $x=5$, leads to an error of $2.2\times 10^{-8}$.

  • $\begingroup$ Thank you Claude, your solution seems to work very well! I am probably going to accept it as the answer. Could you please just comment briefly on how you obtained that expansion? $\endgroup$ – user6376297 Mar 10 at 13:16
  • $\begingroup$ @user6376297. The first one is obtained by composition of the Taylor series for large negative values of $x$. Taking more terms, it is easy to build the corresponding Padé approximant. The one I produced is equivalent to a Taylor series to $O\left(\frac{1}{x^{13}}\right)$. If you want better, let me know since it just need a couple of minutes to build. Cheers :-) $\endgroup$ – Claude Leibovici Mar 11 at 2:21
  • $\begingroup$ Thanks, no, I do not need a more precise approximation, I was just trying to understand how you made yours, so for similar cases I could do it myself. Unfortunately I cannot see how that is a Taylor series, as per definition it should only have positive powers of x. en.wikipedia.org/wiki/Taylor_series . Also, I would not know on which point to 'center' the series, as here we say 'very negative x', which for me means make x go to $-\infty$, but then I cannot write Taylor terms with $(x+\infty)^n$... So I am a bit stuck. I will do more research. And I will accept your answer. $\endgroup$ – user6376297 Mar 11 at 10:05
  • $\begingroup$ OK I found two posts that maybe clarify this. math.stackexchange.com/questions/595426/… and math.stackexchange.com/questions/51770/… . Perhaps there are some substitution and back-substitution steps being taken, not simply a Taylor expansion. The other issue is that the CAS I am using to do this gets the limit wrong for $f(x)$ and its derivatives at $x \to - \infty$. But OK, that can be addressed somehow. $\endgroup$ – user6376297 Mar 11 at 10:14
  • $\begingroup$ No, tried it out and still does not work for me :( I might post another question about that, unless you can comment on the method you used. $\endgroup$ – user6376297 Mar 11 at 10:45

In the negatives, $$\text{erf(x)}$$ very quickly tends to $-1$ and it is no surprise that you soon reach the limits of the floating-point representation.

You will regain full accuracy by computing

$$\text{erf}(x)+1$$ directly. It is possible that you will obtain better results with $\text{erfc}(-x)$, but that depends on the numerical library you use.

Otherwise, you could use the fourth approximation in https://en.wikipedia.org/wiki/Error_function#Approximation_with_elementary_functions, and you will see a nice simplification with your numerator.


$1 + erf(x)$ suffers from catastrophic cancellation in the negative half plane. This is avoided by replacing it with $erfc(-x)$, that is, use of the complementary error function. Most standard math libraries, for example the C++ standard math library, include this as a function erfc.

Many special function libraries also offer the exponentially scaled complementary error function, $e^{x^{2}} erfc(x)$, often as a function called erfcx. With that, one can compute $\frac{e^{-x^{2}}}{1+erf(x)}$ accurately as 1.0 / erfcx (-x). I also provided a double-precision implementation of erfcx on Stackoverflow.

As I became aware belatedly, there is a second instance of catastrophic cancellation in $f(x)$ for negative $x$ large in magnitude: $x$ and the rest of the expression are similar in magnitude but of opposite sign. Around $-10^{7}$ complete loss of accuracy occurs.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.