如何通俗地理解「蒙特卡洛方法」，它解决问题的基本思路是什么，目前主要应用于哪些领域？

您所在的位置：网站首页 › 柱状图怎么表示数值范围和均值的关系 › 如何通俗地理解「蒙特卡洛方法」，它解决问题的基本思路是什么，目前主要应用于哪些领域？

如何通俗地理解「蒙特卡洛方法」，它解决问题的基本思路是什么，目前主要应用于哪些领域？

2023-05-12 09:41| 来源: 网络整理| 查看: 265

这篇回答节选自我在专栏《机器学习中的数学：概率统计》中的文章，一起来聊聊蒙特卡洛方法。

也欢迎关注我的知乎账号 @石溪，将持续发布机器学习数学基础及算法应用等方面的精彩内容。

1.大数定理的经典应用：蒙特卡罗方法

用大样本数据计算出来的频率去估计概率，这就是大数定理的本质，而大数定理思想的一个非常典型的应用就是蒙特卡罗方法。

蒙特卡罗方法，又叫统计模拟方法，名字很洋气，思想很粗暴，真的很管用。它使用随机数来进行场景的模拟或者过程的仿真，其思想核心就是通过模拟出来的大量样本集或者随机过程去近似我们想要研究的实际问题对象，这是一类非常重要的数值计算方法。

该方法的名字来源于世界著名的赌城蒙特卡罗。博彩和概率，二者相视一笑、不谋而合。这种方法最初应用于20世纪40年代美国的曼哈顿原子弹计划，如今在数据分析和机器学习领域中到处都有他的身影。

以下是蒙特卡罗方法的几类非常典型的应用，我们分别来举例介绍

1.近似计算不规则面积/体积/积分

2.模拟随机过程，预测随机过程可能性结果的区间范围

3.结合接受-拒绝采样，对分布的未知参数进行统计推断

要知道，蒙特卡洛方法不仅仅只是一种方法技巧，更是一种思考问题的方式。在这篇回答中，我们主要介绍第一个应用点，让大家进一步理解大数定理的原理和蒙特卡罗方法的要点。第二个和第三个应用如果大家有兴趣可以关注我的专栏。

2.利用蒙特卡罗方法计算不规则面积

蒙特卡罗方法可以近似计算出不规则图形的面积，特别对于那些难以用解析方法计算的图像，非常有效。

这里，我们利用蒙特卡罗方法来近似计算一个圆的面积，然后估计出 $\pi$ 的近似值，选择圆作为例子的原因不是因为圆无法通过解析法进行计算，而是因为大家对他的面积比较熟悉，方便我们进行结果的对比。

代码片段：

import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Circle from scipy.stats import uniform import seaborn seaborn.set() n = 100000 r = 1.0 o_x, o_y = (0., 0.) uniform_x = uniform(o_x-r,2*r).rvs(n) uniform_y = uniform(o_y-r,2*r).rvs(n) d_array = np.sqrt((uniform_x - o_x) ** 2 + (uniform_y - o_y) ** 2) res = sum(np.where(d_array 10]))) print("赌满全场且亏损的人数:{}".format(len(sample_df[sample_df['money']0])))

运行结果：

总轮数:100,总人数:100000 输光赌本提前出局的人数:31148 赌满全场且盈利的人数:44458 赌满全场且亏损的人数:23923

总轮数:1000,总人数:100000 输光赌本提前出局的人数:75154 赌满全场且盈利的人数:23441 赌满全场且亏损的人数:1386

总轮数:10000,总人数:100000 输光赌本提前出局的人数:91902 赌满全场且盈利的人数:8060 赌满全场且亏损的人数:38

从结果中不难发现，这种和庄家1:1的对赌，随着轮数的增加，基本上都破产被收割了。换句话说，哪怕庄家不出千，输赢概率各半，赌的越久，基本上都是输光破产走人，原因是什么？原因是庄家的资金量是无穷的。

3.2.模拟股价的变化

