¿Hay una función de error complementario escalada en python disponible?

Question

En matlab hay un special function que no está disponible en ninguna de las colecciones para el Python que conozco (numpy, scipy, mpmath, ...).

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

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

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 la velocidad máxima posible.

UPD3 Resulta que mpmath da un resultado sorprendentemente inexacto. Incluso cuando python estándar math funciona, 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): 
    """M. M. Shepherd and J. G. 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 = (x-K)/(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 mí, se 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, math da mayor precisión posible, donde funciona mientras mpmath da error par de órdenes de magnitud más grande donde math obras y aún más para los argumentos de mayor tamaño.

Answer 1

Aquí es una aplicación sencilla y rápida dando 12-13 exactitud dígitos a nivel mundial:

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 

Answer 2

No sé que cualquiera de las fuentes estándar incluyen esa función, pero puedo implementarlo de manera sencilla, al menos si se utiliza 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 debe, incluso cuando los cálculos. * rutinas de desbordamiento. El soporte de intervalos de mpmath no está a la altura de la tarea (sin algunos hackers, soy demasiado flojo 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.

Answer 3

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

responsabilidad: Mantengo gmpy2. Estoy trabajando activamente en una nueva versión, pero no debería haber ningún problema con la versión actual para este cálculo.

Editar: Tenía curiosidad por la sobrecarga de hacer dos llamadas a funciones en lugar de una sola, así que implementé un gmpy2.erfcx() completamente en C pero aún usando 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 

Answer 4

A C++ aplicación altamente optimizado de erfcx (para argumentos tanto reales y complejos) fue recientemente merged into SciPy y debe estar en SciPy versión 0.12.

Answer 5

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.