蒙特卡洛与 Python 的集成

Posted 2023-02-16

【中文标题】蒙特卡洛与 Python 的集成

【英文标题】：Monte Carlo Integration with Python

【发布时间】：2020-06-14 03:41:14

【问题描述】：

用蒙特卡洛积分估计以下积分：

我正在尝试对下面的问题进行蒙特卡洛积分，其中 p(x) 是均值为 1，方差为 2 的高斯分布。（见图）。

有人告诉我，一旦我们从正态分布中抽取样本，pdf 就会在积分中消失。请解释这个概念以及如何在 Python 中解决这个问题。下面是我的尝试。

def func(x):
return (math.exp(x))*x

mu = 1
sigma = sqrt(2)
N = 1000
areas = []
for i in range(N):
    xrand = np.zeros(N)

    for i in range (len(xrand)):
        xrand[i] = np.random.normal(mu, sigma)
        integral = 0.0

    for i in range (N):
        integral += func(xrand[i])/N

    answer = integral
    areas.append(answer)

plt.title("Distribution of areas calculated")
plt.hist(areas, 60, ec = 'black')
plt.xlabel("Areas")

integral

【参考方案1】：

蒙特卡洛积分是一种近似复积分的方法，无需计算其闭式解。为了回答您的问题，PDF 消失了，因为您需要做的就是 1）从指定的正态分布中采样一些随机值，2）计算被积函数中函数的值，以及 3）计算这些值的平均值。请注意，PDF 在计算中变得无关紧要；它仅在确保更频繁地采样更可能的值时才有意义。如果这样更直观，您可能会将其理解为采用加权平均值。

这是基于您的原始源代码的 Python 实现。

def func(x):
    return x * math.exp(x) 

def monte_carlo(n_sample, mu, sigma):
    val_lst = []
    for _ in range(n_sample):
        x = np.random.normal(mu, sigma)
        val_lst.append(func(x))
    return mean(val_lst)

您可以将func 更改为您选择的任何函数，以执行该函数的蒙特卡罗近似。如果给定不同的概率分布，您还可以编辑 monte_carlo 函数的参数。

这是一个可用于可视化蒙特卡洛近似逐渐收敛的函数。正如您所料，随着您增加n_sample 的值，这些值将随着更大的迭代而收敛，即。

MAX_SAMPLE = 200 # Adjust this value as you need
x = np.arange(MAX_SAMPLE)
y = [monte_carlo(i, 1, sqrt(2)) for i in x]
plt.plot(x, y)
plt.show()

结果图将显示收敛值，它是从定积分的封闭形式解计算的值的近似值。