博彩的例子结束了，我们再举一个股票的例子，这个让人感觉都很心跳啊，在金融工程中，有下面这样一个公式，他利用目前的股价 S_t 去预测 $\Delta t$ 时间之后的股价 $S_{t+1}$ ：

$S_{t+1}=S_t+\hat{\mu}S_t\Delta t+\sigma S_t \epsilon\sqrt{\Delta t}$

这其中的参数我来解释一下：

$\hat{\mu}$ 表示股票收益率的期望值，这里我们设定为 $15\%$ ，即 $\hat{\mu}=0.15$

$\sigma$ 表示股票的波动率，这里设定为 $\sigma=0.2$

$\Delta t = \frac{T}{n}$ ，其中表示整数年份，表示在整个估算周期内，取的具体步数，就好比说为一年，如果取244，那么 $\Delta t$ 的粒度就是每个交易日了（一年有244个交易日）。

这里面似乎所有的参数都是确定的，唯独除了 $\epsilon$ 之外， $\epsilon$ 是一个服从标准正态分布的随机变量，这是这个 $\epsilon$ ，决定了每日的股价 S_i 是一个随机变量，而由股价构成的序列是一个随机过程。

我们同样的用蒙特卡罗方法，利用大样本来估计一下在目前股价为 S_0=10 的情况下，1年之后股价的概率分布情况。

代码片段：

import scipy import matplotlib.pyplot as plt import seaborn from math import sqrt seaborn.set() s0 = 10.0 T = 1.0 n = 244 * T mu = 0.15 sigma = 0.2 n_simulation = 10000 dt = T/n s_array = [] for i in range(n_simulation): s = s0 for j in range(int(n)): e = scipy.random.normal() s = s+mu*s*dt+sigma*s*e*sqrt(dt) s_array.append(s) plt.hist(s_array, bins=30, normed=True, edgecolor='k') plt.show()

运行结果：

图1.一年之后的股价分布情况

这是我们模拟了10000个样本经过上述随机过程，在一年之后股价的分布情况。

这里的核心就是：

e = scipy.random.normal() s = s+mu*s*dt+sigma*s*e*sqrt(dt)

每一轮 $\Delta t$ ，我们都生成一个服从标准正态分布的随机变量 $\epsilon$ ，不断通过递推公式 $S_{t+1}=S_t+\hat{\mu}S_t\Delta t+\sigma S_t \epsilon\sqrt{\Delta t}$ ，迭代出下一个时间点的股价，循环往复直到生成一年后的最终结果，这样就模拟出了一年过程中股价随机变量序列构成的随机过程。

我们采用蒙特卡罗方法，设置大样本量（这里设置10000个），最终迭代出10000个对应的一年后股价，然后用柱状图就能看出其总体分布特征。

3.3.股价变化曲线的过程展现

上面我们分析的是这10000个样本在一年之后最终股价的整体分布情况，实际上我们还有一个同样重要的过程可以进行监测和展现，那就是从 T_0 时刻起到1年后的这一段时间内，每隔 $\Delta t$ 时间间隔点由实时价格随机变量构成的序列，换句话说就是随机过程的整体展现。

这个在上面的基础上进一步完成，其实不难，就是我们不光要计算出股票最终的价格，还有记录下每个 $\Delta t$ 时间点的价格，并把他记录下来。

代码片段：

import scipy import matplotlib.pyplot as plt import seaborn from math import sqrt import numpy as np seaborn.set() s0 = 10.0 T = 1.0 n = 244 * T mu = 0.15 sigma = 0.2 n_simulation = 100 dt = T/n random_series = np.zeros(int(n), dtype=float) x = range(0, int(n)) for i in range(n_simulation): random_series[0] = s0 for j in range(1,int(n)): e = scipy.random.normal() random_series[j] = random_series[j-1]+mu*random_series[j-1]*dt+sigma*random_series[j-1]*e*sqrt(dt) plt.plot(x, random_series) plt.show()

运行结果：

图2.股价变化随机过程的展现

这里我们清晰的展现出了由244个 $\Delta t$ 时间点的股价数据所构成的序列，这是随着时间变化的随机过程。我们为了从整体把握他的分布特征，设定了100个样本，因此模拟出了100条价格曲线。

