0

在这种情况下，我认为简单的积分应该是 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？那应该怎么做呢？

4

1 回答 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)
于 2019-11-03T19:11:20.853 回答