Vectorize `scipy.integrate.nquad` integrand for use with `qmc_quad`?

Question

I have code that calculates an integral using scipy.integrate.nquad, but I want to perform the integration with scipy.integrate.qmc_quad for speed. qmc_quad requires the integrand to be vectorized such that it can be evaluated at many abscissae in a single call, but my function is written to evaluate the integrand at only one abscissa at a time. How do I vectorize the integrand to be compatible with qmc_quad?

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import nquad, qmc_quad


X = [11.3, 14.8, 7.6, 10.5, 12.7, 3.9, 11.2, 5.4, 8.5, 5.0, 4.4, 7.3, 2.9, 5.7, 6.2, 7.3, 3.3, 4.2, 5.5, 4.2]
N = len(X)
I = np.arange(1, N + 1)


def pdf_normal(m1, mn, s1, sn, i, n, x):
    m = (m1 * (n - i) / (n - 1)) + (mn * (i - 1) / (n - 1))
    s = (s1 * (n - i) / (n - 1)) + (sn * (i - 1) / (n - 1))
    y = (x - m) / s
    f = np.exp(-0.5 * y * y) / np.sqrt(2 * np.pi)
    pdf = f / s
    return pdf


def function_g_vec(args):
    [m1, mn, s1, sn] = args
    xi = X
    i = I
    n = N
    values = [pdf_normal(m1, mn, s1, sn, i_val, n, x_val) for i_val, x_val in zip(i, xi)]
    result = np.prod(values)
    return result


limits = [(-10, 30), (-10, 20), (0, 15), (0, 10)]
K = nquad(lambda *x: function_g_vec(x), ranges=limits)
print(K)
# K = 2.9335981345167932e-18
# error_K = 2.9282172140557324e-18

Answer 1

Comments inline to explain the changes.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import nquad, qmc_quad

# made X an array
X = np.asarray([11.3, 14.8, 7.6, 10.5, 12.7, 3.9, 11.2, 5.4, 8.5, 5.0, 
                4.4, 7.3, 2.9, 5.7, 6.2, 7.3, 3.3, 4.2, 5.5, 4.2])
N = len(X)
I = np.arange(1, N + 1)


def pdf_normal(m1, mn, s1, sn, i, n, x):
    m = (m1 * (n - i) / (n - 1)) + (mn * (i - 1) / (n - 1))
    s = (s1 * (n - i) / (n - 1)) + (sn * (i - 1) / (n - 1))
    y = (x - m) / s
    f = np.exp(-0.5 * y * y) / np.sqrt(2 * np.pi)
    pdf = f / s
    return pdf


def function_g_vec(args):
    # `args` is of shape `(4, M)`, where `M` is the number of
    # abscissae at which the integrand is being evaluated 
    [m1, mn, s1, sn] = args

    # Now, `m1`, `mn`, `s1`, and `sn` are of shape `(M,)`

    # align `xi` and `i` along zeroth axis, 
    # perpendicular to `m1`, `mn`, `s1`, and `sn`. 
    # That is, reshape them to be 2D arrays of shape `(N, 1)`,
    xi = X[:, np.newaxis]
    i = I[:, np.newaxis]
    n = N

    # one vectorized call to `pdf_normal` uses broadcasting to 
    # evaluate the function for every combination of row of xi/i 
    # and column of m1/mn/s1/sn.
    values = pdf_normal(m1, mn, s1, sn, i, n, xi)
    # values has shape `(N, M)`

    # take the product to reduce dimension `N`, resulting
    # in an array of shape `(M,)` 
    result = np.prod(values, axis=0)
    return result


limits = [(-10, 30), (-10, 20), (0, 15), (0, 10)]
a, b = np.asarray(limits).T
# ref = nquad(lambda *x: function_g_vec(x), ranges=limits)

res = qmc_quad(function_g_vec, a, b )
res
# QMCQuadResult(integral=2.885125308944705e-18, standard_error=1.0961630100611955e-18)

Note that the standard error is very high, but not worse than nquad's error estimate.