这种结果图表面上看起来比较凌乱，实际上可以从整体上发现许多端倪，例如股价在运行过程中的整体分布区间，上下界，集中程度等等，都可以有一个整体的把握。

因此我们不仅可以得到最终股价的分布，也可以知道股价运行变化的完整价格路径，这个价格路径代表了蒙特卡罗方法的精髓，通过这个价格路径的可视化呈现，让我们更加直观的从宏观上把握了随机过程。

4.蒙特卡洛方法与接受拒绝采样

4.1.运用蒙特卡洛方法的背景

随机近似方法的核心是蒙特卡洛方法，主要是用采样的方式来进行随机近似，来实现数值积分等目标。

例如我们要求函数 f(Z) 关于分布 p(Z|X) 的期望，从期望的定义我们知道，求期望的本质就是求积分： $E_{z|x}[f(z)]=\int_z p(z|x)f(z)dz$

但是这个积分往往是非常难求的，而随机近似的方法核心就是蒙特卡洛方法，他用采样的方法来实现数值积分的目标：

首先我们从原始分布 p(z|x) 中采出个样本点：

$z^{(1)},z^{(2)},z^{(3)},...,z^{(N)}\sim p(z|x)$

然后依据大数定律，用样本均值来模拟积分的真实结果：

$\frac{1}{N}\sum_{i=1}^Nf(z^{(i)})\approx \int_z p(z|x)f(z)dz$ 如果我们直接令这个函数 f(z)=z ，那么通过上述方法求出的就是分布 p(z|x) 的期望。

这种基于蒙特卡洛方法的近似方法，思想上非常简单直观，但是其中有一个问题目前看上去无法跨越，那就是如果 p(z|x) 比较复杂，我们似乎很难从中采样出服从概率分布的一组样本点。

于是，为了解决这个问题，就派生出了两种方案，一个称之为接受-拒绝采样，一个称之为是重要性采样，我们这里重点谈谈接受-拒绝采样。

他们都是基于一个事实，那就是目标分布 p(z) 无法直接采样，而我们引入一个提议分布 q(z) ，什么叫提议分布？换句话说就是我们好采样的分布，这个分布可以是任意的，比如均匀分布、高斯分布等等，怎么好采样就怎么来。

4.2.接受-拒绝采样的问题背景

这里我们就来重点介绍接受-拒绝采样方法，虽然在最终实际使用过程中，因为效率、维度等原因，不会实际采用他进行采样，但是他所折射出来的思想方法将直接指导马尔科夫链-蒙特卡洛方法的理解和应用，因此很有必要对其进行深入学习。

接下来看一个具体的例子，我们需要采样的目标分布就是下面这个复杂的概率密度函数：

$p(x)=(0.3e^{(-(x-0.3)^2)}+0.7e^{(-(x-2)^2/0.3)})/1.2113$

当然，写的更简洁、专业一点就是：

