Tricking numpy/Python in repräsentieren sehr große und sehr kleine Zahlen

Question

Ich brauche das Integral der folgenden Funktion in Bereichen zu berechnen, die so günstig wie -150 starten:

import numpy as np 
from scipy.special import ndtr 

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

Das Problem ist, dass dieser Teil der Funktion

np.exp(x ** 2) 

gegen unendlich geht - ich inf für Werte von x weniger als etwa -26 bekommen.

Und dieser Teil der Funktion

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

, die

from scipy.special import erf 

1 + erf(x) 

entspricht tendiert zu 0.

Also, eine sehr, sehr große Anzahl mal ein sehr, sehr klein Nummer sollte mir eine vernünftige Größe geben - aber stattdessen gibt python mir nan.

Was kann ich tun, um dieses Problem zu umgehen?

Answer 1

Es gibt bereits eine solche Funktion: erfcx. Ich denke, erfcx(-x) sollte Ihnen den Integrand geben, den Sie wollen (beachten Sie, dass 1+erf(x)=erfc(-x)).

Answer 2

denke ich, eine Kombination aus Lösung des @ askewchan und scipy.special.log_ndtr den Trick:

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 

Für x <= -20, log_ndtr(x)uses a Taylor series expansion of the error function to iteratively compute the log CDF directly, die numerisch viel mehr stabil ist, als nur log(ndtr(x)) nehmen.

---

aktualisieren

Wie Sie in den Kommentaren erwähnt, die exp kann auch Überlauf, wenn x ausreichend groß ist. Während könnte man dieses Problem umgehen mpmath.exp verwenden, eine einfachere und schnellere Methode ist zu einem np.longdouble schütten, die auf meinem Rechner, können Werte repräsentieren bis + 4932 bis 1.189731495357231765e:

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

Nicht sicher, wie nützlich wird diese Seien Sie, aber hier sind ein paar Gedanken, die für einen Kommentar zu lang sind.

Sie müssen das Integral von berechnen, das correctly identified wäre . Wenn Sie die Klammern öffnen, können Sie beide Teile der Summierung integrieren.

Scipy hat dieses imaginary error function implemented

Der zweite Teil ist härter:

Dies ist ein generalized hypergeometric function ist. Leider sieht es aus wie scipy does not have an implementation of it, aber this package behauptet, dass es tut.

Hier verwendete ich unbestimmte Integrale ohne Konstanten, die fromto Werte wissen, es ist klar, wie man bestimmte verwendet.