Не могу понять такое поведение dblquad

Это случай, когда я думаю, что простой интеграл должен быть равен 1, но dblquad упрямо говорит 0,81.

Детали: функция createFun ниже возвращает реальную функцию. Это кусочно-линейная функция, предназначенная для функции плотности случайной величины, которая принимает значения в [-b, 1-a]

def createFun(a, b):
    def fIn(x):
        if x < -a:
            return (b+x)/(b-a)
        elif x > 1-b:
            return 1 - (x-(1-b))/(b-a)
        else:
            return 1
    return fIn

Теперь для конкретных значений a=0 и b=0,2 я создаю одну такую ​​функцию.

a = 0.0
b = 0.2
f1 = createFun(a, b)

Прежде всего, я убеждаюсь, что интеграл в [-b, 1-a] действительно равен 1:

print("Check it is a density, in [-b, 1-a]")
print(quad(f1, -b, 1-a))  # It is 1, as expected

Ладно пока. Теперь я определяю функцию g:RxR --> R как g(y,x) = f1(x) * f1(y).

def g(y, x):
    return f1(y) * f1(x)

Я ожидаю, что двойной интеграл от g в квадрате [-b, 1-a] x [-b, 1-a] должен быть произведением интегралов по каждой переменной (g факторизуется как функция x и a функция y), оба равны 1, поэтому я ожидал, что 1 * 1 = 1. Но код

print("Now, the double integral")
print(dblquad(g, -b, 1-a, lambda _:-b, lambda _:1-a))  # Should be 1, but isn't!

показывает 0,8161 с очень низкой ошибкой порядка 1E-10.

Моя математика неверна, и по какой-то причине этот интеграл не равен 1? Или это так, но я неправильно использую dblquad? Как же тогда быть?


python probability-density integral quad
person zeycus    schedule 26.09.2016    source источник
comment
Вы решили проблему?   -  person Severin Pappadeux    schedule 03.10.2016
comment
@Severin Нет, в конце концов я использовал Mathematica, чтобы получить формулу для интегральной функции в моей исходной задаче (немного более сложной, чем та, что в этом посте). Я до сих пор не понимаю, почему этот результат так далек от 1.0   -  person zeycus    schedule 04.10.2016

Ответы (1)


arrow_upward
0
arrow_downward

На самом деле это 1:

from scipy.integrate import quad, dblquad

a = 0.0
b = 0.2


def f(x):
    if x < -a:
        return (b + x) / (b - a)
    elif x > 1 - b:
        return 1 - (x - (1 - b)) / (b - a)
    else:
        return 1


print(quad(f, -b, 1 - a))


def g(y, x):
    return f(y) * f(x)


print(dblquad(g, -b, 1-a, lambda _:-b, lambda _:1-a))

(1.0, 1.1102230246251565e-15)
(1.0000000000000004, 2.55351295663786e-15)
person Nico Schlömer    schedule 03.11.2019
comment
С текущим scipy это так. Но это было не тогда. - person zeycus; 05.11.2019