Oszukiwanie numpy/Python do reprezentującą bardzo dużych i bardzo małych ilościach

Question

muszę obliczyć całkę następującej funkcji w zakresach, które zaczyna się już od -150:

import numpy as np 
from scipy.special import ndtr 

def my_func(x): 
    return np.exp(x ** 2) * 2 * ndtr(x * np.sqrt(2)) 

Problemem jest to, że ta część funkcji

np.exp(x ** 2) 

dąży do nieskończoności - mam inf dla wartości x mniej niż około -26.

I ta część funkcji

2 * ndtr(x * np.sqrt(2)) 

co jest równoważne

from scipy.special import erf 

1 + erf(x) 

dąży do 0.

tak, to bardzo, bardzo dużych ilość razy bardzo, bardzo małe numer powinien podać mi odpowiednią liczbę - ale zamiast tego python daje mi nan.

Co mogę zrobić, aby obejść ten problem?

Answer 1

Istnieje już taka funkcja: erfcx. Myślę, że erfcx(-x) powinien dać ci integrand, którego potrzebujesz (zauważ, że 1+erf(x)=erfc(-x)).

Answer 2

myślę kombinacja rozwiązania @ askewchan i scipy.special.log_ndtr rade:

from scipy.special import log_ndtr 

_log2 = np.log(2) 
_sqrt2 = np.sqrt(2) 

def my_func(x): 
    return np.exp(x ** 2) * 2 * ndtr(x * np.sqrt(2)) 

def my_func2(x): 
    return np.exp(x * x + _log2 + log_ndtr(x * _sqrt2)) 

print(my_func(-150)) 
# nan 

print(my_func2(-150) 
# 0.0037611803122451198 

Dla x <= -20, log_ndtr(x)uses a Taylor series expansion of the error function to iteratively compute the log CDF directly, która jest o wiele bardziej stabilne numerycznie niż po prostu biorąc log(ndtr(x)).

---

Aktualizacja

Jak wspomniano w uwag exp może również przepełnienia jeśli x jest wystarczająco duża. Chociaż można to obejść stosując mpmath.exp, prostsza i szybsza metoda jest rzucić się do np.longdouble które na moim komputerze, może reprezentować wartości do 1.189731495357231765e + 4932:

import mpmath 

def my_func3(x): 
    return mpmath.exp(x * x + _log2 + log_ndtr(x * _sqrt2)) 

def my_func4(x): 
    return np.exp(np.float128(x * x + _log2 + log_ndtr(x * _sqrt2))) 

print(my_func2(50)) 
# inf 

print(my_func3(50)) 
# mpf('1.0895188633566085e+1086') 

print(my_func4(50)) 
# 1.0895188633566084842e+1086 

%timeit my_func3(50) 
# The slowest run took 8.01 times longer than the fastest. This could mean that 
# an intermediate result is being cached 100000 loops, best of 3: 15.5 µs per 
# loop 

%timeit my_func4(50) 
# The slowest run took 11.11 times longer than the fastest. This could mean 
# that an intermediate result is being cached 100000 loops, best of 3: 2.9 µs 
# per loop 

Answer 3

Nie wiem, jak to będzie pomocne być, ale oto kilka myśli, które są zbyt długie na komentarz.

Musisz obliczyć całkę z , która będzie correctly identified będzie . Otwierając nawiasy, możesz zintegrować obie części sumy.

scipy ma ten imaginary error function implemented

Druga część jest trudniejsze:

To generalized hypergeometric function. Niestety wygląda to na scipy does not have an implementation of it, ale twierdzi, że tak jest.

Tutaj użyłem nieoznaczonych całek bez stałych, znając wartości, które są jasne, jak używać określonych.