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

Question

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

下面代码中的第二个函数尝试解决与上面相同的问题，但使用了多重积分求解器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 的函数给出-0.7127 的答案，而nquad 给出-0.5823 并且明显需要更长的时间。为什么即使我都设置为解决n=2 维问题，解决方案却不同？

Answer 1

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

-0.1931471805597395
0.17055845832017144,

不是-0.7127 和-0.5823。

第一个值（-0.1931471805597395）是正确的（你可以自己检查here）。

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

nquad 根据需要将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