Python实现Monte Carlo高斯积分时ValueError报错问题求解

错误原因 #

• 直接原因：monte_carlo_gaussian函数中调用less_than时传入了整个numpy数组x而非当前遍历的单个元素xi，数组与单值比较返回布尔数组，无法直接用于if判断，因此触发歧义错误。
• 其他逻辑错误：
  1. 高斯函数实现错误：你错误地将自然底数e放到了幂的底数中，标准正态分布的概率密度公式为$\frac{1}{\sqrt{2\pi}}e^{-\frac{x2}{2}}$，对应代码应先算系数再乘以exp的结果。
  2. y的抽样范围错误：标准正态分布的y值恒大于0且最大值约为0.3989，你从[a,b]（含负数）抽样y，会出现大量无意义的负y值，无法正确统计曲线下的点。
  3. 面积计算逻辑错误：蒙特卡洛积分的最终结果是「抽样矩形的总面积」乘以命中比例，矩形面积为(b-a)*y的上限，而非直接乘以b。
  4. 多余的区间判断：你生成x的时候已经指定了范围是[a,b]，不需要再额外判断x是否在区间内。

修正后的完整代码 #

import math
import numpy as np

# 标准正态分布概率密度函数
def gaussian(x):
    return (1 / math.sqrt(2 * math.pi)) * math.exp( (-x ** 2) / 2)

def under_curve(x, y):
    return 0 <= y <= gaussian(x)

def monte_carlo_gaussian(a, b, n):
    # 取标准正态分布y的最大值加少量冗余作为y的抽样上限
    y_max = 1 / math.sqrt(2 * math.pi) + 0.001
    # 分别在x、y范围内均匀抽样
    x = np.random.uniform(low=a, high=b, size=n)
    y = np.random.uniform(low=0, high=y_max, size=n)
    
    ndartsunder = 0
    # 遍历所有抽样点，统计落在曲线下的数量
    for xi, yi in zip(x, y):
        if under_curve(xi, yi):
            ndartsunder += 1
    # 计算矩形总面积乘以命中比例得到积分近似值
    rect_area = (b - a) * y_max
    prob_under_curve = ndartsunder / n
    return rect_area * prob_under_curve

# 测试调用，区间(-1.644854,1.644854)的标准正态积分理论值约为0.9
print(monte_carlo_gaussian(-1.644854, 1.644854, 10000))

实现思路说明 #

Monte Carlo计算定积分的核心逻辑是投点法：

1. 先确定待积分函数在积分区间[a,b]内的取值范围，生成一个刚好覆盖函数曲线的矩形区域
2. 在矩形区域内均匀随机投n个点
3. 统计落在函数曲线下方的点的数量，计算占总投点数的比例
4. 用矩形总面积乘以这个比例，就得到了积分的近似值，投点数量n越大结果越准确。

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