数学建模十大算法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 引例
为了求得圆周率 π 值,在十九世纪后期,有很多人作了这样的试验:将长为 lll 的一根针任意投到地面上,用针与一组相间距离为 a(l<a)a( l<a)a(l<a)的平行线相交的频率代替概率P,再利用准确的关系式:p=2lπap=\frac{2l}{\pi a}p=πa2l ,求出 π 值。(布丰投针)
注意:当针和平行线相交时有,针的中点x与针与直线的夹角φ满足 x≤12sinφx≤\frac{1}{2}sinφx≤21sinφ
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)])
由于一次模拟的结果具有偶然性,因此我们可以重复100次后再来求一个平均的pi,这样子得到的结果会更接近真实的pi值。
result = zeros(100,1); % 初始化保存100次结果的矩阵
l = 0.520; a = 1.314;
n = 1000000;
for num = 1:100 % 重复100次求平均pim = 0; x = rand(1, n) * a / 2 ;phi = rand(1, n) * pi;for i=1:nif x(i) <= l / 2 * sin(phi (i))m = m + 1;endendp = m / n;mypi = (2 * l) / (a * p);result(num) = mypi; % 把求出来的myphi保存到结果矩阵中
end
mymeanpi = mean(result); % 计算result矩阵中保存的100次结果的均值
disp(['蒙特卡罗方法得到pi为:', num2str(mymeanpi)])
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 蒙特卡洛求解积分
θ=∫abf(x)dx\theta=\int_a^b f(x) d x θ=∫abf(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:Nif (x(i)^2+y(i)^2 <= 1)count = count + 1;end
end
PI = 4*count/N
正方形内部有一个相切的圆,它们的面积之比是π/4。现在,在这个正方形内部,随机产生1000000个点(即1000000个坐标对 (x, y)),计算它们与中心点的距离,从而判断是否落在圆的内部。如果这些点均匀分布,那么圆内的点应该占到所有点的 π/4,因此将这个比值乘以4,就是π的值。
【例】 计算定积分∫01x2dx\int_{0}^{1} x^{2} d x ∫01x2dx
计算函数 y =x2x^{2}x2在 [0, 1] 区间的积分,就是求出红色曲线下面的面积。这个函数在 (1,1) 点的取值为1,所以整个红色区域在一个面积为1的正方形里面。在该正方形内部,产生大量随机点,可以计算出有多少点落在红色区域(判断条件 y < x2x^{2}x2)。这个比重就是所要求的积分值。
N = 10000;
x = rand(N,1);
y = rand(N,1);
count = 0;
for i = 1:Nif (y(i) <= x(i)^2)count = count + 1;end
end
result = count/N
蒙特卡洛算法对于涉及不可解析函数 或 概率分布的模拟及计算是个有效的方法。
【例】 套圈圈问题。(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()
设投圈半径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)
注:上图中红圈为物品,散点图为模拟1000次投圈过程中,投圈中心点的位置散布。
然后,我们来计算1000次投圈过程中,投圈套住物品的占比情况。
print(len([xy for xy in points if xy[0] ** 2 + xy[1] ** 2 < (8-5) ** 2]) / N) # 物品半径为5cm,投圈半径为8cm,xy是一个坐标
输出结果为:0.015
代表投1000次,只有15次能够套住物品,这是小概率事件,知道我们为什么套不住了吧~
3.3 书店买书(0-1规划问题)
解:设 i=1,2...,6i=1,2...,6i=1,2...,6 表示ABCDEF六家商城, j=1,2...,5j=1,2...,5j=1,2...,5 表示B1、B2…B5五本书。记 mijm_{i j}mij 为第 jjj 本书在第 iii 家店的售价, qiq_{i}qi 表示第 iii 的运费。引入0-1变量 xijx_{i j}xij 如下:
那么,我们的目标函数就是总花费,包括两个方面:1)五本书总价格;2)运费。
书价=∑j=15[∑i=16(xij⋅mij)]书价 = \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(xij⋅mij)]
书店买书问题的蒙特卡罗的模拟代码实现:
%% 代码求解
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 5924 45 23 54 4422 45 23 53 5328 47 17 57 4724 42 24 47 5927 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); endif money < min_money % 判断刚刚随机生成的这组数据的花费是否小于最小花费,如果小于的话min_money = money % 我们更新最小的花费min_result = result % 用这组数据更新最小花费的结果end
end
循环执行的过程如下所示:
最终得到的最小花费为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:ntext(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_icoord_j = coord(j,:); x_j = coord_j(1); y_j = coord_j(2); % 城市j的横坐标为x_j,纵坐标为y_jd(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; % 初始化走过的路程为0path = randperm(n); % 生成一个1-n的随机打乱的序列for i = 1:n-1 result = d(path(i),path(i+1)) + result; % 按照这个序列不断的更新走过的路程这个值endresult = d(path(1),path(n)) + result; % 别忘了加上从最后一个城市返回到最开始那个城市的距离if result < min_result % 判断这次模拟走过的距离是否小于最短的距离,如果小于就更新最短距离和最短的路径min_path = path;min_result = resultend
end
在运行过程中,我们选择查看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] 数学建模——蒙特卡罗算法(Monte Carlo Method)
[2] 数学建模之蒙特卡洛算法
[3] 蒙特卡洛方法到底有什么用?
[4] 数学建模 | 蒙特卡洛模拟方法 | 详细案例和代码解析(清风课程) ★★推荐
数学建模十大算法01-蒙特卡洛算法(Monte Carlo)相关推荐
- 数学建模_数学模型的分类数学建模十大算法
数学模型的分类 数学建模十大算法 1.蒙特卡罗算法 (该算法又称随机性模拟算法, 是通过计算机仿真来解决问题的算法, 同时可以 通过模拟可以来检验自己模型的正确性,比较好用的算法) 2.数据拟合.参数 ...
- 数学建模十大算法02—插值与拟合(拉格朗日插值、三次样条插值、线性最小二乘法……)
文章目录 引入 一.插值 1.1 分段线性插值 1.2 牛顿插值法 1.3 拉格朗日插值多项式 1.4 样条插值 1.4.1 三次样条插值 1.5 二维插值 1.5.1 插值节点为网格节点 1.5.2 ...
- 数学建模十大算法(收藏)
1.蒙特卡罗算法(该算法又称随机性模拟算法,是通过计算机仿真来解决问题的算法,同时可以通过模拟可以来检验自己模型的正确性,是比赛时必用的方法) 2.数据拟合.参数估计.插值等数据处理算法(比赛中通常会 ...
- 数学建模十大算法04—图论算法(最短路径、最小生成树、最大流问题、二分图)
文章目录 一.最短路径问题 1.1 两个指定顶点之间的最短路径 1.1.1 Dijkstra算法 1.1.2 Matlab函数 1.2 每对顶点之间的最短路径 1.2.1 Dijkstra算法 1.2 ...
- java计算椭圆的面积_java算法3_蒙特卡洛方法(Monte Carlo method)求PI和椭圆面积
蒙特卡洛方法,是一种以概率统计理论为指导的一类非常重要的数值计算方法.是指使用随机数来解决很多计算问题的方法.蒙特卡洛方法的名字来源于摩纳哥的一个城市蒙特卡洛,该城市以×××业闻名,而蒙特卡洛方法正是 ...
- 数学建模-基本算法理解-蒙特卡洛算法
二,蒙特卡洛算法 1.蒙特卡洛算法的来源: 基本思想使用随机数来估算计算的问题. 2.蒙特卡洛算法的基本思想: 3.蒙特卡洛优缺点: 上述图片均来自于百度文库 (觉得讲的很好,这里就不再一一赘述.) ...
- DSt:数据结构的最强学习路线之数据结构知识讲解与刷题平台、刷题集合、问题为导向的十大类刷题算法(数组和字符串、栈和队列、二叉树、堆实现、图、哈希表、排序和搜索、动态规划/回溯法/递归/贪心/分治)总
DSt:数据结构的最强学习路线之数据结构知识讲解与刷题平台.刷题集合.问题为导向的十大类刷题算法(数组和字符串.栈和队列.二叉树.堆实现.图.哈希表.排序和搜索.动态规划/回溯法/递归/贪心/分治)总 ...
- 华为徐文伟:后香农时代,面向数学的十大挑战问题
本文为2020年8月28日徐文伟在长沙由中国工业与应用数学学会举办的"数学促进企业创新发展论坛"上的发言 来源:中国科学院院刊 徐文伟 华为技术有限公司董事,华为战略研究院院长 后 ...
- 2022年第三届MathorCup高校数学建模挑战赛——大数据竞赛(baseline)
教育部<高等学校人工智能创新行动计划>教技[2018]3号,鼓励对计算机专业类的智能科学与技术.数据科学与大数据技术等专业进行调整和整合,鼓励各个领域与大数据进行深度融合,通过大数据技术促 ...
- 2021年MathorCup高校数学建模挑战赛—大数据竞赛A题二手车估价问题解题思路
MathorCup高校数学建模挑战赛-大数据竞赛 A题 二手车估价问题 原题再现: 随着我国的机动车数量不断增长,人均保有量也随之增加,机动车以"二手车"形式在流通环节,包括二 ...
最新文章
- concat mysql sql注入_sql注入-mysql注入基础及常用注入语句
- golang 字符串 去首尾字符
- html4的语法,HTML——语法
- C#学习笔记——通用对话框
- android tabhost黑色背景,android更改FragmentTabHost背景和文本颜色
- Beego框架简介准备搭建分布式爬虫
- 抓鸡 抓服务器 1433 3306 全自动效率抓鸡
- Latex:字体设置
- 新中大连接服务器文件,新中大服务器数据库未能连接
- 对个人来说,最好的记账方法是什么?
- VS2003添加.BMP资源
- 蓝桥杯:合唱队形(C语言)
- 【防火墙流控配置 基于主机的带宽控制】
- kafka第二次课!!!
- python最长的单词判断_Python 找出英文单词列表(list)中最长单词链
- CSS---通向臃肿的道路(关于 “separation of concerns” (SoC)的原则)
- ATM(异步传输模式)是什么?
- 漏损分析与控制技术——漏损分析技术
- 专题7:动态规划 记忆化搜索
- 服务器电脑主板维修,电脑服务器主板维修