查看: 3510|回复: 0

[其他资源] 数学建模十大算法01-蒙特卡洛算法（Monte Carlo）

[复制链接]

字体大小: 正常放大

杨利霞

5273 主题	82 听众	17万积分

TA的每日心情

	开心 2021-8-11 17:59

签到天数: 17 天

[LV.4]偶尔看看III

网络挑战赛参赛者

自我介绍: 本人女，毕业于内蒙古科技大学，担任文职专业，毕业专业英语。

群组: 2018美赛大象算法课程

群组: 2018美赛护航培训课程

群组: 2019年数学中国站长建

群组: 2019年数据分析师课程

群组: 2018年大象老师国赛优

电梯直达

1^#

发表于 2022-9-12 18:20 |只看该作者 |倒序浏览

|招呼Ta 关注Ta

数学建模十大算法01-蒙特卡洛算法（Monte Carlo）
文章目录
一、生成随机数
1.1 rand
1.2 unifrnd
1.3 联系与区别
二、引入
2.1 引例
2.2 基本思想
2.3 优缺点
三、实例
3.1 蒙特卡洛求解积分
3.2 简单的实例
3.3 书店买书（0-1规划问题）
3.4 旅行商问题（TSP）
参考文献

蒙特卡洛方法也称为计算机随机模拟方法，它源于世界著名的赌城——摩纳哥的Monte Carlo(蒙特卡洛)。它是基于对大量事件的统计结果来实现一些确定性问题的计算。使用蒙特卡洛方法必须使用计算机生成相关分布的随机数，Matlab给出了生成各种随机数的命令，常用的有 rand函数和 unifrnd。
一、生成随机数
1.1 rand
rand函数可用于产生由(0,1)之间均匀分布的随机数或矩阵。
Y = rand(n) 返回一个n×n的随机矩阵。
Y = rand(m,n) 或 Y = rand([m n]) 返回一个m×n的随机矩阵。

Y = rand(m,n,p,...) 或 Y = rand([m n p...]) 产生随机数组。

Y = rand(size(A)) 返回一个和A有相同尺寸的随机矩阵。

1.2 unifrnd
unifrnd 生成一组（连续）均匀分布的随机数。
R = unifrnd(A,B) 生成被A和B指定上下端点[A,B]的连续均匀分布的随机数组R。
如果A和B是数组，R(i,j)是生成的被A和B对应元素指定连续均匀分布的随机数。

R = unifrnd(A,B,m,n,...) 或 R = unifrnd(A,B,[m,n,...])
如果A和B是标量，R中所有元素是相同分布产生的随机数。
如果A或B是数组，则必须是mn…数组。

1.3 联系与区别
相同点：

二者都是利用rand函数进行随机值计算。
二者都是均匀分布。
【例】在区间[5,10]上生成400个均匀分布的随机数。

不同点：

unifrnd是统计工具箱中的函数，是对rand的包装。不是Matlab自带函数无法使用JIT加速。
rand函数可以指定随机数的数据类型。
二、引入
2.1 引例
为了求得圆周率 π 值，在十九世纪后期，有很多人作了这样的试验：将长为 l ll 的一根针任意投到地面上，用针与一组相间距离为 a （ l ＜ a ） a（ l＜a）a（l＜a）的平行线相交的频率代替概率P，再利用准确的关系式：p = 2 l π a p=\frac{2l}{\pi a}p=
πa
2l

  ，求出 π 值。（布丰投针）

注意：当针和平行线相交时有，针的中点x与针与直线的夹角φ满足 x ≤ 1 2 s i n φ x≤\frac{1}{2}sinφx≤
2
1

sinφ

l =  0.520;    % 针的长度（任意给的）
a = 1.314; % 平行线的宽度(大于针的长度l即可)
n = 1000000; % 做n次投针试验，n越大求出来的pi越准确
m = 0; % 记录针与平行线相交的次数
x = rand(1, n) * a / 2 ; % 在[0, a/2]内服从均匀分布随机产生n个数， x中每一个元素表示针的中点和最近的一条平行线的距离
phi = rand(1, n) * pi; % 在[0, pi]内服从均匀分布随机产生n个数，phi中的每一个元素表示针和最近的一条平行线的夹角
% axis([0,pi, 0,a/2]); box on;  % 画一个坐标轴的框架，x轴位于0-pi，y轴位于0-a/2，并打开图形的边框
for i=1:n  % 开始循环，依次看每根针是否和直线相交
if x(i) <= l / 2 * sin(phi (i))    % 如果针和平行线相交
      m = m + 1; % 那么m就要加1
