如何在Python蒙特卡洛积分计算中实现多进程加速优化

初始问题与代码 #

我编写了一个使用蒙特卡洛模拟计算给定区间内指定函数积分的Python程序，功能本身可以正常运行，但当为了获得更高精度使用更大的N值时，运行速度极慢。我打算尝试用多进程来提速，却不知道该如何实现，以下是我最初编写的代码：

from scipy import random
import numpy as np
import matplotlib.pyplot as plt
from multiprocessing import Process
import os

# 目标：用蒙特卡洛模拟近似计算函数f(x)在积分上下限a到b的积分值

# 积分区间
a = 0
b = np.pi

# 待积分函数
def f(x):
    return np.sin(x)

N = 10000
areas = []


def mcIntegrate():
    
    for i in range(N):
        
        # 生成区间范围内的随机数数组
        xrand = random.uniform(a, b, N)
        
        # 对每个随机数代入函数求和
        integral = 0.0
        for i in range(N):
            integral += f(xrand[i])
        
        # 用区间差除以样本量缩放得到积分近似值
        ans = integral * ((b - a) / float(N))
        
        # 将近似值加入结果列表
        areas.append(ans)


if __name__ == "__main__":

    processes = []
    numProcesses = os.cpu_count()

    for i in range(numProcesses):
        process = Process(target=mcIntegrate)
        processes.append(process)

    for process in processes:
        process.start()

    for process in processes:
        process.start()

    # 绘制近似值分布图
    plt.title("Distribution of Approximated Integrals")
    plt.hist(areas, bins=30, ec='black')
    plt.xlabel("Areas")
    plt.show()

请问我可以获得该实现的相关帮助吗？ #

优化后实现 #

我采纳了建议，使用了multiprocessing.Pool，同时改用NumPy向量化运算减少了部分运算量。N取值为10000时，运行时长从约5分钟缩短到了现在的6秒左右，以下是最终实现代码：

import scipy
import numpy as np
import matplotlib.pyplot as plt
import multiprocessing
import os

# 目标：用蒙特卡洛模拟近似计算函数f在积分上下限a到b的积分值

a = 0       # 积分下限
b = np.pi   # 积分上限
f = np.sin  # 待积分函数
N = 10000   # 样本量

def mcIntegrate(p):
    xrand = scipy.random.uniform(a, b, N)       # 生成区间范围内的随机数数组
    integral = np.sum(f(xrand))                 # 对每个随机数代入函数求和
    approx = integral * ((b - a) / float(N))    # 用区间差除以样本量缩放得到积分近似值
    return approx


if __name__ == "__main__":
    
    # 并行运行N次模拟，将结果存入数组
    with multiprocessing.Pool(os.cpu_count()) as pool:
        areas = pool.map(mcIntegrate, range(N))

    # 绘制近似值分布图
    plt.title("Distribution of Approximated Integrals")
    plt.hist(areas, bins=30, ec='black')
    plt.xlabel("Areas")
    plt.show()

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