蒙特卡洛模拟(详细,上手,有例题)

布丰投针实验

假设有一组间距为d的平行直线，把一根长度为 l（ l<=d）的针随机地投掷在该平面上，求针与直线相交的概率

概率公式：，

丰投针问题的求解思路蕴含了蒙特卡罗方法的基本思想，即通过大量随机试验，利用随机事件发生的频率来近似求解数学问题。蒙特卡罗方法在现代科学计算、数值模拟、金融工程等众多领域有着广泛应用。

预备的matlab知识

Rand(m,n)生成m*n的0-1均匀分布的随机矩阵

如果只有一个参数默认是方阵

Unifrnd(-2,2,3,2)//生成3*2的矩阵，数值在-2到2之间任意取

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 * si n(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)])

蒙特卡洛方法的概述

概述:统计模拟方法，是一种以概率统计理论为指导的数值计算方法，通过大量随机样本，去模拟和估计复杂系统或问题的行为与结果。

基本原理

随机抽样：依据问题所涉及的概率分布，在计算机上生成大量符合该分布的随机数。例如，在估计不规则图形面积时，可在包含该图形的规则矩形区域内随机生成点。

模拟实验：使用这些随机数模拟实际问题中的各种随机因素，执行大量模拟实验。比如在模拟粒子在介质中的运动时，通过随机数确定粒子每次碰撞的方向和距离。

统计分析：对模拟结果进行统计分析，计算出感兴趣的统计量，如均值、方差等，以此估计问题的解。例如多次模拟投资组合收益后，计算平均收益来评估该投资策略的预期回报。

应用;

FIRST 三门模拟

你参加一档电视节目，节目组提供了ABC三扇门，主持人告诉你，其中一扇门后边有辆汽车，其它两扇门后是空的。

假如你选择了B门，这时，主持人打开了C门，让你看到C门后什么都没有，然后问你要不要改选A门?

预备matlab的知识

randi([1,5],3,2) 生成一个3*2的矩阵,他的值随机均有分布在1到5的整数

字符串连接:[s1,s2]==strcat(s1,s2)

num2str()//将数字转换为字符串

n = 100000; % n代表蒙特卡罗模拟重复次数

a = 0; % a表示不改变主意时能赢得汽车的次数

b = 0; % b表示改变主意时能赢得汽车的次数

c = 0; % c表示没有获奖的次数

for i= 1 : n % 开始模拟n次

x = randi([1,3]); % 随机生成一个1-3之间的整数x表示汽车出现在第x扇门后

y = randi([1,3]); % 随机生成一个1-3之间的整数y表示自己选的门

change = randi([0, 1]); % change =0 不改变主意，change = 1 改变主意

% 下面分为两种情况讨论：x=y和x~=y

if x == y % 如果x和y相同，那么我们只有不改变主意时才能赢

if change == 0 % 不改变主意

a = a + 1;

else % 改变了主意

c= c+1;

end

else % x ~= y ，如果x和y不同，那么我们只有改变主意时才能赢

if change == 0 % 不改变主意

c = c + 1;

else % 改变了主意

b= b + 1;

end

end

end

disp(['蒙特卡罗方法得到的不改变主意时的获奖概率为：', num2str(a/n)]);

disp(['蒙特卡罗方法得到的改变主意时的获奖概率为：', num2str(b/n)]);

disp(['蒙特卡罗方法得到的没有获奖的概率为：', num2str(c/n)]);

最终1/2的概率失败,成功的概率里面1/3的概率属于不改变主意,2/3改变主意

SECOND模拟排队问题

假设某银行工作时间只有一个服务窗口，工作人员只能逐个的接待顾客。当来的顾客较多时，一部分顾客就需要排队等待。

假设:

1)顾客到来的间隔时间服从参数为0.1的指数分布

2)每个顾客的服务时间服从均值为10，方差为4的正态分布(单位为分钟，若服务时间小于1分钟，则按1分钟计算)

3)排队按先到先服务的规则，且不限制队伍的长度，每天工作时长为8小时。

试回答下面的问题:

