Có hàm lỗi bổ sung được chia tỷ lệ trong trăn có sẵn không?

Question

Trong khung chiếu có special function không có sẵn trong bất kỳ bộ sưu tập nào cho Python mà tôi biết (numpy, scipy, mpmath, ...).

Có thể có những nơi khác có thể tìm thấy các chức năng như thế này?

UPD Đối với tất cả những ai nghĩ rằng câu hỏi là tầm thường, hãy thử tính toán hàm này cho đối số ~ 30 trước tiên.

UPD2 Độ chính xác tùy ý là giải pháp tốt, nhưng nếu có thể, tôi muốn tránh điều đó. Tôi cần độ chính xác của máy "chuẩn" (không ít hơn) và tốc độ tối đa có thể.

UPD3 Hóa ra, mpmath cho kết quả không chính xác đáng ngạc nhiên. Ngay cả khi tiêu chuẩn python math hoạt động, kết quả là mpmath tồi tệ hơn. Mà làm cho nó hoàn toàn vô giá trị.

UPD4 Mã để so sánh các cách khác nhau để tính toán 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() 

Đối với tôi, nó mang lại cho

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

Rõ ràng, math cho độ chính xác tốt nhất có thể mà nó hoạt động trong khi mpmath cho vài đơn đặt hàng lỗi của cường độ lớn hơn nơi math công trình, thậm chí nhiều hơn cho các đối số lớn hơn.

Answer 1

Đây là một thực hiện đơn giản và nhanh chóng đưa ra 12-13 chữ số chính xác toàn cầu:

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

Tôi không biết rằng bất kỳ trong những nguồn tiêu chuẩn bao gồm chức năng đó, nhưng bạn có thể thực hiện nó trong thời trang đơn giản, ít nhất là nếu bạn sử dụng mpmath và không quá lo lắng về hiệu suất:

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 

sản xuất

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 

Và nó cư xử như nó nên thậm chí khi việc tính toán. * thói quen tràn. Sự hỗ trợ khoảng thời gian của mpmath không hoàn toàn phụ thuộc vào nhiệm vụ (không có một số hackery tôi quá lười để làm), nhưng với điều này, tôi chắc chắn rằng mpfs sẽ đủ, vì erfcx đơn giản là sản phẩm của hai thứ mà mpmath có thể tính toán độc đáo.

Answer 3

Thư viện gmpy2 cung cấp quyền truy cập vào thư viện nhiều độ chính xác MPFR. Đối với độ chính xác bình thường, nó gần như là 5x nhanh hơn 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 

Cả hai thư viện trả về kết quả tương tự cho 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') 
>>> 

Disclaimer: Tôi duy trì gmpy2. Tôi đang tích cực làm việc hướng tới một bản phát hành mới nhưng không có bất kỳ vấn đề nào với bản phát hành hiện tại cho phép tính này.

Chỉnh sửa: Tôi đã tò mò về chi phí thực hiện hai cuộc gọi hàm thay vì chỉ một vì vậy tôi đã triển khai gmpy2.erfcx() hoàn toàn bằng C nhưng vẫn sử dụng MPFR để thực hiện các phép tính. Sự cải thiện nhỏ hơn tôi mong đợi. Nếu bạn nghĩ rằng erfcx() sẽ hữu ích, tôi có thể thêm nó vào bản phát hành tiếp theo.

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

Answer 4

Một C tối ưu hóa cao ++ thực hiện erfcx (đối với lập luận cả trong thực tế và phức tạp) mới đây đã merged into SciPy và phải ở trong scipy phiên bản 0.12.

Answer 5

Chức năng này hiện khả dụng trong Phiên bản mới nhất của scipy, xem http://docs.scipy.org/doc/scipy/reference/special.html.