Tuesday, April 5, 2011

Is there an easily available implementation of erf() for Python?

I can implement the error function, erf, myself, but I'd prefer not to. Is there a python package with no external dependencies that contains an implementation of this function? I have found http://pylab.sourceforge.net/packages/included_functions.html>this but this seems to be part of some much larger package (and it's not even clear which one!).

I'm sorry if this is a naive question - I'm totally new to Python.

From stackoverflow
  • I would recommend you download numpy (to have efficiant matrix in python) and scipy (a Matlab toolbox substitute, which uses numpy). The erf function lies in scipy.

    >>>import scipy.special.erf as erf
    >>>help(erf)
    

    You can also use the erf function defined in pylab, but this is more intended at plotting the results of the things you compute with numpy and scipy. If you want an all-in-one installation of these software you can use directly the Python Enthought distribution.

    John D. Cook : SciPy is the motherload of numerical software for Python. But getting started using it can be a little challenging. Start by looking at http://www.scipy.org/
    rog : I have to say that I've totally failed to install it. There was a reason that I asked for a package with no external dependencies. Numpy is not the only one. UMFPack is another. It'll be easier to write my own erf()!
    Mapad : try Python Enthought as I mentioned, they've bundled everything you need.
    rog : 363MB of license-encumbered download is further than I'm willing to go for a 30 line function... as it happens it seems that macports do scipy. It's recompiling the entirety of gcc as I speak. There's something wrong here!
    John D. Cook : The example should read "import scipy.special as erf", i.e. get rid of the first "erf" and keep the second. Since erf is a function and not a module, it shouldn't be part of the import path. "import scipy.special as foobar" would work too. The "as" is a convenience.
  • I recommend SciPy for numerical functions in Python, but if you want something with no dependencies, here is a function with an error error is less than 1.5 * 10-7 for all inputs.

    def erf(x):
        # save the sign of x
        sign = 1
        if x < 0: 
            sign = -1
        x = abs(x)
    
        # constants
        a1 =  0.254829592
        a2 = -0.284496736
        a3 =  1.421413741
        a4 = -1.453152027
        a5 =  1.061405429
        p  =  0.3275911
    
        # A&S formula 7.1.26
        t = 1.0/(1.0 + p*x)
        y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x)
        return sign*y # erf(-x) = -erf(x)
    

    The algorithm comes from Handbook of Mathematical Functions, formula 7.1.26.

    rog : This code gives a division-by-zero error for erf(0.0).
    John D. Cook : You're right. I edited my answer to find the sign of x a different way to fix this problem. Now it's OK.
    mariotomo : from wikipedia: "the 'handbook' is the work of US federal government [employees], not protected by copyright". I am putting here a more direct link to the book: http://www.math.sfu.ca/~cbm/aands/frameindex.htm
  • To answer my own question, I have ended up using the following code, adapted from a Java version I found elsewhere on the web:

    # from: http://www.cs.princeton.edu/introcs/21function/ErrorFunction.java.html
    # Implements the Gauss error function.
    #   erf(z) = 2 / sqrt(pi) * integral(exp(-t*t), t = 0..z)
    #
    # fractional error in math formula less than 1.2 * 10 ^ -7.
    # although subject to catastrophic cancellation when z in very close to 0
    # from Chebyshev fitting formula for erf(z) from Numerical Recipes, 6.2
    def erf(z):
     t = 1.0 / (1.0 + 0.5 * abs(z))
         # use Horner's method
            ans = 1 - t * math.exp( -z*z -  1.26551223 +
                 t * ( 1.00002368 +
                 t * ( 0.37409196 + 
                 t * ( 0.09678418 + 
                 t * (-0.18628806 + 
                 t * ( 0.27886807 + 
                 t * (-1.13520398 + 
                 t * ( 1.48851587 + 
                 t * (-0.82215223 + 
                 t * ( 0.17087277))))))))))
            if z >= 0.0:
             return ans
            else:
             return -ans
    
  • A pure python implementation can be found in the mpmath module (http://code.google.com/p/mpmath/)

    From the doc string:

    >>> from mpmath import *
    >>> mp.dps = 15
    >>> print erf(0)
    0.0
    >>> print erf(1)
    0.842700792949715
    >>> print erf(-1)
    -0.842700792949715
    >>> print erf(inf)
    1.0
    >>> print erf(-inf)
    -1.0
    

    For large real x, \mathrm{erf}(x) approaches 1 very rapidly::

    >>> print erf(3)
    0.999977909503001
    >>> print erf(5)
    0.999999999998463
    

    The error function is an odd function::

    >>> nprint(chop(taylor(erf, 0, 5)))
    [0.0, 1.12838, 0.0, -0.376126, 0.0, 0.112838]
    

    :func:erf implements arbitrary-precision evaluation and supports complex numbers::

    >>> mp.dps = 50
    >>> print erf(0.5)
    0.52049987781304653768274665389196452873645157575796
    >>> mp.dps = 25
    >>> print erf(1+j)
    (1.316151281697947644880271 + 0.1904534692378346862841089j)
    

    Related functions

    See also :func:erfc, which is more accurate for large x, and :func:erfi which gives the antiderivative of \exp(t^2).

    The Fresnel integrals :func:fresnels and :func:fresnelc are also related to the error function.

    rog : that's really interesting. presumably this multi-precision implementation is a fair bit slower than using native floating point?
  • I have a function which does 10^5 erf calls. On my machine...

    scipy.special.erf makes it time at 6.1s

    erf Handbook of Mathematical Functions takes 8.3s

    erf Numerical Recipes 6.2 takes 9.5s

    (three-run averages, code taken from above posters).

    meteore : erf called with ctypes from libm.so (standard c math library, 64 bit linux here) goes down to 5.6s.

0 comments:

Post a Comment