当前位置：首页 > news >正文

如何理解蒙特卡洛方法并用python进行模拟

news 2025/11/8 7:53:41

MC，即蒙特卡洛方法，是用来在概率空间，通过随机采样估算兴趣参数的后验分布。

这里先从为什么需要MC开始，然后尝试介绍MC的一般形式和接收-拒绝采样，并通过python对采样过程进行模拟和示例。

这里用到的示例图和演示代码整理和修改自网络资料。

1 采样限制

假设你需要对一维随机变量X进行采样，X的样本空间是{1, 2, 3}，且概率分别是{1/2, 1/4, 1/4}。

首先根据各离散取值的概率大小对[0, 1]区间进行等比例划分，如划分为[0, 0,5], [0.5, 0.75], [0.75, 1]这三个区间，再通过计算机产生[0, 1]之间的伪随机数，根据伪随机数的落点即可完成一次采样。

假如X是连续分布，概率密度是f(X)，就需要用到累积分布函数。

$P(X < t) = \int_{-\infty}^{t} f(x)dx$

即在[0, 1]间随机生成一个数a，求的使P(x<t)=a成立的t，t即可视为从X分布得到的一个采样结果。

1.1 可积和反函数

X能正常采样，依赖两个前提条件:

1）概率密度函数可积

2）累积分布函数有反函数

然而，大量现实世界的分布，大都难以满足这两个条件。

1.2 维度灾难

对于高维随机变量，比如 $\mathbb{R}^{50}$ ，每一维取100点，则一共需要取 $10^{100}$ 个点，已知宇宙基本粒子 $10^{87}$ 个。对于高维连续变量，差不多也是类似情况，同样面对维度灾难的问题。

2 蒙特卡洛

2.1 采样示例

比如对于下图，如何求解圆的面积，假设不用用公式求解，只能用近似的方法求解。

构思一种方法，在这个正方形内随机撒20个米粒，统计有多少个米粒落在圆内。

计算这部分比例，并乘以正方形的面积，即为对圆的面积的估计。

如下图所示，15/20，大约3/4左右。

以下是上述采样过程模拟程序

import random
import math
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from scipy.stats import norm
from matplotlib.patches import Ellipse, Circlefig2 = plt.figure(num=1, figsize=(5, 5))
axes2 = fig2.add_subplot(1, 1, 1)
c = Circle(xy=(0.5, 0.5), radius=0.5, alpha=0.5, color='red')  # 圆for _ in range(20):x = random.random()y = random.random()plt.scatter(x, y, s=30, color='black')axes2.add_patch(c)
plt.show()

这中模拟的方法在无法提前知道方程形状时有用，比如针对如下未知形状。

在一块长方形区域内随机抛掷点，MCMC方法能轻松计算出该区域近似值。

整个过程可以通过如下公式形式化表示

$Area = \int_{a}^{b} f(x)dx = \frac{b-a}{n}\sum_{i=1}^{n}f(x_i)$

即首先求高度的均值，然后乘以宽度(b-a)。任何一个不规则图形的面积都等价于一个矩形的面积。

另外一种理解，假设每一个xi，对应x轴1个基本单位，sum(f(xi))即在这种假设下长度为n的面积。

然而，实际上述不规则形状从左到右的长度不是1，而是b-a，所以变化后面积为

((b-a)/n) * sum(f(xi))。

2.2 MC一般形式

以上采样默认假定x在[a, b]之间是均匀分布，然而大部分情况x在[a, b]之间不是均匀分布。

可以做如下表示

$\int_{a}^{b} \frac{f(x)}{q(x)} q(x) dx$

这里把q(x)看做是x在区间[a, b]内的概率分布函数，

把前面的分数部分看做一个函数，然后在q(x)下抽取n个样本。

当n足够大时，可以用均值来近似。

$\sum_{i=1}^{n} \frac{f(x_i)}{q(x)}$

这个形式就是蒙特卡罗方法的一般形式。

2.3 均匀分布

均匀分布很容易采样的，在计算机生成[0, 1]之间的伪随机数序列，就可以看成是一种均匀分布。

常见概率分布，无论是连续还是离散，都可以基于Uniform(0, 1)样本生成。

例如，正态分布可以通过Box−Muller变换得到。

如果随机变量U1、U2 独立且 U1和U2～Uniform(0, 1)，则

$Z_0 = \sqrt{-2lnU_{1}} cos(2\pi U_{2})$

