эффективное вычисление двойного интеграла

Я вычисляю следующий интеграл с помощью scipy:

from scipy.stats import norm
def integrand(y, x):
    # print "y: %s  x: %s" % (y,x)
    return (du(y)*measurment_outcome_belief(x, 3)(y))*fv_belief(item.mean, item.var)(x)
    return dblquad(
        integrand, norm.ppf(0.001, item.mean, item.var), 
        norm.ppf(0.999, item.mean, item.var),
        lambda x: norm.ppf(0.001, x, 3), 
        lambda x: norm.ppf(0.999, x, 3))[0]

У меня есть предположение о состоянии предмета (нормальное распределение) и измерения, зависящие от реальной стоимости предмета. (Также нормальное распределение). Используя этот интеграл, я вычисляю ценность информации (насколько полезно для измерения этого предмета).

Проблема в том, что вычисление этого интеграла занимает много времени. Есть ли более эффективный способ его вычисления (мне не нужна 100% точность), например, интеграция Монте-Карло или что-то подобное?

Я знаю, что в python есть skmonaco библиотека для интеграции monte carlo, но пределы интеграла должны быть числами, в отличие от scipy, внутренние интегральные пределы зависят от внешних (например, сверху

lambda x: norm.ppf(0.001, x, 3)

) Вот как вычисляется двойной интеграл с помощью skmonaco

>>> from skmonaco import mcquad
>>> mcquad(lambda x_y: x_y[0]*x_y[1], # integrand
...     xl=[0.,0.],xu=[1.,1.], # lower and upper limits of integration
...     npoints=100000 # number of points
...     )

Как видите, пределы внутреннего интеграла здесь не зависят от внешнего интеграла. Может ли кто-нибудь порекомендовать библиотеку или способ эффективного вычисления этого интеграла?


python scipy montecarlo
person Farseer    schedule 24.05.2015    source источник

Ответы (2)


arrow_upward
2
arrow_downward

Самый простой способ справиться с некубическими объемами интеграции в scikit-monaco - это переопределить вашу функцию интеграции так, чтобы она возвращала 0 за пределами области интеграции (см. этот раздел документации):

def modified_integrand(xs):
    x, y = xs
    if norm.ppf(0.001, x, 3) < y < norm.ppf(0.999, x, 3):
        return integrand(xs) # the actual integrand
    else:
        return 0.0

Как говорит Ами Тавори, это будет совершенно неэффективно, если большие участки вашего интеграционного региона равны нулю. Чтобы решить эту проблему, вы можете использовать алгоритмы MISER или VEGAS: оба этих алгоритма "изучают" форму интеграла по мере их выполнения, чтобы более эффективно распределять точки в интересующих областях.

Тем не менее, ваша область интеграции представляет собой повернутый прямоугольник:

import matplotlib.pyplot as plt
from scipy.stats import norm
import numpy as np

xs = np.linspace(-10, 10)
ys = np.linspace(-10, 10)

# Plot integration domain
# Red regions are in the domain of integration, blue
# regions are outside
plt.contourf(xs, ys, 
     [ [ 1 if norm.ppf(0.001, x, 3) < y < norm.ppf(0.999, x, 3) 
         else 0 for x in xs ] for y in ys ])

Домен интеграции

В этом случае вам будет гораздо лучше повернуть координаты интегрирования. Например, определение

r = x - y
R = (x + y)/2.0

Тогда ваше интегральное выражение:

def rotated_integrand(rs):
    R, r = rs
    x = R + r/2.0
    y = R - r/2.0
    return integrand(np.array([x,y]))

Пределы интеграции вдоль r теперь являются постоянными (-9.27..9.27 в приведенном вами примере). Предел интеграции по R по-прежнему (-inf, inf), поэтому вам нужно будет постепенно увеличивать область интеграции по R, пока ваш интеграл не сойдется. Я определенно рекомендую использовать для этого алгоритм MISER (mcmiser в scikit-monaco), а не единообразную выборку.

Наконец, судя по названиям функций, которые вы использовали, похоже, что вы выполняете некую форму байесовского обновления. В этом случае вы можете рассмотреть возможность использования библиотеки PyMC для цепей Маркова Монте-Карло, которая может быть более подходящей, чем библиотека интеграции MC общего назначения.

person Pascal Bugnion    schedule 24.05.2015

arrow_upward
1
arrow_downward

Интеграция методом Монте-Карло - хорошая идея, и ее не так уж сложно реализовать. Однако выборка точек в униформе сходится медленно по сравнению с другими методами, которые осуществляют итеративную выборку, и на каждом шаге выборку на основе результатов объема до этой точки.

Для этого проверьте рекурсивную стратифицированную выборку.

Вероятно, будет не очень сложно адаптироваться к вашему пользовательскому региону (который, по-видимому, не является гиперкубом).

К сожалению, я не думаю, что есть библиотеки Python, которые делают это из коробки.

person Ami Tavory    schedule 24.05.2015
comment
Sckiit-Monaco предлагает рекурсивную стратифицированную выборку из коробки: scikit-monaco.readthedocs.org/en/latest/ - person Pascal Bugnion; 24.05.2015