python中是否有比例缩放的互补误差函数? [英] Is there a scaled complementary error function in python available?
问题描述
在matlab中有一个特殊功能，该功能在我知道的Python的任何集合(numpy，scipy，mpmath等).
In matlab there is a special function which is not available in any of the collections for the Python I know (numpy, scipy, mpmath, ...).
可能在其他地方可以找到类似这样的功能?
Probably there are other places where functions like this one may be found?
UPD 对于所有认为该问题无关紧要的人，请首先尝试为参数〜30计算此函数.
UPD For all who think that the question is trivial, please try to compute this function for argument ~30 first.
UPD2 任意精度是一种不错的解决方法，但如果可能的话，我宁愿避免使用它.我需要标准"的机器精度(不多不少)和可能的最大速度.
UPD2 Arbitrary precision is a nice workaround, but if possible I would prefer to avoid it. I need a "standard" machine precision (no more no less) and maximum speed possible.
UPD3 事实证明，mpmath给出了令人惊讶的不准确结果.即使在标准python math可以工作的地方，mpmath的结果也更糟.这使其绝对一文不值.
UPD3 It turns out, mpmath gives surprisingly inaccurate result. Even where standard python math works, mpmath results are worse. Which makes it absolutely worthless.
UPD4 该代码用于比较计算erfcx的不同方法.
UPD4 The code to compare different ways to compute 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()
对我来说，它给出了
Integral vs exact: 6.81236e-16
math vs exact: 7.1137e-16
mpmath vs math: 4.90899e-14
mpmath vs integral:8.85422e-13
很明显，math在工作时会提供最佳的精度，而mpmath在math工作时给出的误差对数量级要大一些，而对于更大的参数则要大得多.
Obviously, math gives best possible precision where it works while mpmath gives error couple orders of magnitude larger where math works and even more for larger arguments.
推荐答案
这是一个简单快捷的实现，全局精度为12到13位:
Here is a simple and fast implementation giving 12-13 digit accuracy globally:
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
这篇关于python中是否有比例缩放的互补误差函数?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