$Z_1 = \sqrt{-2lnU_{1}} sin(2\pi U_{2})$

则Z0和Z1独立且满足标准正态分布。

其它著名连续分布如指数分布、Gamma分布、t分布、F分布、Beta分布、Dirichlet 分布，都可以通过类似数学变换得到；离散分布通过均匀分布的到。python的numpy，scikit-learn库都有生成这些常用分布样本的函数。

当p(x)形式很复杂，或者p(x)是高维分布时，样本生成就可能困难。

譬如有如下的情况:

1）p(x)形式复杂

$p(x) =\frac{\hat{p}(x)}{\int \hat{p}(x)dx}$

$\hat{p}(x)$ 可以的到，但是位于分母的积分有可能很难得到解析解。

我们可以得到，但是下面的积分我们无法得到显示解（归一化系数很难计算）；

2) 高维分布

p(x, y)是二维分布函数，本身计算很困难，但条件分布p(x|y), p(y|x)的计算相对简单。

如果p(x)也是高维的，这种情形计算就更复杂了。

2.4 接收-拒绝采样

根据某已知分布的概率密度函数f(x)，产生服从此分布的样本X。

1）需要一个辅助建议分布(proposal distribution，已知其概率密度函）G来产生候选样本。

由于是由分布来产生候选样本，因此需要能够从G抽样。即必须能够生成服从此概率分布的样本 Y，比如均匀分布、正态分布。

2）需要另一个辅助的均匀分布 U(0, 1)。

3）计算一个常数值c，即满足c*g(x) >= f(x)不等式的c的最小值，希望c接近于1。

具体采样过程为：

步骤1：从建议分布G抽样，得到样本Y。

步骤2：从分布U(0, 1)抽样，得到样本U。

步骤3：如果U <= (f(Y))/(c*g(Y)) ，则令X=Y（接受Y），否则继续执行步骤1（拒绝）。

直观解释如下:

假设使用U(0, 1)来作为“proposal distribution” ，这样g(x)=1, x\in [0, 1]。

如下图所示，我们每次生成的两个样本U与Y，对应下图中矩形内的一点P(Y, Y*c*g(Y))，横轴坐标为Y，纵轴坐标为 Y*c*g(Y)，这个需要注意，因为通常假设Y为纵轴坐标。

接受条件U <= (f(Y))/(c*g(Y)) ，即U*c*g(Y) <= f(Y)，几何意义极为满足接收条件的点P在f(y)轴的下方，而不满足接收条件的点P在f(y)的上方。

在f(y)下方的点(o形状)满足接受条件，在f(y)上方的点(+形状)不满足接受条件。

2.5 python实现

假如目标概率密度函数是

$f(x) = \frac{(x-0.4)^{4}}{\int_{0}^{1} (x-0.4)^{4}dx}$

针对以上分布生成样本。

1）建议分布函数G选择U(0,1)，概率密度函数为g(x)=1, x \in [0, 1]

通常情况下，需要选择是的c最小的g函数。

2）辅助的均匀分布U(0, 1)

3）计算常数值，c需要选择是的c*g(x) >= f(x)即c>=f(x)的最小的c，相当于选择f(x)的最大值。

结合f(x)的形式，c = max((0.4**4)/sum, (0.6**4)/sum)，其中sum为f(x)的分母，由于sum是一个边界确定的积分，是一个常数项，后面代码实现时，可以省略掉。

python代码示例如下

import random
import math
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np%matplotlib inline
sns.set_style('darkgrid')
plt.rcParams['figure.figsize'] = (12, 8)def AceeptReject(split_val):global cglobal powerwhile True:x = random.uniform(0, 1)y = random.uniform(0, 1)if y*c <= math.pow(x - split_val, power):return xpower = 4
t = 0.4  
sum_ = (math.pow(1-t, power + 1) - math.pow(-t, power + 1)) / (power + 1)  #求积分
x = np.linspace(0, 1, 100)
#常数值c
c = 0.6**4/sum_
cc = [c for xi in x]
plt.plot(x, cc, '--',label='c*f(x)')
#目标概率密度函数的值f(x)
y = [math.pow(xi - t, power)/sum_ for xi in x]
plt.plot(x, y,label='f(x)')
#采样10000个点
samples = []
for  i in range(10000):samples.append(AceeptReject(t))
plt.hist(samples, bins=50, normed=True,label='sampling')
plt.legend()
plt.show()

采样示例图如下所示，绿色部分即为满足f(x)分布的样本。