蒙特卡罗（Monte Carlo）方法——从数学原理到实际案例

Monte Carlo方法简介

Monte Carlo方法是一种应用随机数来进行计算机模拟的方法，通过对所研究系统进行随机观察抽样并对样本值进行统计分析，来得到所研究系统的某些参数。

Monte Carlo方法也称为统计模拟方法，是一种以概率论与数理统计理论为指导的数值计算方法，应用随机数（或伪随机数）来进行模拟试验，适用于一些解析法难以求解甚至不可能求解的问题。

Monte Carlo方法的基本思想是，首先确定求解问题，建立一个概率模型或随机过程，并使问题的参数或数字特征（期望、方差、协方差、矩）等于问题的解，然后通过对模型或过程进行观察或抽样试验来计算这些参数或数字特征，最后给出待求解的近似值。简单来说就是“不断抽样、逐渐逼近”。

Monte Carlo方法的基本原理是，用试验中事件发生的频率来逼近事件发生的概率，当样本容量足够大时，可以认为事件发生的频率等于其概率。

Monte Carlo方法的理论依据是大数定律和中心极限定理。

Monte Carlo方法

Monte Carlo方法适用于本身就具有随机性的问题或者能够转化为概率模型进行求解的确定性问题。 Monte Carlo方法是一种强大的数值计算工具，它通过随机抽样的方式来近似求解复杂的数学问题，尤其适用于难以用传统数学方法求解的问题。

Monte Carlo方法的步骤

①描述或构造概率过程（建立随机试验模型）

根据实际问题的特点，构造简单而又便于实现的概率统计模型，使所求的解恰好是所求问题的概率分布或数学期望。对于本身就具有随机性质的问题，主要是描述和模拟其概率过程；对于本身不具备随即性质的确定性问题（比如计算定积分），就必须先构造一个具有随机性质概率过程，并使它的参数正好是待求问题的解（比如通过投飞镖试验来求解不规则图形的面积） ，即将不具有随机性质的问题转化为随机性质的问题。

确定性数学问题：首先建立一个与待求解有关的概率模型，使待求解就是所建立模型的概率分布或数学期望。然后对这个模型进行随机抽样观察，即产生随机变量。最后用算术平均值作为待求解的近似估计值。比如，求积分、矩阵求逆、解方程（组）、偏微分方程边值问题、计算微分方程特征值等。

随机性问题：一般采用直接模拟，根据实际物理情况的概率法则，用计算机进行抽样试验。

②利用概率分布进行抽样（进行随机试验）

通过计算机产生已知概率分布的随机变量。常用的概率分布包括均匀分布、正态分布、指数分布、泊松分布、卡方分布等。这一步是Monte Carlo方法的核心，因为它涉及到如何准确地模拟随机事件。

在上一步构造了概率模型后，由于各种概率模型都是由各种概率分布构成的，因此产生已知概率分布的随机变量就成了实现Monte Carlo方法模拟试验的基本手段，这也是Monte Carlo方法被称为随机抽样的原因。实现Monte Carlo模拟的基本工具是随机数，以均匀分布为例，随机数可以看成具有均匀分布的随机变量，随机数序列就是具有均匀分布的总体的一个简单子样（服从均匀分布的相互独立的随机变量序列）。因此，在概率分布上抽样的过程实际上就是产生随机数的过程。

随机数的产生有物理方法和数学方法。在计算机上通过物理方法产生的是真正的随机数，但是代价较高。我们一般使用数学递推公式来产生伪随机数序列。大量试验已经证明，伪随机数与随机数具有相近的性质，因此可以作为真正的随机数来使用。

③建立估计量

构造随机概率模型并从中抽样后，需要确定一个随机变量作为所求问题的解，一般称之为无偏估计。通常，使用多次随机抽样结果的算术平均值作为解的近似值。

④重复实验分析结果

为减少估计方差并提高准确性，需要进行大量重复试验。最后对试验结果进行分析，计算估计值的置信区间和误差，评估解的可靠性。

Monte Carlo方法的模拟手段

模拟是Monte Carlo方法的重要手段，一般通过物理或者数学模型来近似类比，模拟现实系统的演变过程，来寻找其规律。模拟的基本思想是建立一个试验模型，这个模型包含了所研究系统的主要特点，通过运行该模型来获取需要的信息。模拟的实现方法有物理模拟和数学模拟。

