¿Existe una función de error complementaria escalada en python disponible?

En matlab hay una función especial que no está disponible en ninguna de las colecciones para Python que conozco (numpy, scipy, mpmath, …).

Probablemente hay otros lugares donde se pueden encontrar funciones como esta?

UPD Para todos los que piensan que la pregunta es trivial, primero trate de calcular esta función para el argumento ~ 30.

UPD2 La precisión arbitraria es una buena solución, pero si es posible preferiría evitarla. Necesito una precisión de máquina “estándar” (ni más ni menos) y una velocidad máxima posible.

UPD3 Resulta que mpmath da un resultado sorprendentemente inexacto. Incluso donde funciona la math estándar de python, los resultados de mpmath son peores. Lo que lo hace absolutamente inútil.

    UPD4 El código para comparar diferentes formas de calcular erfcx.

     import numpy as np def int_erfcx(x): "Integral which gives erfcx" from scipy import integrate def f(xi): return np.exp(-x*xi)*np.exp(-0.5*xi*xi) return 0.79788456080286535595*integrate.quad(f, 0.0,min(2.0,50.0/(1.0+x))+100.0,limit=500)[0] def my_erfcx(x): """MM Shepherd and JG Laframboise, MATHEMATICS OF COMPUTATION 36, 249 (1981) Note that it is reasonable to compute it in long double (or whatever python has) """ ch_coef=[np.float128(0.1177578934567401754080e+01), np.float128( -0.4590054580646477331e-02), np.float128( -0.84249133366517915584e-01), np.float128( 0.59209939998191890498e-01), np.float128( -0.26658668435305752277e-01), np.float128( 0.9074997670705265094e-02), np.float128( -0.2413163540417608191e-02), np.float128( 0.490775836525808632e-03), np.float128( -0.69169733025012064e-04), np.float128( 0.4139027986073010e-05), np.float128( 0.774038306619849e-06), np.float128( -0.218864010492344e-06), np.float128( 0.10764999465671e-07), np.float128( 0.4521959811218e-08), np.float128( -0.775440020883e-09), np.float128( -0.63180883409e-10), np.float128( 0.28687950109e-10), np.float128( 0.194558685e-12), np.float128( -0.965469675e-12), np.float128( 0.32525481e-13), np.float128( 0.33478119e-13), np.float128( -0.1864563e-14), np.float128( -0.1250795e-14), np.float128( 0.74182e-16), np.float128( 0.50681e-16), np.float128( -0.2237e-17), np.float128( -0.2187e-17), np.float128( 0.27e-19), np.float128( 0.97e-19), np.float128( 0.3e-20), np.float128( -0.4e-20)] K=np.float128(3.75) y = (xK) / (x+K) y2 = np.float128(2.0)*y (d, dd) = (ch_coef[-1], np.float128(0.0)) for cj in ch_coef[-2:0:-1]: (d, dd) = (y2 * d - dd + cj, d) d = y * d - dd + ch_coef[0] return d/(np.float128(1)+np.float128(2)*x) def math_erfcx(x): import scipy.special as spec return spec.erfc(x) * np.exp(x*x) def mpmath_erfcx(x): import mpmath return mpmath.exp(x**2) * mpmath.erfc(x) if __name__ == "__main__": x=np.linspace(1.0,26.0,200) X=np.linspace(1.0,100.0,200) intY = np.array([int_erfcx(xx*np.sqrt(2)) for xx in X]) myY = np.array([my_erfcx(xx) for xx in X]) myy = np.array([my_erfcx(xx) for xx in x]) mathy = np.array([math_erfcx(xx) for xx in x]) mpmathy = np.array([mpmath_erfcx(xx) for xx in x]) mpmathY = np.array([mpmath_erfcx(xx) for xx in X]) print ("Integral vs exact: %g"%max(np.abs(intY-myY)/myY)) print ("math vs exact: %g"%max(np.abs(mathy-myy)/myy)) print ("mpmath vs math: %g"%max(np.abs(mpmathy-mathy)/mathy)) print ("mpmath vs integral:%g"%max(np.abs(mpmathY-intY)/intY)) exit() 

    Para mi, da

     Integral vs exact: 6.81236e-16 math vs exact: 7.1137e-16 mpmath vs math: 4.90899e-14 mpmath vs integral:8.85422e-13 

    Obviamente, las math dan la mejor precisión posible donde funciona, mientras que mpmath da un par de errores de magnitud mayor en el caso de las math e incluso más para argumentos más grandes.

    Aquí hay una implementación simple y rápida que da una precisión de 12-13 dígitos globalmente:

     from scipy.special import exp, erfc def erfcx(x): if x < 25: return erfc(x) * exp(x*x) else: y = 1. / x z = y * y s = y*(1.+z*(-0.5+z*(0.75+z*(-1.875+z*(6.5625-29.53125*z))))) return s * 0.564189583547756287 

    La biblioteca gmpy2 proporciona acceso a la biblioteca de precisión múltiple MPFR. Para una precisión normal, es casi 5 veces más rápido que mpmath.

     $ py27 -m timeit -s "import mpmath" -s "def erfcx(x):return mpmath.exp(x**2) * mpmath.erfc(x)" "erfcx(30)" 10000 loops, best of 3: 47.3 usec per loop $ py27 -m timeit -s "import gmpy2" -s "def erfcx(x):return gmpy2.exp(x**2) * gmpy2.erfc(x)" "erfcx(30)" 100000 loops, best of 3: 10.8 usec per loop 

    Ambas bibliotecas devuelven el mismo resultado para 30.

     >>> import mpmath >>> import gmpy2 >>> mpmath.exp(30**2) * mpmath.erfc(30) mpf('0.018795888861416751') >>> gmpy2.exp(30**2) * gmpy2.erfc(30) mpfr('0.018795888861416751') >>> 

    Descargo de responsabilidad: mantengo gmpy2. Estoy trabajando activamente hacia una nueva versión, pero no debería haber problemas con la versión actual para este cálculo.

    Edición: Tenía curiosidad acerca de la sobrecarga de realizar llamadas a dos funciones en lugar de solo una, así que implementé un gmpy2.erfcx () completamente en C, pero aún utilizaba MPFR para realizar los cálculos. La mejora fue menor de lo que esperaba. Si crees que erfcx () sería útil, puedo agregarlo a la próxima versión.

     $ py27 -m timeit -s "import gmpy2" "gmpy2.erfcx(30)" 100000 loops, best of 3: 9.45 usec per loop 

    Una implementación C ++ altamente optimizada de erfcx (tanto para argumentos reales como complejos) se fusionó recientemente en SciPy y debería estar en la versión 0.12 de SciPy.

    No sé si alguna de las fonts estándar incluye esa función, pero puede implementarla de manera directa, al menos si usa mpmath y no está demasiado preocupado por el rendimiento:

     import math import mpmath def erfcx(x): return math.exp(x**2) * math.erfc(x) def erfcx_mp(x): return mpmath.exp(x**2) * mpmath.erfc(x) where = mpmath.linspace(1, 50, 10) + mpmath.linspace(100, 1000, 5) for x in where: try: std = erfcx(x) except OverflowError: std = None new = erfcx_mp(x) approx = (1/(x*mpmath.pi**0.5)) print x, std, new, (new-approx)/approx 

    produce

     1.0 0.427583576156 0.427583576155807 -0.242127843858688 6.44444444444444 0.0865286153111 0.0865286153111425 -0.0116285899486798 11.8888888888889 0.0472890800456 0.0472890800455829 -0.00350053472385845 17.3333333333333 0.032495498521 0.0324954985209682 -0.00165596082986796 22.7777777777778 0.024745497 0.0247454970000106 -0.000960939188986022 28.2222222222222 None 0.0199784436993529 -0.000626572735073611 33.6666666666667 None 0.0167507236463156 -0.000440550710337029 39.1111111111111 None 0.0144205913280408 -0.000326545959369654 44.5555555555556 None 0.0126594222570918 -0.00025167403795913 50.0 None 0.0112815362653238 -0.000199880119832415 100.0 None 0.00564161378298943 -4.99925018743586e-5 325.0 None 0.00173595973189465 -4.73366058776083e-6 550.0 None 0.00102579754728657 -1.6528843659911e-6 775.0 None 0.000727985953393782 -8.32464102161289e-7 1000.0 None 0.000564189301453388 -4.9999925011689e-7 

    Y se comporta como debería incluso cuando las rutinas matemáticas. * Se desbordan. El soporte de intervalo de mpmath no está a la altura de la tarea (sin un poco de piratería, soy demasiado vago para hacerlo), pero para esto estoy bastante seguro de que mpfs será suficiente, ya que erfcx es simplemente el producto de dos cosas que mpmath puede calcular muy bien.

    Esta función ahora está disponible en la versión más reciente de scipy, consulte http://docs.scipy.org/doc/scipy/reference/special.html .