Scipy quad与nquad计算结果不一致问题排查求助

我不确定该问题更适合Math StackExchange还是Stack Overflow，但我认为数学推导无问题，怀疑是代码错误，故在此求助：

公式与对应代码 #

公式1：直接基于定义的五重积分 #

import scipy.integrate as integrate
from scipy.integrate import nquad
from scipy.stats import norm
import math
import numpy as np  # 注：原代码使用了np但未导入，已补充

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：基于二重积分的实现 #

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：基于贝叶斯定理的实现 #

pi_i用于计算$X_i$为所有$\theta_i$对应变量中最小值的概率：

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]

测试代码与运行输出 #

测试代码 #

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则可瞬间完成。
• 三个公式的计算结果完全不同，无法确定是公式推导错误还是代码存在问题，恳请帮忙排查。

内容的提问来源于stack exchange，提问作者Ishigami