Вопрос по calculus, scipy, python, numpy – Точность scipy.integrate.quad на больших числах

2

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

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

Ваш Ответ

2   ответа
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])
@DenisKorzhenkov: это заставляет интегратора выбирать функцию в этих точках. В противном случае, для большого интервала интеграции это будет выборка точекa + eps * (b - a) где eps это небольшое число --- но еслиb - a очень большой, он пропустит пик, близкий к х = 0. pv.
Это не. Интегратор делает выборку функции за пределами x> 100, но, конечно, это элементарный факт, что эта часть интеграла дает очень маленький вклад. pv.
Сравнивая результаты этого дела с тем изnp.quad(g, a=0., b=100)Похоже, что этот подход, по существу, устанавливает верхний предел равным 100, независимо от фактического пользовательского ввода. Конечно, это может быть хорошо для целей ОП. Stelios
@pv после прочтения quad docstring я не могу понять, как помогает твой совет. Точки 1 и 100 не являются точками разрыва Denis Korzhenkov
4

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

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

Одним из подходов, позволяющих избежать этой проблемы, является ограничение верхней границы интегрирования значением, выше которого числовая функция становится очень малой (и, следовательно, вносит незначительный вклад в интегральное значение). Из вашего фрагмента кода, значение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

Как правило, стоимость получения этой точности - увеличенное время вычислений.

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

mpmath также не безошибочен:mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**20]) ->2.20502636520112e-56, Дело в том, что численное интегрирование функций невозможно без некоторых условий «гладкости» - функция не должна иметь слишком резких «пиков» в интервале интегрирования. Когда интервал интегрирования очень большой, функцияexp(-x/2) очень "колючий", который вызывает проблемы. pv.
@pv. Действительно, спасибо за комментарий. Однако, если вы увеличите точность достаточно точно, такой проблемы не будет. Например, попробуйтеmp.mp.dps = 100 перед вызовомmp.quad Stelios
@Stelios ismpmath совместим сscipy, pandas и другие популярные пакеты? Denis Korzhenkov
Увеличивая точность, вы толкаете верхнюю границу вверх, попробуйте10**120, Это также увеличивает стоимость вычислений, что в этом случае не требуется. Проблема не в том, что значения функций настолько малы, что они находятся ниже диапазона с плавающей запятой, а в том, что функция, масштабируемая до интервала интегрирования, является очень колючей, что вводит в заблуждение оценку ошибки алгоритма интегрирования. pv.

Похожие вопросы