如何加快复杂函数的积分过程？

Question

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
import datetime
start = datetime.datetime.now()
plt.rcParams['axes.grid'] = True
XX=[None, 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]
I = [None,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]
N=20
    # Задаем G-ю функцию которая включает в себя произведения плотностей вероятности
    # каждого независимого результата Xi в которых МО и СКО изменяются по линейной зависимости
def G(M1, Mn, S1, Sn):
    def loc (M1, Mn, i, n):
        return (M1*(n-i)/(n-1)) + (Mn*(i-1)/(n-1)) 
    def scale (S1, Sn, i, n):
        return (S1*(n-i)/(n-1)) + (Sn*(i-1)/(n-1))
    def fnorm (x):
        return (np.exp((-x**2)/2))/(np.sqrt(2*np.pi))
    def y (M1, Mn, S1, Sn, i, n, x):
        return (x-loc(M1,Mn,i,n))/scale(S1,Sn,i,n)
    def F (M1, Mn, S1, Sn, i, n, x):
        return (fnorm(y(M1, Mn, S1, Sn, i, n, x)))/scale(S1, Sn, i, n)
    # Распаковка значений x и i из глобальных переменных
    x = XX[1:]
    i = I[1:]
    n=N
    # Вычисляем значение функции F для каждого x и i
    values=[F(M1, Mn, S1, Sn, i_val, n, x_val) for i_val, x_val in zip(i, x)]
    
    # Вычисляем произведение всех значений
    result = np.prod(values)
    
    return result
    # находим сомножитель К для получения общей ПВ оценок
options={'epsrel':1e-20, 'epsabs':1e-20}
K, errorK = integrate.nquad(G, ranges=[[0, 30],[0, 10],[0, 10],[0, 10]], opts=[options, options, options, options])
# K = 2.9613332457351404e-18    errorK = 9.999171231431291e-21
    #  формируем ПВ оценок
def pdf(Mn, S1, Sn,       M1):
    return (G(M1, Mn, S1, Sn) / K)
    #  строим график автономной ПВ оценок для параметра М1 (уменьшаем значения ошибок для оперативности)
def pdf_m1 (M1):
    return [(integrate.tplquad(pdf, 0, 10, 0, 10, 0, 10, args=(m,), epsabs=0.1, epsrel=0.1)[0]) for m in M1]
x = np.arange(4, 16, 0.2)
plt.plot(x, pdf_m1(x), 'k', lw=3)
plt.ylim(bottom=0)
plt.title('Fm1(m1)')
plt.show()
    # находим несмещенную оценку М1 и её ско sm1 (примерные значения по методу ММП М1 = 9.66)
def F1 (M1):
    return integrate.tplquad(pdf, 0, 10, 0, 10, 0, 10, args=(M1,))[0]
M1, _ = integrate.quad(lambda M1: M1 * F1(M1), 4, 16)
print(M1)
Dm1, _ = integrate.quad (lambda x: ((x-M1)**2) * F1(x), 4, 16)
sm1 = np.sqrt(Dm1)
print(sm1)
print("Время вычисления М1 и СКОм1:", datetime.datetime.now()-start)

问题在于效率。这个代码运行得太慢了，竟然要花费15个小时！这个任务主要是关于函数的积分。例子中只考虑了一个变量的积分，但实际上，你需要对PDF中指定的每一个变量都进行积分。

我在这个项目中使用Spyder这个开发环境。

我需要对4个或更多的变量进行积分。

Answer 1

总结一下：要避免设定不可能达到的精度要求，并考虑使用QMC积分。

---

有一点很重要，就是要避免设定那些根本无法实现的精度要求。双精度浮点数的精度大约是1e-16，所以你不可能希望达到1e-20的相对误差，但这会对你的运行时间产生很大的影响。

举个例子：

from time import perf_counter
from scipy.stats import multivariate_normal
from scipy.integrate import nquad
rng = np.random.default_rng(23925340854238936)
options={'epsrel':1e-20, 'epsabs':1e-20}

for d in range(1, 4):
    cov = rng.random((d, d))
    cov = cov @ cov.T
    dist = multivariate_normal(cov=cov)
    x = np.asarray([(0, 1)]*d)
    tic = perf_counter()
    res = nquad(lambda *x: dist.pdf(x), x)
    # res = nquad(lambda *x: dist.pdf(x), x, opts=[options]*d)
    toc = perf_counter()
    print(toc - tic)
    ref = dist.cdf(x[:, 1], lower_limit=x[:, 0])

# 0.0012877020001269557
# 0.08465272499984167
# 0.5809959000000617

数值积分的执行时间随着维度的增加而急剧增加，而严格的精度要求会让这个问题更加严重。当我们使用你的精度要求时，运行时间看起来像：

# 0.004921189000015147
# 0.4714967869999782
# 241.41230427999994

而这些时间几乎没有带来更好的准确性。如果不设定option，三维积分的误差是res[0] - ref = 4.969670549248573e-07，而使用了精度要求后，误差变成了-3.4931814399397076e-07。

实际上，对于你的四重积分，我们已经可以看到： # K = 2.9613332457351404e-18 errorK = 9.999171231431291e-21 第二个数字是一个绝对误差的估计。我不确定真正的误差是多少，但这意味着估计的相对误差大约是3e-3。如果你本来就没有达到很高的准确性，那不如使用scipy.integrate.qmc_quad，它可以在更短的时间内获得合理的准确性。

from time import perf_counter
import numpy as np
from scipy.integrate import nquad, qmc_quad
from scipy.stats import multivariate_normal
rng = np.random.default_rng(23925340854238936)

d = 4
cov = rng.random((d, d))
cov = cov @ cov.T
dist = multivariate_normal(cov=cov)
x = np.asarray([(0, 1)]*d)

ref = dist.cdf(x[:, 1], lower_limit=x[:, 0])

tic = perf_counter()
res = nquad(lambda *x: dist.pdf(x), x)
toc = perf_counter()
print(f"nquad    | Time: {toc - tic}, Error: {abs((res[0] - ref)/ref)})")

tic = perf_counter()
res = qmc_quad(lambda x: dist.pdf(x.T), *x.T)
toc = perf_counter()
print(f"qmc_quad | Time: {toc - tic}, Error: {abs((res[0] - ref)/ref)})")
# nquad    | Time: 13.483493682000244, Error: 2.8983787517192293e-05)
# qmc_quad | Time: 0.0281571279997479, Error: 0.000120747399615769)

如果你希望在相同的时间内从qmc_quad获得更好的准确性（代价是误差估计的准确性），可以减少n_estimates，同时将n_points增加相同的倍数。如果你能花更多时间来提高准确性——因为这已经是一个巨大的提速——那么只需增加n_points即可。