Как вы можете выполнить этот неправильный интеграл, как Mathematica? - PullRequest
Новая
       25

Как вы можете выполнить этот неправильный интеграл, как Mathematica?

2 голосов
/ 20 октября 2019

Возьмите этот код 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? Я не против, какой свободный язык / библиотека используется.

Ответы [ 2 ]

5 голосов
/ 21 октября 2019

В 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

Но разверните журнал (f / g), чтобы войти (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) может дать вам больше десятичных знаков.

3 голосов
/ 21 октября 2019

Очень интересная проблема.

Сначала отметим, что подынтегральное выражение 

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