Existe-t-il une implémentation facilement disponible de erf() pour Python ?

Question

Je peux implémenter la fonction d'erreur, erf, moi-même, mais je préfère ne pas le faire. Existe-t-il un paquetage python sans dépendances externes qui contient une implémentation de cette fonction ? J'ai trouvé http://pylab.sourceforge.net/packages/included\_functions.html>this mais cela semble faire partie d'un paquetage beaucoup plus large (et on ne sait même pas lequel !).

Je suis désolé si cette question est naïve - je suis totalement novice en matière de Python.

Answer 1

Pour répondre à ma propre question, j'ai fini par utiliser le code suivant, adapté d'une version Java que j'ai trouvée ailleurs sur le 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

Answer 2

Une remarque pour ceux qui visent des performances plus élevées : vectorisez, si possible.

import numpy as np
from scipy.special import erf

def vectorized(n):
    x = np.random.randn(n)
    return erf(x)

def loopstyle(n):
    x = np.random.randn(n)
    return [erf(v) for v in x]

%timeit vectorized(10e5)
%timeit loopstyle(10e5)

donne des résultats

# vectorized
10 loops, best of 3: 108 ms per loop

# loops
1 loops, best of 3: 2.34 s per loop