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

Question

下面代码中的第一个函数使用与scipy.integrate.dblquad的二重积分来计算copula密度函数[=的微分熵c*np.log(c) 15=]，它有一个依赖参数，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 的函数给出了 -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 passes lambda to the function as required. The lambda is programmed in such a way that it accepts arbitrary number of arguments, so it happily takes it and puts it in the list of powers and sums. Hence you do not get, e.g. 1/u**t+1/v**t -1 but 1/u**t+1/v**t + 1/t**t -1. The function call just does not match the intended function use. If you instead would write us[0]**() + us[1]**() - 1 it works.

修改后的代码如下：

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

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

Double integral solution using scipy.integrate.nquad doesn't match integrate.dblquad

python

math

mismatch

integral

scipy