scipy.integrate.quad точность на больших чисел

Question

Я пытаюсь вычислить такой интеграл (фактически CDF экспоненциального распределения с его PDF) с помощью scipy.integrate.quad():

import numpy as np 
from scipy.integrate import quad 

def g(x): 
    return .5 * np.exp(-.5 * x) 

print quad(g, a=0., b=np.inf) 
print quad(g, a=0., b=10**6) 
print quad(g, a=0., b=10**5) 
print quad(g, a=0., b=10**4) 

И результат выглядит следующим образом:

(1.0, 3.5807346295637055e-11) 
(0.0, 0.0) 
(3.881683817604194e-22, 7.717972744764185e-22) 
(1.0, 1.6059202674761255e-14) 

Все попытки использовать большой верхний предел интеграции дают неверный ответ, хотя использование np.inf решает проблему.

Аналогичный случай обсуждается в scipy issue #5428 at GitHub.

Что мне делать, чтобы избежать такой ошибки при интеграции других функций плотности?

Answer 1

Я считаю, что проблема связана с np.exp(-x), которая быстро становится очень маленькой, поскольку x увеличивается, что приводит к оценке как нулю из-за ограниченной числовой точности. Например, даже для x всего x=10**2*, np.exp(-x) оценивается в 3.72007597602e-44, тогда как x значения 10**3 или выше результат 0.

Я не знаю особенностей реализации quad, но он, вероятно, выполняет некоторую выборку функции, которая должна быть интегрирована в заданный диапазон интеграции. Для большого верхнего предела интегрирования большинство образцов np.exp(-x) оцениваются в ноль, поэтому значение интеграла недооценивается. (Заметим, что в этих случаях абсолютная погрешность по quad имеет тот же порядок, что и интегральное значение, которое является индикатором того, что последний является ненадежным.)

Один из способов избежать этой проблемы - ограничить верхнюю границу интеграции значение, выше которого числовая функция становится очень малой (и, следовательно, незначительно влияет на интегральное значение). Из вашего кода snipet значение 10**4 представляется хорошим выбором, однако значение 10**2 также дает точную оценку интегралу.

Другой подход, позволяющий избежать числовых проблем точности, заключается в использовании модуля, который выполняет вычисления в произвольной точной арифметике, такой как mpmath. Например, для x=10**5, mpmath оценивает exp(-x) следующим образом (с использованием нативного mpmath экспоненциальной функции)

import mpmath as mp 
print(mp.exp(-10**5)) 

3.56294956530937e-43430

Обратите внимание, как мала эта величина. При стандартной аппаратной числовой точности (используется numpy) это значение становится 0.

mpmath предлагает функцию интегрирования (mp.quad), которая может обеспечить точную оценку интеграла при любых значениях верхней интегральной границы.

import mpmath as mp 

print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, mp.inf])) 
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**13])) 
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**8])) 
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**5])) 

1.0 
0.999999650469474 
0.999999999996516 
0.999999999999997 

Мы также можем получить еще более точные оценки за счет увеличения точности, скажем, 50 десятичной точки (от 15, которая является стандартной точности)

mp.mp.dps = 50; 

print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, mp.inf])) 
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**13])) 
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**8])) 
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**5])) 

1.0 
0.99999999999999999999999999999999999999999829880262 
0.99999999999999999999999999999999999999999999997463 
0.99999999999999999999999999999999999999999999999998 

В общем, стоимость получения этой точности - это увеличенное время вычисления.

P.S .: Само собой разумеется, что если вы в состоянии оценить ваш интеграл аналитически в первую очередь (например, с помощью Sympy), вы можете забыть все вышеперечисленное.

Answer 2

Используйте points аргумент, чтобы сказать алгоритм, где поддержка вашей функции примерно есть:

import numpy as np 
from scipy.integrate import quad 

def g(x): 
    return .5 * np.exp(-.5 * x) 

print quad(g, a=0., b=10**3, points=[1, 100]) 
print quad(g, a=0., b=10**6, points=[1, 100]) 
print quad(g, a=0., b=10**9, points=[1, 100]) 
print quad(g, a=0., b=10**12, points=[1, 100])