蒙特卡洛方法求积分

蒙特卡洛法求积分

如果我们要求解一个函数 $h(x)$ 的积分，我们可以将积分分解成如下的形式：

$\int h(x)dx = \int f(x)p(x)dx$

于是，函数h(x)的积分可以表示为函数f(x)关于概率密度函数p(x)的数学期望，我们只需要取f(x) = h(x)/p(x)。

求解积分举例

求解积分$\int_0^1x^2dx$

解：另 p(x) = 1，表示0-1上的均匀分布，$f(x) = x^2$

求解代码：

import numpy as np

def calc_res(values):
    res = 0
    for one_value in values:
        res += one_value * one_value
    res = res / len(values)
    return res

for i in [100, 1000, 10000, 100000, 1000000]:
    res = np.random.uniform(0, 1, i)
    print(calc_res(res))

结果：
0.3097443410375912
0.33154677188274595
0.3363434070471736
0.33119288078816694
0.33372059640184654