I need to compute the following formula numerically using double precision:

$\frac {\frac{1}{2 \pi^{2} \sigma} \left(\pi^{2} \sigma \operatorname{erf}{\left (\frac{1}{\sqrt 2} \left(\frac{r}{\sigma} + \sigma\right) \right )} - \pi^{2} \sigma \operatorname{erf}{\left (\frac{1}{\sqrt 2} \left(\frac{r - b}{\sigma} + \sigma\right) \right )} + \sqrt{2} \pi^{\frac{3}{2}} \left(-1 + e^{- b}\right) e^{t}\right)} {\left(1 - \frac{1}{2} \operatorname{erf}{\left (\frac{r}{\sqrt 2 \sigma} \right )} - \frac{1}{2} \operatorname{erf}{\left (\frac{\left(b - r\right)}{\sqrt 2 \sigma} \right )}\right) e^{- r - \frac{\sigma^{2}}{2}} + \frac{1}{2} \operatorname{erf}{\left (\frac{\left(r + \sigma^{2}\right)}{\sqrt 2 \sigma} \right )} + \frac{1}{2} \operatorname{erf}{\left (\frac{\sqrt{2}}{2 \sigma} \left(b - r - \sigma^{2}\right) \right )}}$

It is guaranteed that:

  • $\sigma > 0$. This is the parameter I control. It would be nice if the computation worked for $\sigma$ from around 1e-3 to 1e5, though.
  • $r \in (0, 1e7]$.
  • $b \in (0, 1e7]$.
  • $t = - \frac {(b - r)^2} {2 \sigma^2} + (b - r) - \sigma^2/2 \leq 0$

The problem is that for some valid inputs, such as $\sigma=1$, $r=1111$, $b=200$, the formula results in a 0/0 division.

  • In the numerator, $\operatorname{erf}$s in the first and second terms get rounded to 1 and cancel out. The third term underflows.
  • In the denominator, the first term underflows, while the second and third ones cancel out.

How can I avoid these problems and compute this formula stably? Perhaps there are some identities for $\operatorname{erf}(a) - \operatorname{erf}(b)$ that will allow to avoid cancellation?

  • $\begingroup$ I still need to get the numerical answer, so the appropriate action would be to change the formula to work in this particular case, which is the point of my question. $\endgroup$ – Pastafarianist Sep 11 '17 at 15:10
  • $\begingroup$ I'm hoping that all pathological cases will be covered by a simple transformation of the formula, because both numerator and denominator underflow for the same reasons. That's also why I'm asking about identities for $\operatorname{erf}$. $\endgroup$ – Pastafarianist Sep 11 '17 at 15:20

When you have problems with cancellation, the usual prescription is to subtract off the large terms analytically leaving the small ones you care about. $\operatorname{erf}(x)$ gets very close to $1$ as $x$ gets even moderately large. In the numerator, your cancellation happens when $b$ is rather small and $\sigma$ is rather large. Then the term $-\frac b\sigma$ in the argument doesn't matter much. You can test for that and use the expansion $$\operatorname{erf}(x) \approx 1+e^{-x^2} \left(-\frac1{\sqrt\pi x} + \frac 1{2 \sqrt\pi x^3} - \frac3{4 \sqrt\pi x^5} + O(x^{-7})\right)$$ You probably don't need all the terms. The $1$s will cancel, leaving just the terms multiplied by $e^{-x^2}$. The problem is similar in the denominator, with the same solution.


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.