物理模拟：通过与实际系统相近的实物系统模拟。比如，导弹发射试验、抛掷硬币试验等。物理模拟通常代价较大，难以重复进行，甚至无法实现。

数学模拟：数学模拟一般通过计算机实现，也叫做计算机模拟。计算机模拟可以反复进行，但是在实际问题中，建模需要进行种种简化和假设。

Monte Carlo方法的收敛性与误差

Monte Carlo方法作为一种计算方法，有必要讨论其收敛性以及误差。

收敛性

蒙特卡罗方法的收敛性可以由大数定律给出。

蒙特卡罗方法是由随机变量X的简单子样X1,X2,...,Xn的算术平均值来作为所求解的近似值。

由大数定律可知，如果X1,X2,...,Xn独立同分布，且具有有限期望值（E(X)<∞），那么

也就是说，当子样数N足够大时，随机变量X的算术平均值以概率1收敛到它的数学期望值。蒙特卡罗方法的收敛是概率意义上的收敛，虽然不能断言蒙特卡罗方法的误差不会超过某个值，但是能指出它的误差以接近1的概率不超过某个界限。

误差

蒙特卡罗方法的近似值与真值的误差可以由中心极限定理给出。

中心极限定理指出，如果X1,X2,...,Xn独立同分布，并且具有有限非零的方差

其中，f(x)是X的密度函数。那么

当N充分大时，有如下的近似表达

我们把α叫做置信度，1-α叫做置信水平。该近似表明，不等式

近似地以概率1-α成立，且误差收敛速度的阶为

一般将蒙特克罗方法的误差ε定义为

误差ε表达式中的正态差λ与置信度α是一一对应的，根据问题要求确定出置信水平之后，通过查阅正态分布表，就可以确定出λ的值。下面给出常用的λ与α的数值对照表。

当α=0.5时，误差超过ε的概率和误差小于ε的概率相等，即α=1-α=0.5，此时称ε为概然误差。

蒙特卡罗方法的误差是一种概率误差，这和其他数值计算方法是有区别的。并且误差表达式中的σ是未知的，必须使用它的估计值来代替，估计值计算方法如下

效率

通过误差表达式可以看出，在给定置信度α后，误差ε由σ和N决定。想要减小误差ε，要么增大试验次数N，要么减小方差σ²。并且，在标准差σ（标准差也叫均方差）固定的情况下，要想把精度提升一个数量级，试验次数要增加两个数量级（因为N和ε是平方根倒数的关系）。如果通过减小均方差σ，例如σ减小一半，那么误差ε也能减小一半，这相当于试验次数N会增大4倍。

一般来说，降低方差往往会使得观察一个子样的时间增加。在固定时间内，使观察的样本数减少。所以，该方法的优劣，需要由方差和观察一个子样的费用（计算机耗费的时间）综合来衡量。这就是蒙特卡罗方法中效率的概念，效率定义为 σ²·c ，其中c是观察一个子样的平均费用， σ²·c 越小，效率越高。假设在固定误差ε，和产生一个随机变量X的平均费用c不变的条件下，如果均方差σ缩小为原来的1/10，那么试验次数N可以减小为原来的1/100。我们将误差表达式变形可得，N=(λσ/ε)² ，那么总的费用 N·c=(λσ/ε)²·c=(λ/ε)²·σ²c。由此可见要想提高蒙特卡罗方法的效率，应该尽可能的降低 σ²c 的值。

由误差定义可知，在给定置信水平的情况下，蒙特卡罗方法的收敛速度为 O(N^-0.5)，这个值和问题本身的维数无关，问题本身维数的变化，只会引起抽样时间以及估计量计算时间的变化，不会对误差造成影响。也就是说，使用蒙特卡罗方法时，抽取的子样总数N与维数s无关，维数的增加只会引起计算量的增加，不会影响误差。这个特点使得蒙特卡罗方法对解决多维问题有很好的适应性。而一般的数值计算方法，比如计算定积分，计算时间会随着维数的幂次方而增加，并且由于分点数与维数的幂次方成正比，需要占用一定的存储空间。

大数定律和中心极限定理

在完美世界中，我们可以准确的回答任何问题。但是在现实世界中，要精确的计算出某些概率是非常困难甚至不可能的。不过，对于许多现实问题，一般不需要知道确切答案，只要得到一定范围内的答案就足够了。