%       plot(phi(i), x(i), 'r.') % 模仿书上的那个图，横坐标为phi，纵坐标为x , 用红色的小点进行标记
%       hold on  % 在原来的图形上继续绘制
end
end
p = m / n; % 针和平行线相交出现的频率
mypi = (2 * l) / (a * p);  % 我们根据公式计算得到的pi
disp(['蒙特卡罗方法得到pi为：', num2str(mypi)])

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

由于一次模拟的结果具有偶然性，因此我们可以重复100次后再来求一个平均的pi，这样子得到的结果会更接近真实的pi值。

result = zeros(100,1);  % 初始化保存100次结果的矩阵
l =  0.520;    a = 1.314;
n = 1000000;
for num = 1:100  % 重复100次求平均pi
m = 0;
x = rand(1, n) * a / 2 ;
phi = rand(1, n) * pi;
for i=1:n
      if x(i) <= l / 2 * sin(phi (i))
         m = m + 1;
      end
end
p = m / n;
mypi = (2 * l) / (a * p);
result(num) = mypi;  % 把求出来的myphi保存到结果矩阵中
end
mymeanpi = mean(result);  % 计算result矩阵中保存的100次结果的均值
disp(['蒙特卡罗方法得到pi为：', num2str(mymeanpi)])

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
2.2 基本思想
当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率，数学期望有关的量时，通过某种试验的方法，得出该事件发生的概率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。
当随机变量的取值仅为1或0时，它的数学期望就是某个事件的概率。或者说，某种事件的概率也是随机变量（仅取值为1或0）的数学期望。
2.3 优缺点
优点：（可以求解复杂图形的积分、定积分，多维数据也可以很快收敛）
1、能够比较逼真地描述具有随机性质的事物的特点及物理实验过程
2、受几何条件限制小
3、收敛速度与问题的维数无关
4、具有同时计算多个方案与多个未知量的能力
5、误差容易确定
6、程序结构简单，易于实现

缺点：
1、收敛速度慢
2、误差具有概率性
3、在粒子输运问题中，计算结果与系统大小有关

主要应用范围：

1、粒子输运问题（实验物理，反应堆物理）
2、统计物理
3、典型数学问题
4、真空技术
5、激光技术
6、医学
7、生物
8、探矿
……

注：所以在使用蒙特卡罗方法时，要“扬长避短”，只对问题中难以用解析（或数值）方法处理的部分，使用蒙特卡罗方法计算，对那些能用解析（或数值）方法处理的部分，应当尽量使用解析方法。

蒙特卡洛算法，采样越多，越接近最优解。（尽量找好的，但是不保证是最好的），它与拉斯维加斯算法的对比可参考：蒙特卡罗算法与拉斯维加斯算法。

三、实例
3.1 蒙特卡洛求解积分
θ = ∫ a b f ( x ) d x \theta=\int_a^b f(x) d x
θ=∫
a
b

f(x)dx

步骤如下：

在区间[a,b]上利用计算机均匀的产生n个随机数（可用matlab的unifrnd实现）
计算每一个随机数相对应的函数值f(x1)、f(x2)、f(xn)
计算被积函数值的平均值
3.2 简单的实例
【例】求π的值。

N = 1000000; % 随机点的数目
x = rand(N,1);  % rand 生成均匀分布的伪随机数。分布在（0~1）之间
y = rand(N,1);  % 矩阵的维数为N×1
count = 0;
for i = 1:N
if (x(i)^2+y(i)^2 <= 1)
   count = count + 1;
end
end
PI = 4*count/N
1
2
3
4
5
6
7
8
9
10
正方形内部有一个相切的圆，它们的面积之比是π/4。现在，在这个正方形内部，随机产生1000000个点（即1000000个坐标对 (x, y)），计算它们与中心点的距离，从而判断是否落在圆的内部。如果这些点均匀分布，那么圆内的点应该占到所有点的 π/4，因此将这个比值乘以4，就是π的值。

【例】计算定积分
∫ 0 1 x 2 d x \int_{0}^{1} x^{2} d x
∫
0
1

x
2
dx

计算函数 y =x 2 x^{2}x
2
在 [0, 1] 区间的积分，就是求出红色曲线下面的面积。这个函数在 (1,1) 点的取值为1，所以整个红色区域在一个面积为1的正方形里面。在该正方形内部，产生大量随机点，可以计算出有多少点落在红色区域（判断条件 y < x 2 x^{2}x
2
）。这个比重就是所要求的积分值。

N = 10000;
x = rand(N,1);
y = rand(N,1);
count = 0;
for i = 1:N
if (y(i) <= x(i)^2)
   count = count + 1;
end
end
result = count/N
1
2
3
4
5
6
7
8
9
10

蒙特卡洛算法对于涉及不可解析函数或概率分布的模拟及计算是个有效的方法。

【例】套圈圈问题。（Python代码）

