基于Python实现蒙特卡洛积分的技术问题求助

1. 解决C++随机种子与RAND_MAX的对应问题 #

• C++的srand(time(NULL))是用当前时间初始化随机数生成器，Python中可以结合time模块和NumPy的随机种子实现：

  import time
  import numpy as np
  
  # 用当前时间戳作为随机种子
  np.random.seed(int(time.time()))
  

• C中rand() / RAND_MAX是生成[0,1)区间的均匀分布随机数，NumPy里直接用np.random.uniform(0, 1, n_samples)就能生成相同效果的数组，无需手动处理RAND_MAX（NumPy内部已封装均匀分布的生成逻辑）。如果非要模拟C的rand()逻辑，可以用：

  # 模拟C++ rand()生成0到RAND_MAX的整数，再归一化
  RAND_MAX = 2**31 - 1  # 常见的RAND_MAX取值
  random_int = np.random.randint(0, RAND_MAX + 1, n_samples)
  random_float = random_int / RAND_MAX
  

  但实际蒙特卡洛积分中，直接用np.random.uniform更简洁高效。

2. 实现支持自定义被积函数的通用蒙特卡洛积分 #

我们可以写一个通用积分函数，接受被积函数、积分区间、采样数量作为参数，灵活计算各种积分（比如期望、方差相关的积分）。

完整实现代码 #

import numpy as np
import time

def monte_carlo_integrate(func, a, b, n_samples=1000000):
    # 可选：设置随机种子以复现结果，不需要可注释
    # np.random.seed(int(time.time()))
    
    # 生成[a, b]区间的均匀分布随机数
    x = np.random.uniform(a, b, n_samples)
    # 批量计算所有采样点的函数值
    f_vals = func(x)
    # 蒙特卡洛积分核心公式：(区间长度) × 函数值的平均值
    integral = (b - a) * np.mean(f_vals)
    return integral

# 示例1：计算普通函数积分，比如∫(0,1) x² dx
def test_func(x):
    return x**2

result = monte_carlo_integrate(test_func, 0, 1, 1000000)
print(f"∫x² dx from 0 to 1: {result:.6f}")  # 理论值1/3≈0.333333

# 示例2：计算随机变量的期望与方差（假设PDF为[0,2]上的均匀分布，PDF(x)=0.5）
def pdf(x):
    # 定义概率密度函数，仅在[0,2]区间有效
    return np.where((x >= 0) & (x <= 2), 0.5, 0.0)

# 计算E[x] = ∫x*pdf(x) dx
expected_x = monte_carlo_integrate(lambda x: x * pdf(x), 0, 2, 1000000)
# 计算E[x²] = ∫x²*pdf(x) dx
expected_x2 = monte_carlo_integrate(lambda x: x**2 * pdf(x), 0, 2, 1000000)
# 方差Var(x) = E[x²] - (E[x])²
variance_x = expected_x2 - expected_x**2

print(f"期望E[x]: {expected_x:.6f}")  # 理论值1.0
print(f"期望E[x²]: {expected_x2:.6f}")  # 理论值4/3≈1.333333
print(f"方差Var(x): {variance_x:.6f}")  # 理论值1/3≈0.333333

关键说明 #

• 自定义被积函数：通过func参数传入任意支持NumPy数组输入的函数，比如用lambda表达式快速定义x*pdf(x)或x²*pdf(x)。
• 向量化计算：NumPy的数组操作比单样本循环高效得多，批量生成随机数和计算函数值能大幅提升运行速度。
• 种子控制：如果需要复现计算结果，可固定种子；若每次运行要不同结果，用当前时间作为种子或不设置种子（NumPy会自动用系统熵初始化）。

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