极限定理是概率论中最重要的理论结果，而极限定理中最重要的便是大数定律和中心极限定理。通常，当随机变量序列的平均值在某种条件下收敛到数学期望（均值）的时候，就是我们所说的大数定律。当大量随机变量之和的分布在某种条件下逼近于正态分布时，就称之为中心极限定理。

大数定律

大数定律也叫做大数法则，通俗来讲，就是在样本数量很大的时候，样本均值和真实值充分接近。大数定律与中心极限定理结合在一起将成为现代概率统计科学的基石。

切比雪夫不等式

首先介绍两个不等式

我们知道，如果知道了随机变量的概率分布，我们可以直接计算出概率值，但是有时候概率分布是未知的。此时，马尔可夫不等式和切比雪夫不等式表明：当我们只知道随机变量的期望，或者同时知道期望和方差时，可以导出概率的上界。

收敛类型

大数定律有强大数定律和弱大数定律两种类型，二者主要考察的是独立随机变量之和的性质，而在很多情况下，这些性质都与极限有关。所以，有必要给出以下几种收敛类型。

在后面，我们会经常用到依概率收敛的概念，实际上依概率收敛的含义是：随机序列距X的偏差不小于任一给定值的可能性会随着n的增大而越来越小。

弱大数定律（辛钦大数定律）

如果在一个固定分布中进行独立抽样，那么样本均值任意接近于总体均值的概率趋近于1。简言之，样本均值依概率收敛于期望（弱大数定律的收敛是一种依概率收敛）。所以弱大数定律也可以表示为下面形式。

强大数定律

强大数定律表明，独立同分布的随机变量序列，前n个观察值的平均值以概率为1地收敛到分布的平均值。强大数定律中的收敛是指几乎必然收敛。

强大数定律和弱大数定律的区别：

弱大数律只能保证对充分大的 n * ，随机变量(X1+…+Xn * )/n * 靠近u。但它不能保证对一切n>n*，(X1+·+Xn)/n也一定在u的附近。这样，|(X1+·+Xn)/n-u|可以无限多次离开0（尽管出现较大偏离的频率不会很高)。而强大数律能保证这种情况不会出现，强大数律能够以概率为1地保证，

只可能出现有限次。

伯努利大数定律

伯努利大数定律表明， 事件发生的频率依概率收敛于事件发生的概率。该定理以严格的数学形式表达了频率的稳定性，也就是说当n很大时，事件发生的频率与总体概率有较大偏差的可能性很小。（概率是频率的稳定值）。

比如最简单的伯努利试验抛硬币，理论上正面朝上和反面朝上的概率都是0.5。依据伯努利大数定律，我们重复的抛掷一枚均匀的硬币，当抛掷次数n足够大时，正面和反面出现的次数因该是一样的，都是0.5n，并且次数n越大，这个值越接近。

中心极限定理

中心极限定理用粗略的语言来说，大量的独立随机变量之和的分布近似为正态分布。中心极限定理为计算独立随机变量和的有关概率提供了理论依据，同时也解释了现实世界中许多实际分布的频率曲线为钟形曲线（正态概率密度）的原因。中心极限定理揭示了测量误差近似正态分布这一重要的统计规律，因此也被称为误差频率定律。

中心极限定理的陈述

首先来看什么是正态分布。正态分布是概率中最重要的分布，正态分布也称为高斯分布和钟形曲线。

中心极限定理有很多版本，各版本的区别在于假设条件和最终收敛类型的不同。下面给出几种常见形式的中心极限定理。

中心极限定理为统计学提供了正态分布的广泛应用可能性，因为许多统计方法都是基于正态分布的性质来推导的。在实际应用中，中心极限定理使得我们能够通过对样本数据的分析来推断总体的参数，即使我们不知道总体的具体分布形态。这对于数据分析和决策制定有着重要的实际意义。中心极限定理还解释了为什么在实际中可以使用样本均值来估计总体均值，并且可以通过样本均值的分布来估计总体均值的可靠性。中心极限定理表明，在一定条件下，大量独立随机变量的和或平均值会趋向于正态分布，无论这些随机变量本身服从何种分布。这一发现极大地简化了概率计算和统计分析的过程。

中心极限定理与高斯白噪声