在这里，我们设物品中心点坐标为(0,0)，物品半径为5cm。

import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np
import sys
circle_target = mpatches.Circle([0, 0], radius=5, edgecolor='r', fill=False)
plt.xlim(-80, 80)
plt.ylim(-80, 80)
plt.axes().add_patch(circle_target)  # 在坐标轴里面添加圆
plt.show()
1
2
3
4
5
6
7
8
9

设投圈半径8cm，投圈中心点围绕物品中心点呈二维正态分布，均值μ=0cm，标准差σ=20cm，模拟1000次投圈过程。

N = 1000  # 1000次投圈
u, sigma = 0, 20  # 投圈中心点围绕物品中心呈二维正态分布，均值为0，标准差为20cm
points = sigma * np.random.randn(N, 2) + u
plt.scatter([x[0] for x in points], [x[1] for x in points], c=np.random.rand(N), alpha=0.2)
1
2
3
4

注：上图中红圈为物品，散点图为模拟1000次投圈过程中，投圈中心点的位置散布。

然后，我们来计算1000次投圈过程中，投圈套住物品的占比情况。

print(len([xy for xy in points if xy[0] ** 2 + xy[1] ** 2 < (8-5) ** 2]) / N)  # 物品半径为5cm，投圈半径为8cm，xy是一个坐标
1
输出结果为：0.015
代表投1000次，只有15次能够套住物品，这是小概率事件，知道我们为什么套不住了吧~

3.3 书店买书（0-1规划问题）

解：设 i = 1 , 2... , 6 i=1,2...,6i=1,2...,6 表示ABCDEF六家商城， j = 1 , 2... , 5 j=1,2...,5j=1,2...,5 表示B1、B2…B5五本书。记 m i j m_{i j}m
ij

  为第 j jj 本书在第 i ii 家店的售价， q i q_{i}q
i

  表示第 i ii 的运费。引入0-1变量 x i j x_{i j}x
ij

  如下：

那么，我们的目标函数就是总花费，包括两个方面：1）五本书总价格；2）运费。

书价 = ∑ j = 1 5 [ ∑ i = 1 6 ( x i j ⋅ m i j ) ] 书价 = \sum_{j=1}^{5}\left[\sum_{i=1}^{6}\left(x_{i j} \cdot m_{i j}\right)\right]
书价=
j=1
∑
5

[
i=1
∑
6

(x
ij

⋅m
ij

)]

书店买书问题的蒙特卡罗的模拟代码实现：

%% 代码求解
min_money = +Inf;  % 初始化最小的花费为无穷大，后续只要找到比它小的就更新
min_result = randi([1, 6],1,5);  % 初始化五本书都在哪一家书店购买，后续我们不断对其更新
%若min_result = [5 3 6 2 3]，则解释为：第1本书在第5家店买，第2本书在第3家店买，第3本书在第6家店买，第4本书在第2家店买，第5本书在第3家店买
n = 100000;  % 蒙特卡罗模拟的次数
M = [18    39 29 48 59
      24 45 23 54 44
      22 45 23 53 53
      28 47 17 57 47
      24 42 24 47 59
      27 48 20 55 53];  % m_ij  第j本书在第i家店的售价
freight = [10 15 15 10 10 15];  % 第i家店的运费
for k = 1:n  % 开始循环
result = randi([1, 6],1,5); % 在1-6这些整数中随机抽取一个1*5的向量，表示这五本书分别在哪家书店购买
index = unique(result);  % 在哪些商店购买了商品，因为我们等下要计算运费
money = sum(freight(index)); % 计算买书花费的运费
% 计算总花费：刚刚计算出来的运费 + 五本书的售价
for i = 1:5
      money = money + M(result(i),i);
end
if money < min_money  % 判断刚刚随机生成的这组数据的花费是否小于最小花费，如果小于的话
      min_money = money  % 我们更新最小的花费
      min_result = result % 用这组数据更新最小花费的结果
end
end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
循环执行的过程如下所示：

最终得到的最小花费为189元，方案为第一本书在第1家店买，第二本书在第1家店买，第三本书在第4家店买，第四本书在第1家店买，第五本书在第4家店买。

3.4 旅行商问题（TSP）
一个售货员必须访问n个城市，这n个城市是一个完全图，售货员需要恰好访问所有城市一次，并且回到最终的城市。城市与城市之间有一个旅行费用，售货员希望旅行费用之和最少。

如图所示的完全图旅行费用最小时的路径为：城市1→城市3→城市2→城市4→城市1

案例代码实现：

% 只有10个城市的简单情况
coord =[0.6683 0.6195 0.4 0.2439 0.1707 0.2293 0.5171 0.8732 0.6878 0.8488 ;
            0.2536 0.2634 0.4439 0.1463 0.2293 0.761  0.9414 0.6536 0.5219 0.3609]' ;  % 城市坐标矩阵，n行2列
