Python如何使用蒙特卡洛(Monte Carlo)方法对表中列数据进行积分计算

现有代码存在的问题 #

• 采样逻辑前提错误：random.choices直接随机采样原数据集的Y_k值，默认假设所有Y_k对应的X在[a,b]区间均匀分布，该逻辑仅在原始表的X是[a,b]上均匀采样的前提下成立；如果原始X分布不均，计算结果完全错误。
• 边界兼容问题：如果积分区间[a,b]和原始Y_k对应的X的取值范围不一致，直接采样原始Y_k会丢失区间外的函数信息，结果完全偏离真实值。
• 工程实现缺陷：没有把data1作为参数传入而是直接调用全局变量，复用性极低；硬编码采样量等于原始数据长度，无法灵活调整模拟样本数平衡精度和速度；用for循环累加效率远低于向量化运算。

正确实现方案 #

蒙特卡洛积分的核心是计算被积函数在积分区间的期望，乘以区间长度。根据原始表数据的采样特征，分两种场景实现：

场景1：原始data1的X是[a,b]区间的均匀采样结果 #

该场景下原代码的核心逻辑成立，仅需优化工程实现：

import numpy as np
import random

def MonteCarlo(a, b, data, sample_num=None):
    if sample_num is None:
        sample_num = len(data)
    y_samples = random.choices(data['Y_k'], k=sample_num)
    # 用均值替代循环累加，计算效率更高
    I = np.mean(y_samples) * (b - a)
    return I

调用示例：result = MonteCarlo(a=0, b=1, data=data1, sample_num=10000)

场景2：原始data1的X是任意非均匀采样结果 #

该场景需要先通过插值得到连续被积函数，再对积分区间做均匀采样计算：

from scipy.interpolate import interp1d
import numpy as np

def MonteCarlo(a, b, data, sample_num=10000):
    # 构造线性插值函数得到连续被积函数
    f = interp1d(data['X'], data['Y_k'], kind='linear', fill_value="extrapolate")
    # 在积分区间均匀采样X值
    x_samples = np.random.uniform(low=a, high=b, size=sample_num)
    # 计算Y值的期望乘以区间长度得到积分结果
    y_samples = f(x_samples)
    I = np.mean(y_samples) * (b - a)
    return I

调用示例：result = MonteCarlo(a=0, b=2, data=data1, sample_num=20000)

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