python - 对数组执行双积分
问题描述
我正在这里和这里描述的数组上执行一维积分。正如这些答案中所述,我不能用来scipy.integrate.quad
对数组上的积分进行矢量化,因为它采用了自适应算法,这就是我numpy.trapz
在下面使用的原因(我也可以使用scipy.integrate.simps
or scipy.integrate.romb
)
import numpy as np
# Define some random data
Ndata = 500
data1 = np.random.uniform(.1, .8, Ndata)
data2 = np.random.uniform(.01, .2, Ndata)
# 1D function to integrate
def distFunc(x, c1=1., c2=.1):
B1 = ((data1 - (1. / x)) / data2)**2
B2 = ((x - c1) / c2)**2
f = np.exp(-.5 * B1) * np.exp(-.5 * B2) / c2
return f
# Values in x to evaluate the integral.
x = np.linspace(.1, 10, 100).reshape(-1, 1)
# Integral in x for each of the Ndata values defined above.
int_exp = np.trapz(distFunc(x), x, axis=0)
这适用于单维,但现在我想执行双积分,用c2
变量替换常量:
# 2D function to integrate
def distFunc(x, y, c1=1.):
B1 = ((data1 - (1. / x)) / data2)**2
B2 = ((x - c1) / y)**2
f = np.exp(-.5 * B1) * np.exp(-.5 * B2) / y
return f
对于我所能收集到的信息,唯一可用的函数是scipy.integrate.dblquad
bu,这意味着我不能再一次将积分应用于整个数组,我必须使用一个for
速度相当慢的循环。
有什么解决办法吗?只要性能合理,我几乎对任何事情都持开放态度(我将这个双积分插入到 MCMC 中,它需要被评估数百万次)
添加
scipy.integrate.quad
这是我在循环内使用一维积分的尝试for
(即:一次在数组中的一个数据值)。np.trapz
该过程比在整个阵列上使用要慢 50 倍以上。
import numpy as np
from scipy.integrate import quad
# Define some random data
Ndata = 500
data1 = np.random.uniform(.1, .8, Ndata)
data2 = np.random.uniform(.01, .2, Ndata)
# Function to integrate
def distFunc2(x, data1_i, data2_i, c1=1., c2=.1):
B1 = ((data1_i - (1. / x)) / data2_i)**2
B2 = ((x - c1) / c2)**2
f = np.exp(-.5 * B1) * np.exp(-.5 * B2) / c2
return f
s = t.time()
int_exp = np.zeros(Ndata)
for i in range(Ndata):
int_exp[i] = quad(distFunc2, .1, 10., args=(data1[i], data2[i]))[0]
print(t.time() - s)
添加 2
测试下面给出的答案,它有点工作,但需要注意的是,与dblquad
(慢得多但更精确)相比,它有时会失败得非常严重。我猜这与np.trapz
.
# Define some random data
Ndata = 10
data1 = np.random.uniform(.1, 10., Ndata)
data2 = np.random.uniform(.1, .2, Ndata)
c1 = .1
print(integ_dblquad(c1, data1, data2))
print(integ_trapz(c1, data1, data2))
def integ_dblquad(c1, data1, data2):
def distFunc(y, x, d1_i, d2_i, c1):
B1 = ((d1_i - (1. / x)) / d2_i)**2
B2 = ((x - c1) / y)**2
return (np.exp(-.5 * B1) / d2_i) * np.exp(-.5 * B2) / y
int_exp = np.zeros(data1.size)
for i in range(data1.size):
int_exp[i] = dblquad(
distFunc, .1, 10., lambda x: 0, lambda x: 5.,
args=(data1[i], data2[i], c1))[0]
return np.sum(np.log(int_exp))
def integ_trapz(c1, data1, data2):
def distFunc2d(x, y):
B1 = ((data1 - (1. / x)) / data2)**2
B2 = ((x - c1) / y)**2
return (np.exp(-.5 * B1) / data2) * np.exp(-.5 * B2) / y
# Values in x to evaluate the integral.
x = np.linspace(.1, 10, 1000)
y = np.linspace(.1, 5., 1000)
# Integral in x for each of the Ndata values defined above.
int_exp2d = np.trapz(np.trapz(distFunc2d(x[:, np.newaxis], y[:, np.newaxis, np.newaxis]), y, axis=0), x, axis=0)
return np.sum(np.log(int_exp2d))
解决方案
如果我正确理解了您的问题,您可以拨打trapz
两次:
import numpy as np
# Define some random data
Ndata = 500
data1 = np.random.uniform(.1, .8, Ndata)
data2 = np.random.uniform(.01, .2, Ndata)
# 1D function to integrate
def distFunc(x, c1=1., c2=.1):
B1 = ((data1 - (1. / x)) / data2)**2
B2 = ((x - c1) / c2)**2
f = np.exp(-.5 * B1) * np.exp(-.5 * B2) / c2
return f
def distFunc2d(x, y, c1=1.):
B1 = ((data1 - (1. / x)) / data2)**2
B2 = ((x - c1) / y)**2
f = np.exp(-.5 * B1) * np.exp(-.5 * B2) / y
return f
# Values in x to evaluate the integral.
x = np.linspace(.1, 10, 100)
y = np.linspace(.1, 10, 100)
# Integral in x for each of the Ndata values defined above.
int_exp = np.trapz(distFunc(x[:,np.newaxis]), x, axis=0)
int_exp2d = np.trapz(np.trapz(distFunc2d(x[:,np.newaxis],y[:,np.newaxis,np.newaxis]), y, axis=0), x, axis=0)
推荐阅读
- node.js - Dialogflow V2 并将请求正文中的参数传递给 webhook
- python - 如何使用 numpy 将一个小矩阵添加到一个大矩阵中?
- docusignapi - Docusign API - 锁定字段
- c++ - Whats 表示符号“2”。在 C++ 中?
- html - 如何为元素设置正确的高度?
- css - 在另一个元素被转换后,如何使 div 占用 100%
- azure - 如何在 azure vm 中托管 neo4j 服务并从任何地方(本地计算机)访问它?
- python - python圆形函数
- hibernate - Hibernate 中 LockMode.PESSIMISTIC_WRITE 和 LockMode.UPGRADE_NOWAIT 的区别
- python-3.x - 如何在定义的代码中编辑变量?