1)模拟一个工作日，在这个工作日共接待了多少客户，客户平均等待的时间为多少?

2)模拟100个工作日，计算出平均每日接待客户的个数以及每日客户的平均等待时长。

Matlab预备知识

normrnd(0,1)生成一个服从正态分布的随机数(均值为0,标准差为1,方差为1^2)

exprnd(5) 生成一个均值为5的质数分布随机数(对应的参数是1/5)

mean(求平均值的函数)

tic和toc函数,是一组计算运行时长的函数

clear

tic %计算tic和toc中间部分的代码的运行时间

i = 1; % i表示第i个客户，最开始取i=1

w = 0; % w用来表示所有客户等待的总时间，初始化为0

e0 = 0; c0 = 0; % 初始化e0和c0为0

x(1) = exprnd(10); % 第0个客户(假想的)和第1个客户到达的时间间隔

c(1) = c0 + x(1); % 第1个客户到达的时间

b(1) = c(1); % 第1个客户的开始服务的时间

while b(i) <= 480 % 开始设置循环，只要第i个顾客开始服务的时间(时刻)小于480，就可以对其服务（银行每天工作8小时，折换为分钟就是480分钟）

y(i) = normrnd(10,2); % 第i个客户的服务持续时间，服从均值为10方差为4(标准差为2)的正态分布

if y(i) < 1 % 根据题目的意思：若服务持续时间不足一分钟，则按照一分钟计算

y(i) = 1;

end

e(i) = b(i) + y(i); % 第i个客户结束服务的时间 = 第i个客户开始服务的时间 + 第i个客户的服务持续时间

wait(i) = b(i) - c(i); % 第i个客户等待的时间 = 第i个客户开始服务的时间 - 第i个客户到达银行的时间

w = w + wait(i); % 更新所有客户等待的总时间

i = i + 1; % 增加一名新的客户

x(i) = exprnd(10); % 这位新客户和上一个客户到达的时间间隔

c(i) = c(i-1) + x(i); % 这位新客户到达银行的时间 = 上一个客户到达银行的时间 + 这位新客户和上一个客户到达的时间间隔

b(i) = max(c(i),e(i-1)); % 这个新客户开始服务的时间取决于其到达时间和上一个客户结束服务的时间

end

n = i-1; % n表示银行一天8小时一共服务的客户人数

t = w/n; % 客户的平均等待时间

disp(['银行一天8小时一共服务的客户人数为: ',num2str(n)])

disp(['客户的平均等待时间为: ',num2str(t)])

toc %计算tic和toc中间部分的代码的运行时间

第二题就是while外面套一个day的循环,就省略了

THIRD:解决有约束的非线性规划问题

从约束条件中找出每个变量受限的大致范围,然后在这个范围内生成若干的点,并且验证是否满足约束,从满足约束的组种找到函数的最大值和最小值,这样只要次数足够多,就可以在一定程度上解决这道题

max f(x) = x1*x2*x3

s.t.

(1) -x1+2*x2+2*x3>=0

(2) x1+2*x2+2*x3<=72

(3) x2<=20 & x2>=10

(4) x1-x2 == 10

Matlab预备知识:

format long g

可以得到一个精确的解,而不是一个科学计数法的解

clc,clear;

tic %计算tic和toc中间部分的代码的运行时间

n=10000000; %生成的随机数组数

x1=unifrnd(20,30,n,1); % 生成在[20,30]之间均匀分布的随机数组成的n行1列的向量构成x1

x2=x1 - 10;

x3=unifrnd(-10,16,n,1); % 生成在[-10,16]之间均匀分布的随机数组成的n行1列的向量构成x3

fmax=-inf; % 初始化函数f的最大值为负无穷（后续只要找到一个比它大的我们就对其更新）

for i=1:n

x = [x1(i), x2(i), x3(i)]; %构造x向量, 这里千万别写成了：x =[x1, x2, x3]

if (-x(1)+2*x(2)+2*x(3)>=0) & (x(1)+2*x(2)+2*x(3)<=72) % 判断是否满足条件

