使用Scipy积分计算scipy stats分布期望E(X)的代码问题修复

• 对数正态分布仅在x>0时有非零概率密度，x≤0时pdf恒为0，理论上全区间积分等价于从0到+∞的积分。你使用的参数下，分布的取值集中在1e7~1e8量级，quad默认从-∞到+∞积分时采样点集中在原点附近，无法覆盖分布的有效取值范围，因此返回了错误的极小值。
• 你定义的被积函数本身逻辑没有错误，不需要修改被积函数，仅需要调整积分区间即可。

from scipy.integrate import quad
import math
import numpy as np
from scipy import stats

lognorm_mu = 17.79
lognorm_sigma = 0.81
scale = math.exp(lognorm_mu)

# 计算∫xf(x)dx（期望）的被积函数
def integrand_expectation(x):
    return stats.lognorm.pdf(x, s=lognorm_sigma, scale=scale) * x

# 计算∫x²f(x)dx（二阶原点矩）的被积函数
def integrand_second_moment(x):
    return stats.lognorm.pdf(x, s=lognorm_sigma, scale=scale) * x ** 2

# 积分计算，下限设为0
expectation = quad(integrand_expectation, 0, np.inf)[0]
second_moment = quad(integrand_second_moment, 0, np.inf)[0]

# 和理论公式对比
print("积分得到的期望：", expectation)
print("理论期望：", math.exp(lognorm_mu + lognorm_sigma ** 2 / 2))
print("积分得到的二阶原点矩：", second_moment)
print("理论二阶原点矩：", math.exp(2 * lognorm_mu + 2 * lognorm_sigma ** 2))

运行后可以看到积分得到的期望和二阶矩与理论公式计算结果完全一致，误差在数值计算精度范围内。
补充：如果不需要自定义积分逻辑，也可以直接调用scipy统计分布的内置方法快速得到结果：stats.lognorm.mean(s=lognorm_sigma, scale=scale) 直接返回期望，stats.lognorm.moment(2, s=lognorm_sigma, scale=scale) 直接返回二阶原点矩。

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