如何用Scipy.integrate实现复杂嵌套正态积分函数？

你的代码核心问题是被积函数的参数顺序与积分变量不匹配，导致外层积分错误地把v当成了积分变量，而非目标的u。以下是快速修正后的代码：

import scipy.integrate as integrate
from scipy.stats import norm
import numpy as np

def f(v, u, t1, ti, t2, *theta):
    prod = 1
    # 乘积项：遍历除ti外的theta元素（对应s≠1,i,2）
    for t in theta:
        if t == ti:
            continue
        prod *= (1 - norm.cdf(u + t2 - t))
    return norm.pdf(v) * norm.cdf(v + ti - t1) * prod

def g(u, t1, t2, *theta):
    S = 0
    # 求和项：遍历theta中的每个ti（对应i≠1,2）
    for ti in theta:
        # 对v积分，上限是u + t2 - ti
        integral_v = integrate.quad(f, -np.inf, u + t2 - ti, args=(u, t1, ti, t2, *theta))[0]
        S += integral_v
    return S * norm.pdf(u)

def P(t1, t2, *theta):
    # 对u积分，范围是(-inf, inf)
    return integrate.quad(g, -np.inf, np.inf, args=(t1, t2, *theta))[0]

关键修改点 #

1. 调整g函数的参数顺序：把积分变量u放在第一个位置，符合integrate.quad的传参规则（积分变量作为被积函数的第一个参数）。
2. 修正乘积项的逻辑：原代码的乘积项包含了当前的ti，不符合题目中“s≠1,i,2”的要求，新增判断跳过当前ti，确保乘积只针对除θ₁、θᵢ、θ₂外的元素。
3. 清理冗余导入：移除了没用的nquad和math导入，简化代码。

验证结果 #

运行以下测试代码：

print(P(0.2, 0.1, 0.3, 0.4))       # 输出：0.08856347190679764
print(P(0.2, 0.1, 0.4, 0.3, 0.5))  # 输出：0.06094233268837703

结果与你提供的正确示例完全一致。

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