矢量化`scipy.integrate.nquad`被积函数以与`qmc_quad`一起使用？

Question

我有使用

scipy.integrate.nquad

计算积分的代码，但我想使用

scipy.integrate.qmc_quad

执行积分以提高速度。

qmc_quad

要求对被积数进行向量化，以便可以在一次调用中在多个横坐标上对其进行求值，但我编写的函数一次仅在一个横坐标上求值。如何矢量化被积函数以与

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

内嵌注释来解释更改。

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)

请注意，标准误差非常高，但并不比

nquad

的误差估计更差。