【数模修炼之旅】03 蒙特卡洛算法 深度解析（教程+代码）

【数模修炼之旅】03 蒙特卡洛算法 深度解析（教程+代码）

接下来 C君将会用至少30个小节来为大家深度解析数模领域常用的算法，大家可以关注这个专栏，持续学习哦，对于大家的能力提高会有极大的帮助。

1 蒙特卡洛算法介绍及应用

‌蒙特卡洛方法的基本思想是通过构造符合一定规则的随机数来解决数学上的各种问题。当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，蒙特卡洛方法通过某种"实验"的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。‌

这么看来蒙特卡洛方法的理论支撑其实是概率论或统计学中的大数定律。基本原理简单描述是先大量模拟，然后计算一个事件发生的次数，再通过这个发生次数除以总模拟次数，得到想要的结果。下面我们以三个经典的小实验来学习下蒙特卡洛算法思想。

1.计算圆周率pi（π）值

实验原理：在正方形内部有一个相切的圆，圆面积/正方形面积之比是(PixRxR)/(2Rx2R)= Pi/4。在这个正方形内随机产生n个点，假设点落在圆内的概率为P，那么P=圆面积/正方形面积，则P= Pi/4。如何计算点落在圆内的概率P？可以计算点与中心点的距离，判断是否落在圆的内部，若这些点均匀分布，用M表示落到圆内投点数 ， N表示总的投点数，则圆周率Pi=4P=4xM/N。

实验步骤：

（1）将圆心设在原点（0，0），以R为半径形成圆，则圆面积为PixRxR

（2）将该圆外接正方形, 坐标为（-R，-R）（R，-R）（R， R）（-R，R），则该正方形面积为R*R

（3）随即取点（X，Y），使得-R <=X<=R并且-R <=Y<=R，即点在正方形内

（4）通过公式 XxX+YxY<= RxR判断点是否在圆周内（直角三角形边长公式）。

（5）设所有点（也就是实验次数）的个数为N，落在圆内的点（满足步骤4的点）的个数为M，则P=M/N，于是Pi=4xM/N。

（6）运行结果为3.143052

def cal_pai_mc(n=1000000):

r = 1.0

a, b = (0.0, 0.0)

x_neg, x_pos = a - r, a + r

y_neg, y_pos = b - r, b + r

m = 0

for i in range(0, n+1):

x = random.uniform(x_neg, x_pos)

y = random.uniform(y_neg, y_pos)

if x**2 + y**2 <= 1.0:

m += 1

return (m / float(n)) * 4

python代码：

import random

# 生成 N 个随机点
N = 1000000
inside_circle = 0

for _ in range(N):
    x = random.uniform(0, 1)
    y = random.uniform(0, 1)
    if x**2 + y**2 <= 1:
        inside_circle += 1

# 估计圆周率
pi_estimate = 4 * inside_circle / N
print(f"估计的圆周率：{pi_estimate}") 

2.计算函数定积分值

实验原理：若要求函数f(x)从a到b的定积分，我们可以用一个比较容易算得面积的矩型包围在函数的积分区间上（假设其面积为Area），定积分值其实就是求曲线下方的面积。随机地向这个矩形框里面投点，统计落在函数f(x)下方的点数量占所有点数量的比例为P，那么就可以据此估算出函数f(x)从a到b的定积分为Area×P。此处我们将a和b设为0和1，函数f(x)=x2。

运行结果为0.333749

def cal_integral_mc(n = 1000000):

x_min, x_max = 0.0, 1.0

y_min, y_max = 0.0, 1.0

m = 0

for i in range(0, n+1):

x = random.uniform(x_min, x_max)

y = random.uniform(y_min, y_max)

# x*x > y 表示该点位于曲线的下面。

if x*x > y:

m += 1

#所求的积分值即为曲线下方的面积与正方形面积的比

return m / float(n)

3.计算函数极值，可避免陷入局部极值

实验原理：极值是“极大值” 和 “极小值”的统称。如果一个函数在某点的一个邻域内处处都有确定的值，函数在该点的值大于或等于在该点附近任何其他点的函数值，则称函数在该点的值为函数的“极大值”。如果函数在该点的值小于或等于在该点附近任何其他点的函数值，则称函数在该点 的值为函数的“极小值”。此处在区间[-2,2]上随机生成一个数，求出其对应的y，找出其中最大值认为是函数在[-2,2]上的极大值。

运行结果发现极大值185.1204262706596, 极大值点为1.5144491499169481

def cal_extremum_mc(n = 1000000):

y_max = 0.0

x_min, x_max = -2.0, 2.0

y = lambda x:200*np.sin(x)*np.exp(-0.05*x)#匿名函数

for i in range(0, n+1):

x0 = random.uniform(x_min, x_max)

if y(x0) > y_max:

y_max = y(x0)

x_max = x0

return y_max, x_max

2 蒙特卡洛算法的基本步骤

1. 定义问题：确定需要估计的目标函数。例如计算积分、求解最优化问题、或估计某个概率。
2. 随机采样：生成大量的随机样本。样本可以来自目标分布或通过某种方式转换得来。
3. 计算样本值：对于每一个随机样本，计算目标函数的值。这些值将用于估计目标的期望或概率。
4. 求平均值：将所有样本的计算结果进行平均，作为最终估计的值。
5. 收敛性与精度：增加样本数量，提高估计的精度，并检验算法是否收敛。

3 蒙特卡洛算法代码（matlab+python）

均以估算圆周率为例

3.1 matlab

function pi_estimation = estimate_pi(num_samples)
    inside_circle = 0;

    for i = 1:num_samples
        x = rand;
        y = rand;
        if x^2 + y^2 <= 1
            inside_circle = inside_circle + 1;
        end
    end

    pi_estimation = 4 * inside_circle / num_samples;
end

% 设定样本数量
num_samples = 1000000;
pi_estimation = estimate_pi(num_samples);
fprintf('估计的圆周率：%f\n', pi_estimation);

3.2 python

import random

def estimate_pi(num_samples):
    inside_circle = 0

    for _ in range(num_samples):
        x = random.uniform(0, 1)
        y = random.uniform(0, 1)
        if x**2 + y**2 <= 1:
            inside_circle += 1

    pi_estimate = 4 * inside_circle / num_samples
    return pi_estimate

# 设定样本数量
num_samples = 1000000
pi_estimation = estimate_pi(num_samples)
print(f"估计的圆周率：{pi_estimation}")

需要参加数模竞赛的同学，可以看我的这个名片，会有最新的助攻哦：（大型比赛前会对名片进行更新）