Правильно вычисляя двойной интеграл в питоне - PullRequest
Новая
       27

Правильно вычисляя двойной интеграл в питоне

0 голосов
/ 27 октября 2010

Я пытаюсь вычислить определенный двойной интеграл, используя scipy. Интегральная функция немного сложна, так как содержит некоторые вероятностные распределения, чтобы определить вес вероятности каждого значения x и y (например, смешанной модели). Следующий код оценивается как отрицательное число, но оно должно быть связано с [0,1]. Кроме того, для вычисления потребовалось около получаса.

У меня есть два вопроса.

1) Есть ли лучший способ вычислить этот интеграл?

2) Откуда это отрицательное значение? Большой вопрос для меня - как ускорить вычисления, так как я могу найти ошибку в своем коде, которая позже приведет к отрицательному результату самостоятельно. 

from scipy import stats
from scipy.integrate import dblquad
import itertools

p= [list whose entries are each different stats.beta(a,b) distributions]

def integrand(x,y):
        delta=x-y
        marg=0
        for distA,distB in itertools.permutations(p,2):
                first=distA.pdf(x)
                second=distB.pdf(y)
                weight1=0
                weight2=0
                for distC in p:
                        if distC == distA:
                                continue
                        w1=distC.cdf(x)-distC.cdf(y)
                        if weight1 == 0:
                                weight1=w1
                        else:
                                weight1=weight1*w1
                marg+=(first*weight1*second)
        I=delta*marg
        return I

expect=dblquad(integrand,0,1,lambda x: 0, lambda x: x)

Это, по сути, вопрос, почему ожидаемое значение максимального расстояния между двумя точками находится в векторе распределений. Пределы интегрирования: y ∊ [0, x] и x ∊ [0,1]. Это дало мне около -49, с оценочной погрешностью интеграла порядка 10e-10, поэтому это не должно происходить из-за метода интегрирования.

Я боролся с этим некоторое время и ценю любую помощь. Спасибо.

редактировать: исправлена ​​опечатка

Ответы [ 2 ]

1 голос
/ 07 ноября 2010

Существует несколько способов увеличить скорость вычислений.

  1. Вы можете использовать параметры epsabs и epsrel до dblquad, чтобы увеличить толерантность вашей интеграции.,Конечно, ваши результаты будут менее точными, но для отладки это нормально.

  2. Вы можете значительно сократить количество вычислений функций в integrand, переупорядочив код как (предупреждение,непроверенный код) 

    def integrand(x, y):
    marg = 0.0
    cdf = dict((id(distC), distC.cdf(x) - distC.cdf(y)) for distC in p)
    for distA in p:
        weight = numpy.prod(cdf[id(distC)]
                            for distC in p if distC is not distA)
        marg += weight * distA.pdf(x) * sum(
            distB.pdf(y) for distB in p if distB is not distA)
    return (x-y) * marg

    Но учтите, что у Python есть большие накладные расходы на вызовы функций, поэтому написание этого на чистом Python вас не зашло слишком далеко (для этого нужно использовать что-то вроде Cython проблема может немного помочь).

Я не знаю, почему интеграл становится отрицательным.Возможно, я мог бы сказать вам, если бы вы дали пример для p - это позволило бы нам на самом деле попробовать ваш код.

0 голосов
/ 27 октября 2010

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

Кстати: вы интегрируете PDF? Если да: вы уверены в своих лимитах интеграции?

...