Python scipy quad и nqaud дают разные ответы

Я не уверен, подходит ли это больше для математического обмена стеком или переполнения стека, но, поскольку математические доказательства кажутся мне приемлемыми, я подозреваю, что проблема в коде, и поэтому я опубликую его здесь:

Первая формула (формула 1) по определению является прямой:

и вот мой код для формулы 1:

import scipy.integrate as integrate
from scipy.integrate import nquad
from scipy.stats import norm
import math

def normalcdf(x):
    return (1+math.erf(x/math.sqrt(2)))/2

def normalpdf(x): 
    return math.exp(-x*x/2)/math.sqrt(2*math.pi)

def integrand12345(x1,x2,x3,x4,x5,theta1,theta2,theta3,theta4,theta5):
  return normalpdf(x1 - theta1) * normalpdf(x2 - theta2) * normalpdf(x3 - theta3) * normalpdf(x4 - theta4) * normalpdf(x5 - theta5)

def range_x1(theta1,theta2,theta3,theta4,theta5):
  return [-np.inf, np.inf]

def range_x2(x1,theta1,theta2,theta3,theta4,theta5):
  return [x1, np.inf]

def range_x3(x2,x1,theta1,theta2,theta3,theta4,theta5):
  return [x2, np.inf] 

def range_x4(x3,x2,x1,theta1,theta2,theta3,theta4,theta5):
  return [x3, np.inf] 

def range_x5(x4,x3,x2,x1,theta1,theta2,theta3,theta4,theta5):
  return [x4, np.inf]

def pi_12345(theta1,theta2,theta3,theta4,theta5):
  return (nquad(integrand12345, [range_x5, range_x4, range_x3, range_x2, range_x1], args=(theta1,theta2,theta3,theta4,theta5)))[0]

Вторая формула (формула 2) использует двойное интегрирование: и вот код формулы 2:

def integrandforpi_12(x1, x2, t1, t2, *theta): 
  prod = 1
  for ti in theta:
    prod = prod * (1 - normalcdf(x2 - ti))
  return prod * normalpdf(x1 - t1) * normalpdf(x2 - t2)

def range_x1(t1, t2, *theta):
  return [-np.inf, np.inf]

def range_x2(x1, t1, t2, *theta):
  return [x1, np.inf]

def pi_12(t1, t2, *theta):
  return (nquad(integrandforpi_12, [range_x2, range_x1], args=(t1, t2, *theta)))[0]

Моя третья формула (формула 3) основана на теореме Байеса:

и мой код для формулы 3:

def integrandforpi_i(xi, ti, *theta):
  prod = 1
  for t in theta:
    prod = prod * (1 - normalcdf(xi - t))
  return prod * normalpdf(xi - ti)

def pi_i(ti, *theta):
  return integrate.quad(integrandforpi_i, -np.inf, np.inf, args=(ti, *theta))[0]

Таким образом, pi_i вычисляет вероятность того, что X_i является наименьшим среди theta_i.

Однако когда я запускаю свой код, используя три разные формулы, все они дают разные ответы, и я понятия не имею, почему. Я не уверен, есть ли ошибка в выводе формулы или ошибка в моем коде. Любая помощь будет оценена по достоинству.

Вот пример:

t1,t2,t3,t4,t5 = 0.83720022,0.61704171,1.21121701,0,1.52334078

p12345 = pi_12345(t1,t2,t3,t4,t5)
p12354 = pi_12345(t1,t2,t3,t5,t4)
p12435 = pi_12345(t1,t2,t4,t3,t5)
p12453 = pi_12345(t1,t2,t4,t5,t3)
p12534 = pi_12345(t1,t2,t5,t3,t4)
p12543 = pi_12345(t1,t2,t5,t4,t3)

print('p12345=',p12345)
print('p12354=',p12354)
print('p12435=',p12435)
print('p12453=',p12453)
print('p12534=',p12534)
print('p12543=',p12543)

print('formula 1 gives', p12345+p12354+p12435+p12453+p12534+p12543)

print('formula 2 gives', pi_12(t1,t2,t3,t4,t5))

print('formula 3 gives', pi_i(t2,t3,t4,t5) * pi_i(t1,t2,t3,t4,t5))

и выход

p12345= 0.0027679276698449086
p12354= 0.008209750140618218
p12435= 0.0016182955786153714
p12453= 0.001921206801273682
p12534= 0.009675713474375739
p12543= 0.003904872716765966
formula 1 gives 0.028097766381493885
formula 2 gives 0.21897431741874426
formula 3 gives 0.0418669679120933

Реплика: Формула 1 очень медленная, ее запуск на моем бедном старом ноутбуке занимает около 3 часов. Формулы 2 и 3 являются мгновенными.