$p(x)=\frac{1}{1.2113}(0.3exp(-(z-0.3)^2)+0.7exp(-(z-2)^2/0.3)$

其中， 1.2113 是归一化参数，目的是为了让概率密度函数 p(x) 在整个取值域上的积分为，以满足归一化的概率基本定义。

4.3.采样方法和步骤

看上去这个概率密度函数非常复杂，那么针对这个概率密度函数的分布，我们如何进行采样，才能使得样本能够满足这个分布要求，从而支撑我们的后续的统计分析呢？这里，我们介绍最为基本、也最为经典的接受-拒绝采样方法，他的根本目标就是根据给定的概率密度函数 p(x) ，产生服从目标分布的样本集。

在开始采样之前，我们需要做一些准备工作。首先需要找到一个辅助的建议分布，我们将其记作是，他的概率密度函数为 g(y) ，这个分布在原理上可以是任意分布，既然是可以任意选择，那么我们就应该去选择一个易于处理、易于通过计算机程序去采样的分布类型，可以说是越简单越好。显然在这里，均匀分布、正态分布是最为合适的。

同时我们还需获取一个常数值，使得对于任意变量，都能够保证 $C*g(x)\ge f(x)$ 成立，当然为了提高采样的效率，选取的常数在满足条件的前提下，当然是越小越好了。

接下来就是真正采样的过程：

第一步：从建议分布中进行采样，获取一个采样样本；

第二步：从 [0,1] 的均匀分布中进行采样，获取一个采样样本；

第三步：判断，如果 $\frac{f(Y)}{C*g(Y)}\ge U$ ，那么我们就接受这个采样值，作为我们记录下来的采样值，如果不满足这个不等关系，我们就拒绝这个采样值，舍弃他，重新采样。这就是所谓的接受-拒绝的关键一步。

4.4.接受-拒绝采样的实践

我们在这个例子中，选择均值为1.4，方差为1.2的正态分布作为我们的建议分布：分布，取常数值为2.5。

我们先把目标分布的概率密度函数 p(x) 、建议分布的概率密度函数的倍： g(x)*C 绘制在一张图上进行比较：

代码片段：

import numpy as np import matplotlib.pyplot as plt from scipy.stats import norm import seaborn seaborn.set() # 目标采样分布的概率密度函数 def p(x): return (0.3 * np.exp(-(x - 0.3) ** 2) + 0.7 * np.exp(-(x - 2.) ** 2 / 0.3)) / 1.2113 # 建议分布G norm_rv = norm(loc=1.4, scale=1.2) # C值 C = 2.5 x = np.arange(-4., 6., 0.01) plt.plot(x, p(x), color='r', lw=5, label='p(x)') plt.plot(x, C*norm_rv.pdf(x), color='b', lw=5, label='C*g(x)') plt.legend() plt.show()

运行结果：

图1.建议分布和概率密度函数

从图中可以看出，我们所选取的建议分布的概率密度函数的倍，确实是恒大于目标分布的概率密度函数 p(x) ，因此满足我们接受-拒绝采样的前提条件。

那么接下来，我们就在此基础上，利用上文中所介绍的接受-拒绝采样算法的具体步骤，来进行实际的样本采样。

代码片段：

import numpy as np import matplotlib.pyplot as plt from scipy.stats import uniform, norm import seaborn seaborn.set() # 目标采样分布的概率密度函数 def p(x): return (0.3 * np.exp(-(x - 0.3) ** 2) + 0.7 * np.exp(-(x - 2.) ** 2 / 0.3)) / 1.2113 # 建议分布 norm_rv = norm(loc=1.4, scale=1.2) # C值 C = 2.5 uniform_rv = uniform(loc=0, scale=1) sample = [] for i in range(100000): Y = norm_rv.rvs(1)[0] U = uniform_rv.rvs(1)[0] if p(Y) >= U * C * norm_rv.pdf(Y): sample.append(Y) x = np.arange(-3., 5., 0.01) plt.gca().axes.set_xlim(-3, 5) plt.plot(x, p(x), color='r') plt.hist(sample, color='b', bins=150, normed=True, edgecolor='k') plt.show()

运行结果：

图2.实际采样效果图

我们通过接受-拒绝采样方法，进行了10万次的数据采样，最后把采样结果绘制成柱状图。通过运行结果我们发现，柱状图的形状几乎和原始目标分布的概率密度函数曲线完全拟合。由此通过实践证明了这个采样的方法确实是有效的，他高度近似的拟合了原始样本的分布特性。

有了高度近似于实际分布的采样值之后，我们就能够很容易的利用这组样本来对原始的分布进行统计特征的分析，比如均值、方差等数字特征的计算。

此内容节选自我的专栏《机器学习中的数学：概率统计》，前三节可免费试读，欢迎订阅！机器学习中的数学：概率统计_机器学习数学基础系列专栏-CSDN博客

当然还有《机器学习中的数学（全集）》系列专栏，欢迎大家阅读，配合食用，效果更佳~

机器学习中的数学（全集）_机器学习数学基础系列专栏-CSDN博客

有订阅的问题可咨询微信：zhangyumeng0422

【本文地址】

公司简介

联系我们