您如何像Mathematica一样执行这种不适当的积分? [英] How can you perform this improper integral as Mathematica does?

问题描述

使用此Mathematica代码:

Take this Mathematica code:

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示例，由于需要无限的精度而失败:

Mathematica computes this without problem and gives the answer 0.010089328699390866240. I would like to be able to perform similar integrals but I don't have a copy of Mathematica. Just naively implementing it in scipy, for example, using a standard quadrature library fails sadly because f(x) and g(x) get arbitrarily close to 0. Here is an example in Python using standard quadrature which fails due to the infinite precision needed.:

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一样执行这种不正确的积分?我不介意使用哪种免费语言/库.

How, in code, can one perform this improper integral as Mathematica does? I don't mind which free language/library is used.

推荐答案

一个非常有趣的问题.

首先要注意的是被积物

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]上的积分.经过一番操作，您会发现

has a singularity at 0 that is integrable, and the integral over [0, infty] can be evaluated analytically. After some manipulation, you'll find

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，然后:

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, a project of mine). Cut off the domain at some large value, where the function is almost 0 anyway, then:

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

---

现在，对于其他一些事情，也许令人惊讶的是，它们不是工作得很好.

---

Now for some other things that, perhaps surprisingly, do not work so well.

看到类型为exp(-x) * f(x)的整数时，首先想到的是 quadpy (我的项目之一):

When seeing an integral of the type exp(-x) * f(x), the first thing that should come to mind is Gauss-Laguerre quadrature. For example with quadpy (one of my projects):

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:

which is a surprisingly bad approximation for the actual value despite the fact that we were using 100 integration points. This is due to the fact that the integrand cannot be approximated very well by polynomials, especially the terms log(x) and 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

这篇关于您如何像Mathematica一样执行这种不适当的积分?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！