result = x(1)*x(2)*x(3); % 如果满足条件就计算函数值

if result > fmax % 如果这个函数值大于我们之前计算出来的最大值

fmax = result; % 那么就更新这个函数值为新的最大值

X = x; % 并且将此时的x1 x2 x3保存到一个变量中

end

end

end

disp(strcat('蒙特卡罗模拟得到的最大值为',num2str(fmax)))

disp('最大值处x1 x2 x3的取值为：')

disp(X)

toc %计算tic和toc中间部分的代码的运行时间

FOURTH:解决0-1规划问题

某同学要从六家线上商城选购五本书籍 B1,B2,B3,B4,B5，每本书籍在不同商家的售价以及每个商家的单次运费如下表所示，请给该同学制定最省钱的选购方案。(注:在同一个店买多本书也只会收取一次运费)

Matlab预备知识:

unique函数:删除一个矩阵.or 向量的重复值,并且按照从小到大排序

代码:

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

min_money % 18+39+48+17+47+20

min_result

FIFTH:导弹追踪

位于坐标原点的A船向位于其正东方20个单位的B船发射导弹，导弹始终对准B船，B船以时速V单位(常数)沿东北方向逃逸。若导弹的速度为3V，导弹的射程是50个单位，画出导弹运行的曲线，导弹是否能在射程内击中B船?

matlab预备知识:

mod(8,3)  求8/3的余数

绘图:

x=1:0.01:3

y=x.^2

plot(x,y)//绘图

axis([0,3,0,10])//设置展示的范围

matlab代码;

clear;clc

v=200; % 任意给定B船的速度（后期我们可以再改的）

dt=0.0000001; % 定义时间间隔

x=[0,20]; % 定义导弹和B船的横坐标分别为x(1)和x(2)

y=[0,0]; % 定义导弹和B船的纵坐标分别为y(1)和y(2)

t=0; % 初始化导弹击落B船的时间

d=0; % 初始化导弹飞行的距离

m=sqrt(2)/2; % 将sqrt(2)/2定义为一个常量，使后面看起来很简洁

dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 导弹与B船的距离

while(dd>=0.001) % 只要两者的距离足够大，就一直循环下去。（两者距离足够小时表示导弹击中，这里的临界值要结合dt来取，否则可能导致错过交界处的情况）

t=t+dt; % 更新导弹击落B船的时间

d=d+3*v*dt; % 更新导弹飞行的距离

x(2)=20+t*v*m; y(2)=t*v*m; % 计算新的B船的位置 （注：m=sqrt(2)/2）

dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 更新导弹与B船的距离

tan_alpha=(y(2)-y(1))/(x(2)-x(1)); % 计算斜率，即tan(α)

cos_alpha=sqrt(1/(1+tan_alpha^2)); % sec(α)^2 = (1+tan(α)^2)

sin_alpha=sqrt(1-cos_alpha^2); % sin(α)^2 +cos(α)^2 = 1

x(1)=x(1)+3*v*dt*cos_alpha; y(1)=y(1)+3*v*dt*sin_alpha; % 计算新的导弹的位置

if d>50 % 导弹的有效射程为50个单位

disp('导弹没有击中B船');

break; % 退出循环

end

if d<=50 & dd<0.001 % 导弹飞行的距离小于50个单位且导弹和B船的距离小于0.001（表示击中）

disp(['导弹飞行',num2str(d),'单位后击中B船'])

disp(['导弹飞行的时间为',num2str(t*60),'分钟'])

end

end

SIXTH:旅行商问题

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

完全图:完全图是一个简单的无向图，其中每对不同的顶点之间都恰连有一条边相连

matlab的预备知识:

plot([1,2],[5,10],'-o') % 画出一条线段，x范围是[1, 2] ，y范围是[5,10]

text(1.5,7.5,'xxx') % 在坐标(1.5,7.5)处标上文本：xxx

close

randperm(5) % 生成1-5组成的一个随机序列(类似于洗牌的操作)

matlab代码:

clear;clc

% 只有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 = 10000000; % 蒙特卡罗模拟的次数

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

min_path

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