1

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?

Pastafarianist
  • 365
  • 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. – Pastafarianist Sep 11 '17 at 15:10
  • 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}$. – Pastafarianist Sep 11 '17 at 15:20

1 Answers1

1

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.

Ross Millikan
  • 374,822