在使用KF、EKF、UKF进行滤波估计时，我们一般会假设传感器的观测噪声为高斯白噪声，并且在很多场合都会假设噪声是服从正态分布的。那么为什么要假设噪声为正态分布，这种假设的理论依据又是什么呢？

首先，高斯白噪声在不同时刻下的观测值相互独立。在实际系统中，一般存在众多噪声源（比如通过传感器测量角度时，影响测量结果的因素有很多），实际上这些噪声源大多满足相互独立的假设（比如声呐每次测量角度时，测量结果都是相互独立的），当噪声源的数量足够多且每个噪声源对于总体的贡献可忽略不计时，根据中心极限定理可知，这些噪声源的累加结果服从正态分布（无论这些噪声本身服从什么分布，只要他们相互独立，那么他们的累加和也就是最终的噪声就是服从正态分布的）。 简单来说，根据中心极限定理，多个独立随机变量的和趋近于高斯分布，无论这些变量本身的分布如何。在实际环境中，噪声往往是多种不同来源的随机变量的叠加，因此其整体分布趋向于高斯分布。

中心极限定理的应用

在微积分中，我们可以对一个函数在某个区间上积分，并得到积分值。但实际上，求积分的前提是，可以找到一个给定函数的原函数的完美解析表达式。在很多时候，这是一件非常困难甚至不可能实现的事。在很多时候为了计算积分，需要借助无穷级数展开（误差函数）来实现。于是，就有了中心极限定理最重要的应用之一：蒙特卡罗积分。

假设有一个面积为R的正方形，在正方形内有一个不规则的区域A，我们现在要求这个区域A的面积。如果我们可以用函数来较好的描述区域A，那么很自然的我们会想到通过对函数积分来求A的面积，此时需要知道原函数，即：若A可以描述为

那么A的面积可以通过以下积分来求

如果函数f和g的原函数难以得到，那么积分法就很难实现了。这时就可以采用蒙特卡罗法，其可以描述为：在正方形R内随机取点，那么A的面积应该接近于落在A中的点数占总共抽取的点数的比例乘以R。以投飞镖为例，向R中投掷N次飞镖，那么

这就是蒙特卡罗法。再比如下面图形求阴影面积，使用蒙特卡罗将是一个很好的选择。

随机数的产生和随机数的意义

在前面说过，随机数（RNs）是Monte Carlo的基石，随机试验是Monte Carlo的实现核心。最初，随机抽样试验的形式是投掷飞镖、掷硬币、投骰子等，并观察理论概率和频率。随着计算机技术的出现和发展，这种多次重复进行的随机抽样试验特别适合计算机来执行。而随机数的产生也经历了大量的研究，比如各种随机抽样数字表，但是这种随机数表中的随机数总是有限的，而有时候我们要做的抽样试验次数要远大于随机数表中给出的随机数个数。利用物理现象产生的随机数是完全随机的，比如RAND公司曾经以随机脉冲源为信息源，用电子旋转轮产生随机数表，但是物理方法的代价太过昂贵。通过数学方法，不断迭代，可以产生一系列伪随机数，但是显然伪随机数并不是真正随机的。不过，在迭代过程开始之前，每一项都是不可预测的，对产生的伪随机数，只要可以通过一系列局部随机性检验（比如均匀性、独立性检验等）就可以作为随机数使用，因此通过数学方法产生的随机数称为伪随机数。我们所熟知的matlab中的rand()函数就是用于产生0~1之间服从均匀分布的伪随机数。

随机数的意义

下面依然通过蒙特卡罗积分的例子来说明随机数的重要意义。

假设我们现在要估计一个积分值

首先，我们引入一个在(0,1)上服从均匀分布的随机变量x并构造一个新的随机变量 y=g(x) 。进一步，通过

我们把待求的未知量 I （积分值）表示成了一个随机变量y的期望值。假设随机变量x代表某个具体实验中的物理量，那么，我们便可以用期望的频率解释去估计 I ，通过多次重复试验，然后观察每次试验x的取值xi，并计算相应的yi值 yi=g(xi) ，然后再求平均值，以此来求出 I 的一个近似值。

我们将具有某些性质的数据xi称为随机数，这样，只要我们能够产生这些随机数，就可以确定出待求积分值 I 。这就是随机数的意义。

