Journal de l'expérience 1 dans scipy

Question

Le suivant

from scipy.special import gamma

gamma(x)

déborde pour de grands x. C'est pourquoi scipy fournit gammaln, qui est équivalent à np.log(gamma(x)) et nous permet de travailler dans l'espace logarithmique et d'éviter le débordement.

Y a-t-il quelque chose de similaire pour la fonction exp1 de scipy? J'aimerais avoir quelque chose qui retourne la même chose que ci-dessous mais sans avoir de sous-débordement pour de grands x:

import numpy as np
from scipy.special import exp1

def exp1ln(x):
    return np.log(exp1(x))

(Ma raison de penser que ceci serait similaire à gammaln est que exp1 est de la même famille de fonctions, voir ici: Fonction Gamma incomplète dans scipy .)

Merci

Answer 1

Pour un x réel, vous pouvez utiliser le développement en série du logarithme de Ei(x) pour x (voir ici) qui est très précis pour de grandes valeurs de x.

D'après quelques expérimentations rapides, 18 termes sont suffisants pour une précision totale en virgule flottante lorsque x >= 50 (c'est à peu près là que scipy commence à perdre sa précision totale). De plus, le développement en série est vraiment sympa, les coefficients étant les nombres factoriels, donc nous pouvons utiliser une évaluation précise sans annulation catastrophique :

def _exp1ln_taylor(x, k):
    r = 0; s = -1
    for i in range(k, 0, -1):
        r = (r + s) * i / x
        s = -s
    return r*s

def exp1ln(x):
    x = np.array(x)
    use_approx = x >= 50

    # Éviter les infinis/sous-débordements en dehors du domaine.
    ax = np.where(use_approx, x, 100)
    sx = np.where(use_approx, 1, x)
    approx = -x + -np.log(x) + np.log1p(_exp1ln_taylor(ax, 18))
    sp = np.log(scipy.special.exp1(sx))
    return np.where(use_approx, approx, sp) * 1

Comparaison avec WolframAlpha (qui peut fournir des centaines de chiffres) :

x           exp1ln(x)           WolframAlpha
0.1         0.6004417824948862  0.6004417824948862
1          -1.5169319590020440 -1.5169319590020456
10         -12.390724371937408 -12.390724371937408
100        -104.61502435050535 -104.61502435050535
1000       -1006.9087537832978 -1006.9087537832978
1000000000 -1000000020.7232659 -1000000020.7232658

Answer 2

Pour ceux qui veulent une version de exp1ln qui ne traite que des nombres simples au lieu de tableaux, j'ai adapté la solution de orlp pour obtenir ceci :

import numpy as np
from scipy.special import exp1

def _exp1ln_taylor(x, k):
    r = 0
    s = -1
    for i in range(k, 0, -1):
        r = (r + s) * i / x
        s *= -1
    return r * s

def exp1ln(x):
    if x <= 50:
        return np.log(exp1(x))
    else:
        return -x - np.log(x) + np.log1p(_exp1ln_taylor(x, 18))