Python实现蒙特卡洛中点求和法求解1-8维积分咨询

你想要实现的是D维积分的蒙特卡洛中点求和法，对应公式如下：
#

现有代码的问题 #

你当前实现的是确定性均匀网格中点求积法，属于传统确定性数值积分方案，全程没有用到随机采样逻辑，因此和蒙特卡洛中点求和法的定义完全不符。
另外你当前的实现本质上仅处理了一维积分，仅因为你使用的被积函数可分离、积分结果刚好等于一维积分结果乘以维度D，才得到了正确的输出，并没有真正实现D维积分的计算逻辑。

蒙特卡洛中点求和法的核心逻辑 #

蒙特卡洛类积分方法的核心是通过随机采样的统计期望近似积分值，D维蒙特卡洛中点求和的执行流程如下：

• 确定D维积分域为[0,1]^D，该域的体积为1
• 将每个维度的区间[0,1]均匀划分为K个小段，每个小段的长度为1/K
• 随机生成N组D维的区间索引（每组索引对应D维空间中的一个网格单元）
• 取每个网格单元的中点作为采样点，计算所有采样点的被积函数值的平均值，即为积分的近似结果

正确实现代码 #

import numpy as np
from time import process_time

# 定义D维被积函数：输入x是形状为(N, D)的数组，N是采样数，D是维度
def f(x):
    return np.sum(x**2, axis=1)

a = 0
b = 1
N = 10**6  # 蒙特卡洛采样数，可根据精度需求调整

for D in range(1,9):
    t_start = process_time()
    # 随机生成N组D维的区间索引
    idx = np.random.randint(0, N, size=(N, D))
    # 计算每个采样点的D维坐标（网格中点）
    dx = (b - a) / N
    x_mid = (idx + 0.5) * dx
    # 计算积分近似值：函数值平均 * 积分域体积（D维单位立方体体积为1）
    mc_midpoint = np.mean(f(x_mid))
    t_stop = process_time()
    
    print(f"维度D={D}时CPU耗时：{t_stop-t_start:.6f}秒")
    print(f"积分近似结果：{mc_midpoint:.8f}\n")

补充说明 #

• 蒙特卡洛积分的结果存在随机误差，多次运行结果会有小幅波动，误差量级约为1/√N，你可以通过增大采样数N来降低误差
• 如果要达到你原确定性方法的高精度结果，可将N调整到1e7以上，误差会降到1e-4以内
• 你原代码中引入的scipy随机数、sympy积分、matplotlib、scipy.stats.norm等依赖在当前积分场景中没有用到，可以按需移除

内容的提问来源于stack exchange，提问作者J.Snowden