scipy.integrate.quad точность по большим числам

Я пытаюсь вычислить такой интеграл (фактически 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 решает проблему.

Сибирский вопрос обсуждается в скупой проблеме № 5428 в GitHub .

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

    Я считаю, что проблема связана с тем, что np.exp(-x) быстро становится очень маленьким с ростом x , что приводит к оценке как нулю из-за ограниченной числовой точности. Например, даже для x , np.exp(-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 

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

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

    Используйте аргумент 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])