f[x_] := Exp[-x];
c = 0.9;
g[x_] := c*x^(c - 1)*Exp[-x^c];
SetPrecision[Integrate[f[x]*Log[f[x]/g[x]], {x, 0.001, \[Infinity]}],20]

from scipy.integrate import quad
import numpy as np

def f(x):
    return sum([ps[idx]*lambdas[idx]*np.exp(- lambdas[idx] * x) for idx in range(len(ps))])

def g(x):
    return scipy.stats.weibull_min.pdf(x, c=c)

c = 0.9
ps = [1]
lambdas = [1]
eps = 0.001  # weibull_min is only defined for x > 0
print(quad(lambda x: f(x) * np.log(f(x) / g(x)), eps, np.inf)) # Output

一个非常有趣的问题。

首先要注意的是

from numpy import exp

def f(x):
    return exp(-x)

def g(x):
    c = 0.9
    return c * x**(c - 1) * exp(-x ** c)

def integrand(x):
    return f(x) * log(f(x) / g(x))

在0处具有 可积分 的奇点，并且[0，infty]上的积分可以通过 分析 求值。经过一些操作后，您会发现

import numpy
import scipy.special

c = 0.9

# euler_mascheroni constant
gamma = 0.57721566490153286060
val = scipy.special.gamma(c + 1) - 1 - numpy.log(c) + (c - 1) * gamma

print(val)



0.0094047810750603

wolfram-
alpha
正确地将其值赋予许多数字。为了用数值方法重现这一点，一个好的尝试总是总是tanh-
sinh正交
（例如，来自我的一个项目Quadpy）。以某个较大的值切断该域，无论如何该函数几乎都为0，然后：

from numpy import exp, log
import quadpy


def f(x):
    return exp(-x)


def g(x):
    c = 0.9
    return c * x**(c - 1) * exp(-x ** c)


def integrand(x):
    return f(x) * log(f(x) / g(x))


val, err = quadpy.tanh_sinh(integrand, 0.0, 100.0, 1.0e-8)
print(val)



0.009404781075063085

---

现在，对于其他一些事情，也许令人惊讶的是，它做得 不太 好。

当看到这种类型的整数时exp(-x) * f(x)，首先想到的是高斯-
拉格瑞正交
。例如，使用Quadpy（我的一个项目）：

import numpy
import quadpy

c = 0.9


def f(x):
    return numpy.exp(-x)


def g(x):
    return c * x ** (c - 1) * numpy.exp(-x ** c)


scheme = quadpy.e1r.gauss_laguerre(100)
val = scheme.integrate(lambda x: numpy.log(f(x) / g(x)))

print(val[0])

这给

0.010039543105755215

尽管我们使用了100个积分点，但对于实际值而言却是一个令人惊讶的差值。这是由于以下事实：多项式，尤其是项log(x)和不能很好地逼近被积数x ** c：

import numpy
from numpy import exp, log, ones
from scipy.special import gamma
import quadpy


c = 0.9


def integrand(x):
    return exp(-x) * (-x - log(c) - (c - 1) * log(x) - (-x ** c))


scheme = quadpy.e1r.gauss_laguerre(200)
val = scheme.integrate(lambda x: -x - log(c) - (c - 1) * log(x) - (-x ** c))[0]

vals = numpy.array([
    - scheme.integrate(lambda x: x)[0],
    -log(c) * scheme.integrate(lambda x: ones(x.shape))[0],
    -(c - 1) * scheme.integrate(lambda x: log(x))[0],
    scheme.integrate(lambda x: x ** c)[0]
])
euler_mascheroni = 0.57721566490153286060
exact = numpy.array([
    -1.0,
    -log(c),
    euler_mascheroni * (c-1),
    gamma(c + 1)
])
print("approximation, exact, diff:")
print(numpy.column_stack([vals, exact, abs(vals - exact)]))
print()
print("sum:")
print(sum(vals))



approximation, exact, diff:
[[-1.00000000e+00 -1.00000000e+00  8.88178420e-16]
 [ 1.05360516e-01  1.05360516e-01  6.93889390e-17]
 [-5.70908293e-02 -5.77215665e-02  6.30737142e-04]
 [ 9.61769857e-01  9.61765832e-01  4.02488825e-06]]

sum:
0.010039543105755278