如何解决使用 scipy.integrate.nquad 的双积分解决方案与integrate.dblquad 不匹配
下面代码中的第一个函数使用与scipy.integrate.dblquad
的双重积分来计算copula密度函数c*np.log(c)
的微分熵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)),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
维问题,但解决方案却不同?
解决方法
使用您提供的 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)),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
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。