您如何像 Mathematica 那样执行这种不正确的积分？

Question

使用这个 Mathematica 代码：

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]

Mathematica 毫无问题地计算了这个并给出了答案 0.010089328699390866240。我希望能够执行类似的积分，但我没有 Mathematica 的副本。只是天真地在 scipy 中实现它，例如，使用标准正交库可悲地失败了，因为 f(x) 和 g(x) 任意接近 0。这是 Python 中使用标准正交的示例由于所需的无限精度而失败。:

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

应该大于 0

如何在代码中像 Mathematica 那样执行这种不正确的积分？我不介意使用哪个免费 language/library。

Answer 1

一个很有趣的问题。

首先注意被积函数

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 gives its value correctly to many digits. To reproduce this with numerical methods, a good first try is always tanh-sinh quadrature (e.g., from 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)类型的积分，首先想到的应该是Gauss-Laguerre quadrature. For example with 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

Answer 2

在julia中，QuadGK包可以做这些积分。直接这样做你会遇到问题，正如你所注意到的：

f(x) = exp(-x)
g(x; c=0.9) = c*x^(c - 1)*exp(-x^c)
h(x) = f(x) * log(f(x)/g(x))
using QuadGK
a,b = 0.001, Inf
quadgk(h, a, b)  # errors

但是将 log(f/g) 扩展为 log(f) - (log(c) + (c-1)log(x) + x^c) 我们可以得到每个项进行积分：

c = 0.9
quadgk(x -> f(x) * -x, a,b)
quadgk(x -> -f(x)*log(c), a,b)
quadgk(x -> -f(x)*(c-1)*log(x), a,b)
quadgk(x -> f(x) * x^c, a,b)

将这些值相加得出答案。

您也可以通过过滤掉 NaN 值来获得答案，这可能效率低得多：

h1(x) = isnan(h(x)) ? 0.0 : h(x)
quadgk(h1, a,b) # (0.010089328699390816, 9.110982026738999e-11)

使用big(a)和big(b)可以获得更多的小数点。

您如何像 Mathematica 那样执行这种不正确的积分？

How can you perform this improper integral as Mathematica does?

python

r

math

julia

numerical-integration