如何加快复杂功能的集成过程？

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 IDE。

我需要集成具有 4 个或更多变量的函数。

Answer 1

有帮助的一件事是避免指定不可能实现的公差。双精度浮点数的 epsilon 约为

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

。