但是，实际上“什么是随机数？随机数如何产生？能否得到真正的随机数？”这一系列问题并没有一个标准的答案。因为随机数就像概率一样，在理论上和经验上的意义是完全不同的。在理论上，随机数是按照一个抽象模型定义的抽象概念抽象概念。在经验上，随机数是一个实数序列，它要么来自一个具体实验的物理数据，要么是一个确定的计算机程序的输出。

富兰克林（Franklin）对随机数的定义：如果一个序列具有一个均匀分布的随机变量独立样本构成的所有无穷序列所具有的所有性质，称这个序列是随机的。

莱莫（Lehmer）对随机数的定义：随机序列是一个模糊的概念，包含了这样一些思想，即序列的每一项对非初始项来说是不可预测的，同时其数字通过了一定数目的检验，这一概念与统计学家的习惯一致，并在某种程度上取决于这些序列的用途。

下面是从实际用途考虑给出的随机数的两种定义。

①概念上的定义：如果一个序列 xi 是在重复试验空间上独立同分布的随机变量序列 Xi 的样本，即 xi=Xi(ξ) ，则被称为是随机的。

这个定义和富兰克林的定义类似。区别在于，在富兰克林定义中，数列 xi 具有独立同分布随机变量所具有的所有性质，而在本定义中，xi 等于独立随机变量序列 Xi 的样本。

②经验上的定义：如果一个数列 xi 的统计性质与从随机试验中得到的随机数据是相同的，则被称为是随机的。

随机数的产生

在蒙特卡罗计算中，随机数主要通过计算机程序产生，当然，随机数也可以是随机试验产生的观测数据。比如，投掷一个骰子可以产生1-6的随机序列，放射性物质粒子发射之间的时间间隔生成的指数分布的观测样本。

目前计算机并不能直接产生任意分布的随机数序列，已有的算法只能生成服从 (0,m) 上均匀分布的随机整数序列，（但是任意分布的序列可以在此基础上间接得到）。

下面是一个最简单的产生随机数序列的算法

其中，f 是一个由最近的已得到的 r 个数确定的函数，z 是 f 除以 m 的余数。序列 z 的生成是非线性的，递归关系由整数 m ，函数 f ，以及初值 z1确定。这样设计的随机数发生器所生成的随机数的质量依赖于函数 f 的形式，f 越复杂，产生的随机数质量越高。

莱莫算法：一种最简单的递归产生随机数的算法

其中，m 是一个很大的素数（莱莫建议大小为2^31-1），a 是一个整数。

随机数的产生算法很多，我们不需要去深究，因为现在的编程语言大都已经封装好了产生随机数的函数，我们只需要去调用即可。并且也不建议自己去生成随机数序列（确切的说是伪随机数）。

Monte Carlo方法实例

投硬币模拟试验

硬币投掷试验是最经典的随机试验之一，下面通过计算机程序来模拟硬币投掷试验。随机投掷一枚均匀的硬币，正面朝上的次数应该是服从参数为 (1,p) 的二项分布，p 为正面朝上的概率。

matlab中产生二项分布伪随机数的函数为 r=binornd(n,p) ，参数与返回值如下：

n - Number of trials positive integer | array of positive integersp - Probability of success for each trial scalar value | array of scalar values

r - Random numbers from binomial distribution scalar value | array of scalar values

matlab试验代码如下

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % file:	硬币投掷模拟试验
% % author:	mindtechnist
% %
% % log:
% % 
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%end-of-header%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; clear; clc;
N = 2000; % 仿真试验次数
count = 0; % 结果为正面的次数
freq = zeros(1,N); % 频率
prob = 0.5*ones(1,N); % 概率
for i = 1:N
    count = count + binornd(1, 0.5); % 均匀硬币，正反面出现的概率都是0.5
    freq(i) = count/i;
end
figure
hold on;
plot(1:N, freq, '-r');
plot(1:N, prob, '-k');
legend("频率", "概率");

仿真结果如下

该试验证明了，频率稳定于概率。

古典概率试验

在一个袋子中有10个完全相同的球，分别标号为1~10。每次任取一球，记录号码后将球放回，然后继续取球。即有放回的抽取。现在有放回的抽取3个球，则3个球的号码均为偶数的概率是多少？

通过数学计算可以得到理论概率为0.125。下面做500次重复的仿真试验，并取其均值。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % file:   古典概率模拟试验
% % author: mindtechnist
% %
% % log:
% % 
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%end-of-header%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; clear; clc;

