不能理解dblquad的这种行为

问题描述 投票:0回答:1

在这种情况下,我认为一个简单的整数应为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)

我期望在平方[-b,1-a] x [-b,1-a]中g的双积分应该是每个变量中积分的乘积(g是x的函数)和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 integral probability-density quad
1个回答
0
投票

实际上是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)
© www.soinside.com 2019 - 2024. All rights reserved.