Python scipy quad and nqaud giving different answers

I am not sure if this is more suitable for math stackexchange or stack overflow, but since the mathematical proofs look alright to me, I suspect its the code that's the issue and hence I will post it here:

The first formula (formula 1) is straight by definition:

and here's my code for formula 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]

The second formula (formula 2) is using double integration: and here is the code for formula 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]

My third formula (formula 3) is based on Bayes' theorem:

and my code for formula 3 is:

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]

So pi_i computes the probability that X_i is the smallest among the theta_i’s.

However, when I run my code using the 3 different formulas, they all give different answers and I have no idea why. I am not sure if there is a flaw in my derivation of the formula or if there's a mistake in my code. Any help would be appreciated.

Here is an example:

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))

and the output is

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

Remark: Formula 1 is extremely slow, it takes about 3 hours to run on my poor old laptop. Formula 2 and 3 are instant.

Since the nquad ranges list is reversed, the integrand X argument list should be reversed too:

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

and

def integrandforpi_12(x2, x1, t1, t2, *theta):

Then my results become:

formula 1 gives 0.04444581969881939
formula 2 gives 0.04444465903549115
formula 3 gives 0.0418669679120933

Also I significantly decreased the calculation time of formula 1 passing the opts argument to nquad: opts={'epsabs':0.001}