Python中是否有易于获取的erf()实现？

Question

我能够自己实现误差函数erf，但我不希望这样做。有没有一个Python包含有这个函数的实现且不需要外部依赖？我找到了这个，但它似乎是某个较大包的一部分（而且不清楚是哪个！）。

- rog

Answer 1

自从版本2.7起，标准的math模块包含erf函数。这应该是最简单的方法。 http://docs.python.org/2/library/math.html#math.erf

- golobor

---

有没有一个Python模块提供erf⁻¹(x)函数？ - Lori

@Lori - 是的，math.erfc。 - Matthew

根据文档，我认为@Matthew的说法是不正确的。这会计算 1.0 - erf(x)。而erf()的逆函数在SciPy中可以找到：https://dev59.com/Q10Z5IYBdhLWcg3wnBVS - tmn

Answer 2

我建议使用SciPy进行Python中的数值函数，但如果您想要没有依赖项的内容，这里有一个函数，其误差小于所有输入的1.5 * 10^-7。

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)

该算法来自于数学函数手册，公式7.1.26。

- John D. Cook

---

你说得对。我编辑了我的答案，以不同的方式找到x的符号来解决这个问题。现在没问题了。 - John D. Cook

这段话的意思是：“'手册'是由美国联邦政府[雇员]编写的，没有版权保护。”下面提供了一本更直接的链接：http://www.math.sfu.ca/~cbm/aands/frameindex.htm - mariotomo

我试图使用这个函数，但是将numpy数组传递给x并返回numpy数组作为答案没有成功。你能否重写该函数的数组版本？ - Bob

这个近似只适用于实数x吗？ - poisson

Answer 3

我建议您下载NumPy（以在Python中拥有高效的矩阵）和SciPy（Matlab工具箱替代品，使用NumPy）。 erf函数位于SciPy中。

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

你也可以使用在pylab中定义的erf函数，但这更适用于绘制使用numpy和scipy计算出来的结果。如果你想要一个包含这些软件的完整安装包，你可以直接使用Python Enthought distribution。

- Mapad

---

SciPy是Python的数值软件的宝库。但是，要开始使用它可能有点具有挑战性。请从http://www.scipy.org/开始查看。 - John D. Cook

我必须说，我完全没能成功安装它。我要求一个没有外部依赖的软件包是有原因的。Numpy并不是唯一的一个。UMFPack也是如此。写自己的erf()会更容易！ - rog

尝试使用我提到的Python Enthought，他们已经捆绑了你所需要的一切。 - Mapad

363MB的受许可限制的下载超出了我为了30行函数愿意承受的范围...碰巧，似乎macports可以使用scipy。正在我说话的同时，它正在重新编译整个gcc。这里有些问题！ - rog

示例应为 "import scipy.special as erf"，即去掉第一个 "erf" 保留第二个。因为 erf 是一个函数而不是模块，所以它不应该是导入路径的一部分。"import scipy.special as foobar" 也可以工作。"as" 是一种方便的写法。 - John D. Cook

最好使用“from scipy.special import erf”来导入erf函数，这样会更清晰明了。 - Lanny

Answer 4

可以在 mpmath 模块中找到纯 Python 实现（http://code.google.com/p/mpmath/）

从文档字符串中：

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

对于较大的实数x，\mathrm{erf}(x)会非常迅速地接近1：

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

错误函数是一个奇函数：

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

:func:erf 实现任意精度计算并支持复数：

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

相关函数

另请参阅 :func:erfc，它对于大的x更加精确， 以及 :func:erfi，它给出了\exp(t^2)的反导数。

Fresnel积分 :func:fresnels和 :func:fresnelc 也与误差函数相关。

- Charles McCreary

---

这真的很有趣。大概这个多精度实现比使用本地浮点数要慢得多吧？ - rog

Answer 5

为了回答我的问题，我最终使用了以下代码，该代码改编自我在网上找到的Java版本：

# 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

- rog

---

不错，但从2.7开始，我们应该使用 from math import erf（为了可移植性、准确性、速度等）。 - smci

Answer 6

我有一个函数需要调用10^5次erf。在我的机器上...

scipy.special.erf的时间为6.1秒。

erf数学函数手册需要8.3秒。

erf Numerical Recipes 6.2需要9.5秒。

（三次运行平均值，代码取自以上帖子）。

- meteore

---

使用 ctypes 从 libm.so（标准的 C 数学库，在这里是 64 位 Linux）调用 erf 函数的执行时间降至 5.6 秒。 - meteore

我还需要成千上万次的erf调用。根据你的数字，我选择使用scipy.special.erf。自那以后你有找到更快的方法吗？我一直在考虑使用查找表.. - jtlz2

顺便说一句，你是如何使用scipy的erf函数的？通过以下设置：from scipy.special import erf; import numpy as np; data = np.random.randn(10e5)，我从result = erf(data)中获得非常快的运行时间。特别是在这种情况下，每个循环需要32毫秒。如果我天真地循环遍历numpy数组中的所有元素，唯一会导致运行时间> 1秒的方法就是这样。 - 8one6

Answer 7

对于那些追求更高性能的人，有一个提示：如果可能的话，进行向量化。

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)

提供结果

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

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

- 8one6

Answer 8

从Python的math.erf函数文档中可以看出，它在近似计算中使用了高达50个术语：

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.

- rllb

Answer 9

SciPy有一个erf函数的实现，参见scipy.special.erf。

- teekarna