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 this but this seems to be part of some much larger package (and it's not even clear which one!).
9 Answers
Since v.2.7. the standard math module contains erf function. This should be the easiest way.
-
-
2
-
5@Matthew I don't think that's correct according to documentation. That calculates 1.0 - erf(x). The inverse of erf() is in SciPy: stackoverflow.com/questions/31266249/…– tmnCommented Oct 26, 2019 at 17:45
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 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)
The algorithm comes from Handbook of Mathematical Functions, formula 7.1.26.
-
2You're right. I edited my answer to find the sign of x a different way to fix this problem. Now it's OK. Commented Jan 19, 2009 at 16:01
-
6from 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: math.sfu.ca/~cbm/aands/frameindex.htm Commented Nov 13, 2009 at 7:56
-
I am trying to use this function but passing numpy array to x and return numy array as answer with no success. Can you rewrite array version of the function?– BobCommented Aug 10, 2018 at 3:51
-
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.
>>>from scipy.special import 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.
-
2I 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()!– rogCommented Jan 19, 2009 at 13:35
-
1try Python Enthought as I mentioned, they've bundled everything you need.– MapadCommented Jan 19, 2009 at 13:51
-
1363MB 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!– rogCommented Jan 19, 2009 at 14:45
-
1The 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. Commented Jan 20, 2009 at 20:07
-
2It would be less confusing to "from scipy.special import erf" as that imports only erf function.– LannyCommented Nov 27, 2010 at 13:20
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.
-
that's really interesting. presumably this multi-precision implementation is a fair bit slower than using native floating point?– rogCommented Jan 22, 2009 at 17:28
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
-
2Nice, but as of 2.7 we should do
from math import erf(for portability, accuracy, speed, etc.)– smciCommented May 20, 2013 at 23:54
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).
-
1erf called with ctypes from libm.so (standard c math library, 64 bit linux here) goes down to 5.6s.– meteoreCommented Oct 26, 2010 at 14:44
-
I also require 1000s of erf calls. Based on your numbers I'm going with scipy.special.erf. Have you found anything faster since? I've wondered about using a lookup..– jtlz2Commented Aug 14, 2013 at 13:47
-
1FWIW, how are you using scipy's
erffunction? With the following setup:from scipy.special import erf; import numpy as np; data = np.random.randn(10e5), I get very fast runtimes from:result = erf(data). In particular, I get 32ms per loop in that case. The only way I get runtimes > 1s is if I naively loop over all the elements in anumpyarray.– 8one6Commented Mar 4, 2014 at 16:15
One note for those aiming for higher performance: vectorize, if 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)
gives results
# vectorized
10 loops, best of 3: 108 ms per loop
# loops
1 loops, best of 3: 2.34 s per loop
-
10e5 is really misleading, as it actually means 1000000 (10^6). Also, it has a type of float, but
np.random.randnaccepts only integers (at least in recent versions of NumPy; probably that was not the case 10 years ago). Commented Feb 12 at 17:05
From Python's math.erf function documentation, it uses up to 50 terms in the approximation:
Implementations of the error function erf(x) and the complementary error
function erfc(x).
Method: we use a series approximation for erf for small x, and a continued
fraction approximation for erfc(x) for larger x;
combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x),
this gives us erf(x) and erfc(x) for all x.
The series expansion used is:
erf(x) = x*exp(-x*x)/sqrt(pi) * [
2/1 + 4/3 x**2 + 8/15 x**4 + 16/105 x**6 + ...]
The coefficient of x**(2k-2) here is 4**k*factorial(k)/factorial(2*k).
This series converges well for smallish x, but slowly for larger x.
The continued fraction expansion used is:
erfc(x) = x*exp(-x*x)/sqrt(pi) * [1/(0.5 + x**2 -) 0.5/(2.5 + x**2 - )
3.0/(4.5 + x**2 - ) 7.5/(6.5 + x**2 - ) ...]
after the first term, the general term has the form:
k*(k-0.5)/(2*k+0.5 + x**2 - ...).
This expansion converges fast for larger x, but convergence becomes
infinitely slow as x approaches 0.0. The (somewhat naive) continued
fraction evaluation algorithm used below also risks overflow for large x;
but for large x, erfc(x) == 0.0 to within machine precision. (For
example, erfc(30.0) is approximately 2.56e-393).
Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and
continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) <
ERFC_CONTFRAC_CUTOFF. ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the
numbers of terms to use for the relevant expansions.
#define ERF_SERIES_CUTOFF 1.5
#define ERF_SERIES_TERMS 25
#define ERFC_CONTFRAC_CUTOFF 30.0
#define ERFC_CONTFRAC_TERMS 50
Error function, via power series.
Given a finite float x, return an approximation to erf(x).
Converges reasonably fast for small x.