使用scipy.integrate.nquad的双积分解与integrate.dblquad不匹配

Question

下面的代码中的第一个函数使用与scipy.integrate.dblquad的双重积分来计算copula密度函数c的微分熵c*np.log(c)，它有一个依赖参数theta，通常是正的。

​

​

下面代码中的第二个函数试图解决与上面相同的问题，但使用多重积分求解器scipy.integrate.nquad

from scipy import integrate
import numpy as np

def dblquad_(theta):
    "Double integration"
    c = lambda v, u: ((1+theta)*(u*v)**(-1-theta)) * (u**(-theta)+v**(-theta)-1)**(-1/theta-2)
    return -integrate.dblquad(
        lambda u,v: c(v,u)*np.log(c(v,u)), 
        0, 1, lambda u: 0, lambda u: 1
        )[0]

def nquad_(n,theta):
    "Multiple integration"
    c = lambda *us: ((1+theta)*np.prod(us)**(-1-theta)) * (np.sum(np.power(us,-theta))-1)**(-1/theta-2)
    return -integrate.nquad(
        func   = lambda *us : c(*us)*np.log(c(*us)), 
        ranges = [(0,1) for i in range(n)],
        args   = (theta,) 
        )[0] 

n=2
theta = 1
print(dblquad_(theta))
print(nquad_(n,theta))

dblquad-based函数给出了-0.7127的答案，而nquad给出了-0.5823，并且明显地花费了更长的时间。为什么解决方案是不同的，即使我已经将两者都设置为解决n=2-dimensional问题？

Answer 1

使用您提供的n和theta值，代码的输出如下：

-0.1931471805597395
0.17055845832017144,

不是-0.7127和-0.5823。

第一个值(-0.1931471805597395)是正确的(您可以自己检查它这里 )。

nquad_中的问题在于theta参数的处理。谢谢@mikuszefski提供了解释；为了清楚起见，我在这里复制它：

nquad根据需要将lambda传递给函数。lambda的编程方式是它接受任意数量的参数，所以它很高兴地接受它，并将它放在幂和列表中。因此，您没有得到，例如1/u**t+1/v**t -1但是1/u**t+1/v**t + 1/t**t -1。函数调用与预期的函数使用不匹配。如果您想编写us[0]**() + us[1]**() - 1，它就能工作。

以下是修改后的代码：

from scipy import integrate
import numpy as np

def dblquad_(theta):
    "Double integration"
    c = lambda v, u: ((1+theta)*(u*v)**(-1-theta)) * (u**(-theta)+v**(-theta)-1)**(-1/theta-2)
    return -integrate.dblquad(
        lambda u,v: c(v,u)*np.log(c(v,u)),
        0, 1, lambda u: 0, lambda u: 1
        )[0]

def nquad_(n,theta):
    "Multiple integration"
    c = lambda *us: ((1+theta)*np.prod((us[0], us[1]))**(-1-theta)) * (np.sum(np.power((us[0], us[1]),-theta))-1)**(-1/theta-2)
    return -integrate.nquad(
        func   = lambda *us : c(*us)*np.log(c(*us)),
        ranges = [(0,1) for i in range(n)],
        args=(theta,)
        )[0]

n=2
theta = 1
print(dblquad_(theta))
print(nquad_(n,theta))

输出：

-0.1931471805597395
-0.1931471805597395