scipy中exp1的日志

Question

以下

from scipy.special import gamma

gamma(x)

大型x的溢出。这就是为什么scipy提供了gammaln，它相当于np.log(gamma(x))，允许我们在日志空间中工作并避免溢出

scipy的^{}函数有类似的功能吗？我想要返回与下面相同的内容，但是对于大的x没有下溢：

import numpy as np
from scipy.special import exp1

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

（我认为这类似于gammaln的原因是exp1属于同一个函数族，请参见此处：Incomplete Gamma function in scipy。）

谢谢

Answer 1

对于任何想要只处理单个数字而不是数组的exp1ln版本的人，我修改了orlp中的解决方案以获得以下结果：

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))

Answer 2

对于realx，您可以使用log的系列扩展（Ei（x）来表示x→ ∞ （参见here），这对于大型x非常准确

通过一些快速实验，当x>；=50（及 这就是scipy开始失去全部精度的地方）。而且系列的扩展也很不错， 系数是阶乘数，因此我们可以使用精确的计算，而不会产生灾难性的后果 取消：

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

    # Avoid infinity/underflows from outside domain.
    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

与WolframAlpha（可提供数百位数字）的比较：

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