scipy 中 exp1 的日志

问题描述

以下

f_ejer1 <- function(x,y) {
  sqrt(144 + x^2) + sqrt(25+(y - x)^2) + sqrt(4 + (7 - y)^2)
}

x1_ejer1 <- seq(-1,1,len=50)
x2_ejer1 <- seq(-1,len=50)

z_ejer1 <- outer(x1_ejer1,x2_ejer1,f_ejer1)

persp(x1_ejer1,z_ejer1,theta=-30,phi=15,ticktype="detailed",col="lightblue",shade=0.2)

image(x1_ejer1,z_ejer1)
contour(x1_ejer1,nlevels=15,col="black",add=T)

大 from scipy.special import gamma gamma(x) 溢出。这就是 scipy 提供 x 的原因，它相当于 gammaln 并允许我们在日志空间中工作并避免溢出。

scipy 的 exp1 函数有没有类似的东西？我想要一些返回与下面相同但没有下溢的大 np.log(gamma(x)) 的东西：

（我认为这类似于 import numpy as np from scipy.special import exp1 def exp1ln(x): return np.log(exp1(x)) 的原因是因为 gammaln 属于同一函数系列，请参见此处：Incomplete Gamma function in scipy。）

谢谢

解决方法

对于实数 x，您可以使用 log(Ei(x)) 的级数展开式来计算 x → ∞（参见 {{3 }}) 这对于大 x 非常准确。

根据一些快速实验，当 x >= 50（和这就是 scipy 开始失去完全精度的地方）。而且系列扩展真的很好，系数是阶乘数，因此我们可以使用精确的评估而不会造成灾难性的后果取消：

def _exp1ln_taylor(x,k):
    r = 0; s = -1
    for i in range(k,-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

对于想要只处理单个数字而不是数组的 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,-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))

python scipy scipy