Apologies if the question is too elementary here.

For a certain computational application I need to approximate Gaussian distribution $e^{-x^2}$ with specific absolute precision (within $10^{-7}$ over $\mathbb{R}$), preferably with rational functions.

Alas, I'm not familiar with approximation theory. Google pointed me toward Pade approximation as the way to go. Alas, I still don't know how to derive Pade approximation for a given function, much less how to ensure the approximation would fit to the prescribed precision. Could you point me towards the relevant information?


Here is some background. For some engineering computations I need analytical expressions for certain integrals; in this particular case for $e^{-x^2} \operatorname{erf}(x+a)$, where $a$ is a parameter that changes throughout the computation. Mathematica was unable to provide an analytical expression of $\int e^{-x^2}\operatorname{erf}(x+a)\,dx$ and neither was I.

That caused me thinking how can I approximate well known functions such as $e^{-x^2}$ and $\operatorname{erf}(x+a)$ (and possibly a few others I may need) so that their products would be analytically integrable?

It's very well known that one can obtain an analytic expression for an integral of a rational function, and that products of rational functions are rational. Thus the idea: approximate the relevant functions with rational ones using Pade approximation and integrate the products of Pade approximations. That would basically replace the integral I want with an expression containing a bunch of rational functions and $\log$s and $\arctan$s.

I've got to tightly control the precision of the above approximations though to have the model perform as expected. Thus the question.

  • $\begingroup$ Do you need to compute the Gaussian probability density or the comulative desnity (integral of the pdf from $-\infty$ to $x$) The former is trivial using an $\exp$ library function. The second is actually somewhat of a challenge. $\endgroup$ – Brian Borchers Jul 22 '16 at 0:13
  • $\begingroup$ @BrianBorchers: edited the question with the reason why I need that. $\endgroup$ – Michael Jul 22 '16 at 0:30
  • $\begingroup$ You'd do well to explain exactly what integral you actually want to evaluate. It seems unlikely that the strategy you're suggesting would be optimal. Is it just $\int_{-\infty}^{c} e^{-x^2 }\mbox{erf}(x+a)dx$? $\endgroup$ – Brian Borchers Jul 22 '16 at 0:48
  • 4
    $\begingroup$ If all you need is the definite integral $\int_{-\infty}^{\infty} \exp(-x^2)\,\text{erf}(x+a)\,dx$, it has the closed form $\sqrt{\pi}\, \text{erf}(a/\sqrt{2})$. Mathematica does not know this, but it is 4.3.13 in Geller and Ng's paper at nvlpubs.nist.gov/nistpubs/jres/73B/jresv73Bn1p1_A1b.pdf. You can prove it by expanding erf as an integral of $y$, and then changing variables with a 45-degree rotation of the $xy$-plane. $\endgroup$ – Matt F. Jul 22 '16 at 4:23
  • $\begingroup$ @BrianBorchers: Things like that, but possibly with different values of $c$ in the kernel of $e^{-x^c-y^c}$. The case $c=2$, a.k.a. Gaussian distribution, is the one that led to the question and is the most useful for this application. However, the case $c=1.5$, a.k.a. Holtsmark distribution, is also very useful in that context. $\endgroup$ – Michael Jul 22 '16 at 5:46

It would be worth looking at Schraudolph's sneaky trick if performance is your aim.

| cite | improve this answer | |

In general, approximation with least maximal error is possible with the Remez Algorithm.

In your case I would suggest giving approximation with Chebyshev rational functions a try.

| cite | improve this answer | |

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.