主页 / coding / 问题 / 78874684
Accepted
Ishigami
Asked: 2024-08-15 18:33:15 +0800 CST

Python scipy quad 和 nqaud 给出不同的答案

  • 772

我不确定这是否更适合数学 stackexchange 或 stack overflow，但由于数学证明对我来说看起来没问题，我怀疑问题出在代码上，因此我将它发布在这里：

在此处输入图片描述

第一个公式（公式 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 中最小的概率。

但是，当我使用 3 个不同的公式运行我的代码时，它们都给出了不同的答案，我不知道为什么。我不确定我的公式推导是否存在缺陷，或者我的代码是否存在错误。任何帮助都将不胜感激。

以下是一个例子：

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 是即时的。

python

2 个回答

  1. Best Answer

    由于nquad范围列表被反转，被积函数参数X列表也应该反转：

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

    和

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

    那么我的结果就变成：

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

    opts此外，通过将参数传递给nquad：我还显著减少了公式 1 的计算时间：opts={'epsabs':0.001}

    • 3

  2. 这并不能完全回答这个问题，但了解这些是否正确可能会有所帮助。为此，我们进行了模拟：

    import numpy as np
from scipy import stats
rng = np.random.default_rng(43258234589253)

ti = np.asarray([0.83720022, 0.61704171, 1.21121701, 0, 1.52334078])

n_trials = 1000000
n_obs = len(ti)
desired_ranks = np.arange(1, n_obs+1)

Xi = rng.standard_normal(size=(n_trials, n_obs)) + ti
ranks = stats.rankdata(Xi, axis=1)
in_order = np.all(ranks[:, :2] == desired_ranks[:2], axis=1)
p = np.mean(in_order)
print(p)  # 0.044489

    如果此代码是正确的（如果我误解了请纠正我，我会更正），那么所有积分都不正确。

    实际上，我还模拟了 P(x1 是所有 xi 中最小的) * P(x2 是 x2、x3、x4、x5 中最小的)。它给出的结果与积分 3 相同，但与完整模拟或积分 1/2（校正后）不同。在您的推导中，您将 P(x2 是第二小的 | x1 是最小的) 视为等同于 P(x2 是 x2、x3、x4、x5 中最小的)。那部分可能不正确，因为它忽略了 x2 > x1 的事实。

    • 1

相关问题