I trying to calculate this function:

Which is the result of multiples integrals.

Functions F_A_, F_B_, F_C_, F_D_ and F_FPS_ is a truncated normal function, contained between intervals seted. I tried it with scipy.nquad, but I can't get in the end because it's too slow.

My last step is this below, in which I tried to consider Fs as a normal functions (not truncated) and adjust the limits of integration to get what I want.

%%time

muA, sigmaA = 303, 1
muB, sigmaB = 517, 2
muC, sigmaC = 1524, 1
muD, sigmaD = 1784, 2
muF, sigmaF = 30, 1

mu1 = muA
sigma1 = sigmaA
mu2 = muB
sigma2 = sigmaB
mu3 = muC
sigma3 = sigmaC
mu4 = muD
sigma4 = sigmaD
mu5 = muF
sigma5 = sigmaF

#overhead
factor = (sigma3*np.sqrt(2*np.pi))**(-1)*\
         (sigma1*np.sqrt(2*np.pi))**(-1)*\
         (sigma5*np.sqrt(2*np.pi))**(-1)*\
         (sigma1*np.sqrt(2*np.pi))**(-1)*\
         (sigma2*np.sqrt(2*np.pi))**(-1)*\
         (sigma2*np.sqrt(2*np.pi))**(-1)*\
         (sigma4*np.sqrt(2*np.pi))**(-1)*\
         (sigma4*np.sqrt(2*np.pi))**(-1)*\
         (sigma3*np.sqrt(2*np.pi))**(-1)

def pdfV(w, p, m, k, h, q, g, y, n, \
         mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, mu5, sigma5):
    return \
           np.exp(-0.5*((w-mu3)/sigma3)**2-0.5*((p-mu1)/sigma1)**2-0.5*((m-mu5)/sigma5)**2-0.5*((k-mu1)/sigma1)**2-0.5*((h-mu2)/sigma2)**2-0.5*((q+p-mu2)/sigma2)**2-0.5*((g+h-mu4)/sigma4)**2-0.5*((y/q+w-mu4)/sigma4)**2-0.5*((n/g+k-mu3)/sigma3)**2)

def lim1(p, m, k, h, q, g, y, n, mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, mu5, sigma5):
    return [mu3 - sigma3, mu3 + sigma3]

def lim2(m, k, h, q,  g, y, n, mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, mu5, sigma5):
    return [mu1 - sigma1, mu1 + sigma1]

def lim3(k, h, q,  g, y, n, mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, mu5, sigma5):
    return [mu5 - sigma5, mu5 + sigma5]

def lim4(h, q,  g, y, n, mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, mu5, sigma5):
    return [mu1 - sigma1, mu1 + sigma1]

def lim5(q,  g, y, n, mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, mu5, sigma5):
    return [mu2 - sigma2, mu2 + sigma2]

def lim6(g, y, n, mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, mu5, sigma5):
    return [mu2-mu1 - sigma2, mu2-mu1 + sigma2]

def lim7(y, n, mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, mu5, sigma5):
    return [mu4-mu2 - sigma2, mu4-mu2 + sigma2]

def lim8(n, mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, mu5, sigma5):
    return [(mu2-mu1)*(mu4-mu3) - sigma2, (mu2-mu1)*(mu4-mu3) + sigma2]

def lim9(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, mu5, sigma5):
    return [(mu4-mu2)*(mu3-mu1) - sigma2, (mu4-mu2)*(mu3-mu1) + sigma2]

options={'epsabs':0.01,'epsrel':0.01}

result = integrate.nquad(pdfV, [lim1, lim2, lim3, lim4, lim5, lim6, lim7, lim8, lim9], args=(muA, sigmaA, muB, sigmaB, muC, sigmaC, muD, sigmaD, muF, sigmaF), opts=options)

factor*result[0], result[1]

But I can't get the result because the algorithm didn't reach the end, rs.