Pfthb

I'm performing a one dimensional integral over an array as described here and here. As stated in those answers, I can not use scipy.integrate.quad to vectorize the integrals over the array because it employs an adaptive algorithm, which is why I use numpy.trapz below (I could also have used 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)

This works fine for a single dimension, but now I'd like to perform a double integral, replacing the c2 constant by a variable:

# 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

For what I could gather, the only available function is scipy.integrate.dblquad bu that means I could no longer apply the integral to an entire array in a single pass, and I would have to use a for loop which is considerably slower.

Is there any solution to this? I'm open to almost anything as long as it is reasonable performance-wise (I'm plugging this double-integral into an MCMC and it needs to be be evaluated millions of times)

ADD

This is my attempt with a 1-dimensional integral using scipy.integrate.quad inside a for loop (ie: one data value in the array at a time). The process is more than 50x slower than using np.trapz over the entire array.

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)

ADD 2

Testing the answer given below, it sort of works, with the caveat that at times it can fail super hard compared to dblquad (which is far slower but much more precise). I'm guessing this is related to the algorithm used by 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))

If I've understood your problem correctly, you can just call trapz twice:

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)