使用scipy.nquad计算联合概率异常，修正后遇类型错误求助

问题分析 #

1. 语法与引用错误：
   • 代码中多处使用spi.nquad，但未导入spi模块，正确调用应为直接使用nquad（已从scipy.integrate导入）。
   • 积分上下限中调用prob.function，但当前脚本未定义prob对象，应直接调用自定义的function函数，否则传递的是函数对象而非计算后的数值，触发TypeError: must be real number, not function。
2. 数值积分稳定性问题：
   直接对对数正态联合分布积分时，因变量范围设置过大（如10**8）易导致数值精度丢失，出现结果趋近于0的情况。通过变量替换转换为标准二元正态分布的积分，能大幅提升计算稳定性与准确性。

修正方案 #

1. 修复语法错误

将spi.nquad统一改为nquad，把所有prob.function替换为直接调用function。

2. 变量替换优化积分

将对数正态变量转换为标准正态变量：

• 令 ( Z_1 = \frac{\ln X_1 - \ln \mu_1}{\sigma_1} )，( Z_2 = \frac{\ln X_2 - \ln \mu_2}{\sigma_2} )
• 此时 ( Z_1, Z_2 ) 服从均值为0、方差为1、相关系数为correlation的二元正态分布，联合概率密度更简洁，积分计算更稳定。

修正后的代码 #

import numpy as np
from scipy.integrate import nquad
from scipy.stats import multivariate_normal

def get_z_threshold(x, mean, deviation):
    # 计算对应的标准正态变量阈值
    return (np.log(x) - np.log(mean)) / deviation

def Joint(x1, mean1, deviation1, sign1, x2, mean2, deviation2, sign2, correlation):
    # 转换为标准正态变量的阈值
    z1 = get_z_threshold(x1, mean1, deviation1)
    z2 = get_z_threshold(x2, mean2, deviation2)
    
    # 定义二元正态分布的参数
    mean = [0, 0]
    cov = [[1, correlation], [correlation, 1]]
    rv = multivariate_normal(mean, cov)
    
    # 根据符号确定积分区间（1对应<=阈值，2对应>=阈值）
    if sign1 == "1" and sign2 == "1":
        limits = [(-np.inf, z1), (-np.inf, z2)]
    elif sign1 == "2" and sign2 == "1":
        limits = [(z1, np.inf), (-np.inf, z2)]
    elif sign1 == "1" and sign2 == "2":
        limits = [(-np.inf, z1), (z2, np.inf)]
    elif sign1 == "2" and sign2 == "2":
        limits = [(z1, np.inf), (z2, np.inf)]
    else:
        raise ValueError("sign1和sign2只能是'1'或'2'")
    
    # 计算积分
    result, error = nquad(rv.pdf, limits)
    return result

额外说明 #

• 使用scipy.stats.multivariate_normal直接调用概率密度函数，避免手动编写公式出错。
• 用-np.inf和np.inf替代手动设置的10**-4和10**8，更符合理论积分区间，且nquad支持无穷区间积分。
• 变量替换后，积分的数值稳定性大幅提升，结果更接近Matlab的计算值。

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