0

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

Thanks

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?

user6376297
  • 791

3 Answers3

2

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

Claude Leibovici
  • 289,558
  • 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? – user6376297 Mar 10 '21 at 13:16
  • @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 :-) – Claude Leibovici Mar 11 '21 at 02:21
  • 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. https://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. – user6376297 Mar 11 '21 at 10:05
  • OK I found two posts that maybe clarify this. https://math.stackexchange.com/questions/595426/series-expansion-of-a-function-at-infinity/595437 and https://math.stackexchange.com/questions/51770/taylor-expansion-at-infinity . 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. – user6376297 Mar 11 '21 at 10:14
  • 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. – user6376297 Mar 11 '21 at 10:45
  • @user6376297. Where is your problem ? In the derivative ? Do you want to go to a chat room ? If yes, open a room and I join you. – Claude Leibovici Mar 11 '21 at 11:00
  • @user6376297. I you prefer "Laurent series" – Claude Leibovici Mar 11 '21 at 11:07
  • Indeed, I submitted 'Taylor series exp(-x^2)/(1+erf(x)) of degree 6 for x -> -infinity' to Wolfram Alpha and it gave me the same result as yours, but with a '(Laurent series)' note next to it. I had found 'Laurent series' when searching for Taylor series with negative powers of the variable, but when I saw those circular integrals and it said that it applied to complex functions, I thought it had little to do with my example. This post says something about Laurent series. Still very unclear to me. – user6376297 Mar 11 '21 at 11:37
  • I created the chat room. Not sure why it assigned it to chemistry, but OK. https://chat.stackexchange.com/rooms/120743/laurent-series-workshop – user6376297 Mar 11 '21 at 11:44
1

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

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

njuffa
  • 2,168
  • 1
  • 20
  • 21