¿Hay una implementación fácilmente disponible de erf() para Python?

Question

Puedo implementar la función de error, erf, yo mismo, pero preferiría no hacerlo. ¿Hay un paquete de Python sin dependencias externas que contenga una implementación de esta función? He encontrado this pero esto parece ser parte de un paquete mucho más grande (¡y ni siquiera está claro cuál!).

Answer 1

Desde v.2.7. el módulo estándar math contiene erf función. Esta debería ser la forma más fácil.

http://docs.python.org/2/library/math.html#math.erf

Answer 2

Yo recomendaría que descargues numpy (para tener una matriz efficiant en python) y scipy (un sustituto de Matlab toolbox, que usa numpy). La función erf se encuentra en scipy.

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

También puede utilizar la función erf se define en pylab, pero esto es más destinados a trazar los resultados de las cosas que calcular con numpy y scipy. Si desea una instalación todo en uno de estos software, puede usar directamente el Python Enthought distribution.

Answer 3

Recomiendo SciPy para funciones numéricas en Python, pero si quiere algo sin dependencias, aquí hay una función con un error de menos de 1.5 * 10 ^-7 para todas las entradas.

def erf(x): 
    # save the sign of x 
    sign = 1 if x >= 0 else -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) 

El algoritmo procede de Handbook of Mathematical Functions, fórmula 7.1.26.

Answer 4

para responder a mi propia pregunta, he terminado usando el siguiente código, adaptado de una versión de Java que encontré en otro lugar en la 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 5

Una aplicación Python puro se puede encontrar en el módulo mpmath (http://code.google.com/p/mpmath/)

Desde la cadena de documentación:

>>> 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 

Para grandes verdadera x, \mathrm{erf}(x) approache s 1 muy rápidamente ::

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

la función de error es una función impar ::

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

: func: erf implementa la evaluación de precisión arbitraria y soporta números complejos ::

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

Funciones relacionadas

Consulte también: func: erfc, que es más preciso para grandes x, y : func: erfi que da la primitiva de \exp(t^2).

Las integrales de Fresnel: func: fresnels y: func: fresnelc también están relacionadas con la función de error.

Answer 6

Tengo una función que hace 10^5 llamadas erf. En mi máquina ...

scipy.special.erf hace la hora en 6.1s

Handbook of Mathematical Functions erf toma 8.3s

ERF Numerical Recipes 6.2 toma (promedios de tres carreras, código tomadas desde arriba posters) 9.5s

.

Answer 7

Una nota para aquellos que quieran para un mayor rendimiento: vectorizar, si es posible.

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) 

da resultados

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

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