N = 500;  

M = 1000;

freq = zeros(1,N);

for k = 1:N
    freq(k) = subfun(M);
end

freq_mean = mean(freq);

figure
hold on; box on;
plot(freq, '-r');
line([0,N],[freq_mean,freq_mean],'LineWidth',2,'Color','k');
legend("估计概率", "均值");

function p = subfun(M)

    count = 0; % 全为偶数的次数
     
    for k=1:M
     
        ball1 = unidrnd(10); % r = unidrnd(n) 从由最大值 n 指定的离散均匀分布中生成随机数。
        ball2 = unidrnd(10);
        ball3 = unidrnd(10);
     
        if( (mod(ball1,2)==0) && (mod(ball2,2)==0) && (mod(ball3,2)==0) )
            count = count + 1;
        end
    end
     
    p = count/M; % 计算频率

end

仿真结果如下

蒙特卡罗方法求圆周率

如下图所示，假设圆的半径为1，正方形边长为2，均匀的产生随机变量点 (x,y) ，x和y的取值都在[-1,1]内，并以圆心为原点。共产生N个随机点，计算距离原点不大于1的点的个数（即落在圆内的点的个数）n，那么圆周率的值近似为4n/N。( πr²/4 = n/N )

python代码如下

import numpy as np

print('投点次数：')
n = eval(input())

x = np.random.uniform(-1, 1, n) # 均匀分布
y = np.random.uniform(-1, 1, n)

d = np.sqrt(x**2 + y**2)
res = sum(np.where(d <= 1.0, 1, 0))

pi = 4.0 * res / n

print("圆周率的估计值为：")
print(round(pi, 6))

运行结果

蒙特卡罗法计算定积分

蒙特卡罗方法在计算多重积分的问题中，有着良好的应用。因为在多重积分运算中，蒙特卡罗方法的误差与积分重数无关。尽管蒙特卡罗方法的精度较低，但是它能够很快的给出一个低精度的模拟结果。

随机投点法

假设待求的定积分表达式，以及函数图像如下所示

求函数 f(x) 的定积分实际上就是求图中阴影部分的面积。假设矩形面积为 S=M(b-a) ，待求阴影部分面积为 θ ，现在我们向矩形区域随机投点，那么落在阴影部分的概率可以通过下面的数学公式计算出来

这就是随机投点法求定积分的原理。

回到上面的定积分计算问题，我们可以通过数学方法求出这个定积分的解为7.2432.下面我们编写matlab程序，使用随机投点法求解该定积分。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % file:   随机投点法计算定积分
% % author: mindtechnist
% %
% % log:
% % 
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%end-of-header%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; clear; clc;

N = [1000, 10000, 100000, 1000000]; % 投点总数
value = zeros(1,length(N));

for i = 1:length(N)
    value(i) = simulation(N(i));
end

% figure
% hold on;
% plot(value);
value

function theta = simulation(N)

    % 积分上下限
    a = 0; b = 4; 
    
    M = 4; % f(x)的上界（矩形的上界）
    
    count = 0; % 落在阴影部分的点数
    
    for i = 1:N
      
        % 随机投点
        x = unifrnd(a,b);
        y = unifrnd(0,M);
     
        if (cos(x)+2.0) >= y % y <= f(x) 说明落入阴影区
            count = count+1; 
        end
    end
     
    freq = count/N;
     
    theta = freq*(b-a)*M;

end

通过仿真结果可以看到，随着仿真次数的增加，结果越来越接近数学上算出来精确解。

平均值法

假设要计算的积分形式为

matlab代码如下

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % file:   平均值法计算定积分
% % author: mindtechnist
% %
% % log:
% % 
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%end-of-header%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; clear; clc;

% 积分上下限
a = 0; b = 1;

N = 100000; 

value = zeros(1,10);

for k = 1:10
    value(k) = experiment(a,b,N);
end

value

function result = experiment(a, b, N)

    count = 0;
    xrandnum = unifrnd(a,b,1,N);

    for k = 1:N
        count = count + exp(xrandnum(1,k));
    end    

    result = count/N;

end

运行结果如下

计算定积分是蒙特卡罗方法的重要应用之一，特别是在复杂积分运算中，通过蒙特卡罗方法模拟是一个有效的解决手段。