% 38个城市，TSP数据集网站(http://www.tsp.gatech.edu/world/djtour.html) 上公测的最优结果6656。
% coord = [11003.611100,42102.500000;11108.611100,42373.888900;11133.333300,42885.833300;11155.833300,42712.500000;11183.333300,42933.333300;11297.500000,42853.333300;11310.277800,42929.444400;11416.666700,42983.333300;11423.888900,43000.277800;11438.333300,42057.222200;11461.111100,43252.777800;11485.555600,43187.222200;11503.055600,42855.277800;11511.388900,42106.388900;11522.222200,42841.944400;11569.444400,43136.666700;11583.333300,43150.000000;11595.000000,43148.055600;11600.000000,43150.000000;11690.555600,42686.666700;11715.833300,41836.111100;11751.111100,42814.444400;11770.277800,42651.944400;11785.277800,42884.444400;11822.777800,42673.611100;11846.944400,42660.555600;11963.055600,43290.555600;11973.055600,43026.111100;12058.333300,42195.555600;12149.444400,42477.500000;12286.944400,43355.555600;12300.000000,42433.333300;12355.833300,43156.388900;12363.333300,43189.166700;12372.777800,42711.388900;12386.666700,43334.722200;12421.666700,42895.555600;12645.000000,42973.333300];

n = size(coord,1);  % 城市的数目

figure(1)  % 新建一个编号为1的图形窗口
plot(coord(:,1),coord(:,2),'o'); % 画出城市的分布散点图
for i = 1:n
text(coord(i,1)+0.01,coord(i,2)+0.01,num2str(i)) % 在图上标上城市的编号（加上0.01表示把文字的标记往右上方偏移一点）
end
hold on % 等一下要接着在这个图形上画图的

d = zeros(n); % 初始化两个城市的距离矩阵全为0
for i = 2:n
for j = 1:i
      coord_i = coord(i,

; x_i = coord_i(1); y_i = coord_i(2); % 城市i的横坐标为x_i，纵坐标为y_i
coord_j = coord(j,

; x_j = coord_j(1);    y_j = coord_j(2);  % 城市j的横坐标为x_j，纵坐标为y_j
      d(i,j) = sqrt((x_i-x_j)^2 + (y_i-y_j)^2); % 计算城市i和j的距离
end
end
d = d+d'; % 生成距离矩阵的对称的一面

min_result = +inf;  % 假设最短的距离为min_result，初始化为无穷大，后面只要找到比它小的就对其更新
min_path = [1:n]; % 初始化最短的路径就是1-2-3-...-n
N = 10000;  % 蒙特卡罗模拟的次数，清风老师设的次数为10000000，这里我为了快速得到结果，改为10000
for i = 1:N  % 开始循环
result = 0;  % 初始化走过的路程为0
path = randperm(n);  % 生成一个1-n的随机打乱的序列
for i = 1:n-1
      result = d(path(i),path(i+1)) + result;  % 按照这个序列不断的更新走过的路程这个值
end
result = d(path(1),path(n)) + result;  % 别忘了加上从最后一个城市返回到最开始那个城市的距离
if result < min_result  % 判断这次模拟走过的距离是否小于最短的距离，如果小于就更新最短距离和最短的路径
      min_path = path;
      min_result = result
end
end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
在运行过程中，我们选择查看min_result的变化：

最终得到的路径（不一定是最优的路径）为：

图中显示最短路径：

min_path = [min_path,min_path(1)]; % 在最短路径的最后面加上一个元素，即第一个点（我们要生成一个封闭的图形）
n = n+1;  % 城市的个数加一个（紧随着上一步）
for i = 1:n-1
   j = i+1;
coord_i = coord(min_path(i),

; x_i = coord_i(1); y_i = coord_i(2);
coord_j = coord(min_path(j),

; x_j = coord_j(1); y_j = coord_j(2);
plot([x_i,x_j],[y_i,y_j],'-') % 每两个点就作出一条线段，直到所有的城市都走完
pause(0.5) % 暂停0.5s再画下一条线段
hold on
end
1
2
3
4
5
6
7
8
9
10

参考文献
[1] 数学建模——蒙特卡罗算法（Monte Carlo Method）
[2] 数学建模之蒙特卡洛算法
[3] 蒙特卡洛方法到底有什么用？
[4] 数学建模 | 蒙特卡洛模拟方法 | 详细案例和代码解析（清风课程） ★★推荐
————————————————
版权声明：本文为CSDN博主「美式咖啡不加糖x」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/Serendipity_zyx/article/details/126592916

zan