蒙特卡罗(洛)模拟

Posted pxlsdz

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了蒙特卡罗(洛)模拟相关的知识,希望对你有一定的参考价值。


title: 蒙特卡罗(洛)模拟
date: 2020-02-27 21:26:53
categories: 数学建模
tags: [ MATLAB, 模拟]
mathjax: true


引例

布丰投针实验

技术图片

法国数学家布丰(1707-1788)最早设计了投针试验。

这一方法的步骤是:

1) 取一张白纸,在上面画上许多条间距为a的平行线。

2) 取一根长度为l(l≤a) 的针,随机地向画有平行直线的纸上掷n次,观察针与直线相交的次数,记为m。

3)计算针与直线相交的概率.

18世纪,法国数学家布丰提出的“投针问题”,记载于布丰1777年出版的著作中:“在平面上画有一组间距为a的平行线,将一根长度为l(l≤a)的针任意掷在这个平面上,求此针与平行线中任一条相交的概率。”

布丰本人证明了,这个概率是:
[ p=frac{2l}{pi a} ]

(其中π为圆周率

由于它与π有关,于是人们想到利用投针试验来估计圆周率的值。

布丰惊奇地发现:有利的扔出与不利的扔出两者次数的比,是一个包含π的表示式.如果针的长度等于a/2,那么有利扔出的概率为1/π.扔的次数越多,由此能求出越为精确的π的值。

下面是利用这个公式,用概率的方法得到圆周率的近似值的一些资料。

试验者 时间 投掷次数 相交次数 圆周率估计值
Wolf 1850年 5000 2532 3.1596
Smith 1855年 3204 1218.5 3.1554
C.De Morgan 1860年 600 382.5 3.137
Fox 1884年 1030 489 3.1595
Lazzerini 1901年 3408 1808 3.1415929
Reina 1925年 2520 859 3.1795

证明

技术图片

代码

%%  蒙特卡罗用于布丰投针实验

%% (1)预备知识
%  rand(m,n)函数产生由在[0,1]之间均匀分布的随机数组成的m行n列的矩阵(或称为数组)。
rand(5,4)
%     0.8300    0.1048    0.2396    0.4398
%     0.5663    0.1196    0.8559    0.5817
%     0.9281    0.2574    0.3013    0.9355
%     0.3910    0.3173    0.2108    0.1676
%     0.3645    0.4372    0.8819    0.9232
rand(3) % 若只给一个输入,则会生成一个方阵
%     0.1709    0.4951    0.0541
%     0.9374    0.8500    0.6155
%     0.2400    0.3156    0.5741
% a + rand(m,n)*(b-a) 可以输出在[a,b]之间均匀分布的随机数组成的m行n列的矩阵。
-2 + rand(3,2) * (2 - (-2))  % 输出在[-2,2]之间均匀分布的随机数组成的3行2列的矩阵。
%    -1.2743    0.6013
%    -1.3084    0.0766
%     1.5075    0.7563
% a + rand(m,n)*(b-a)等价于unifrnd(a,b,m,n)
unifrnd(-2,2,3,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 * 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)])


%% (3) 由于一次模拟的结果具有偶然性,因此我们可以重复100次后再来求一个平均的pi
result = zeros(100,1);  % 初始化保存100次结果的矩阵
l =  0.520;     a = 1.314;
n = 1000000;    
for num = 1:100
    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)])

蒙特卡罗方法概述

定义

? 蒙特卡罗?法?称统计模拟法,是?种随机模拟?法,以概率和统计理论?法为基础的?种计算?法,是使?随机数 (或更常?的伪随机数)来解决很多计算问题的?法。将所求解的问题同?定的概率模型相联系,?电?计算机实现统计模拟或抽样,以获得问题的近似解。为象征性地表明这??法的概率统计特征,故借?赌城蒙特卡罗命名。

提出

? 蒙特卡罗?法于20世纪40年代美国在第?次世界?战中研制原?弹的“曼哈顿计划”计划的成员S.M.乌拉姆和J. 冯·诺伊曼?先提出。数学家冯·诺伊曼?驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种?法,为它蒙上了? 层神秘?彩。在这之前,蒙特卡罗?法就已经存在。1777年,法国Buffon提出?投针实验的?法求圆周率,这被认 为是蒙特卡罗?法的起源。

原理

由?数定理可知,当样本容量?够?时,事件的发?频率即为其概率。

讨论

  1. 蒙特卡罗是?种算法吗 ?

    算法(Algorithm)是指解题?案的准确?完整的描述,是?系列解决问题的清晰指令。蒙特卡罗准确的来说只是? 种思想,或者是是?种?法。如果我们所求解的问题与概率模型有?定的关联,那么我们就可以使?计算机多次模拟事 件发?,以获得问题的近似解。从数学建模?度来看,?家千万别认为蒙特卡罗有?个通?的代码。每个问题对应的代 码都是不同的,我们分析清楚题?后,就要??进?编写适?于这个题?的代码。

  2. 蒙特卡罗与计算机仿真的关系

    计算机仿真(模拟)早期称为蒙特卡罗?法,是??利?随机数实验求解随机问题的?法,其主要应?在复杂问题的数 值模拟上。但随着计算机的性能的提?以及各种新兴产业的发展,?前计算机仿真涉及的内容要?得多,例如过程控制、动 画仿真、图像静态模拟等都属于计算机仿真的应?(例如?计算机模拟原?弹爆炸的过程、使?计算机?成特效??等)。 在数学建模中,我们不?刻意的去区分两者的区别,?家只需要知道如何使?计算机对问题进?模拟即可。

  3. 蒙特卡罗与枚举法

    枚举法是我们中学就接触的算法,就是把所有可能发?情况都考虑进去,最终计算出来?个确定结果。这就与蒙特卡罗 ?法的想法很类似,蒙特卡罗法模拟的次数越多,计算的就越准确。由于?活中有许多事件发?的结果都有?限种可能(例如?个连续分布的取值),因此我们不可能枚举出所有的结果,这时候就只能通过蒙特卡罗模拟,将?个不确定性的问题转 化成很多个确定性问题,并得到?个近似解,因此蒙特卡罗算法也可以看成是枚举法的?种变异。

已经学过的例?

第?期 视频 (AOL) 第七讲 : 多元回归分析

探究内?性对回归结果的影响

第?期 视频 (AOL) 番外篇 : 基于熵权法对印的模型的修正

信息熵和标准差的关系

蒙特卡罗?法 的应?实例

三?问题

? 你参加?档电视节?,节?组提供了ABC三扇?, 主持?告诉你,其中?扇?后边有辆汽?,其它两扇? 后是空的。 假如你选择了B?,这时,主持?打开了C?,让你 看到C?后什么都没有,然后问你要不要改选A??

有约束的?线性规划问题

技术图片

技术图片

技术图片

代码

%%  蒙特卡罗求解有约束的非线性规划问题
% 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

%% (1)预备知识
%  (1) format long g  可以将Matlab的计算结果显示为一般的长数字格式(默认会保留四位小数,或使用科学计数法)
5/7
5895*514100
format long g
5/7
5895*514100
%  (2)unifrnd(a,b,m,n)可以输出在[a,b]之间均匀分布的随机数组成的m行n列的矩阵。(等价于 a + rand(m,n)*(b-a))
unifrnd(0,5,4,3)
%           4.07361843196589          3.16179623112705          4.78753417717149
%            4.5289596853781         0.487702024997048          4.82444267599638
%           0.63493408146753          1.39249109433524         0.788065408387741
%            4.5668792806951          2.73440759602492          4.85296390880308

%% (2)代码部分
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中间部分的代码的运行时间


%% (3)缩小范围重新模拟得到更加精确的取值
clc,clear;
tic %计算tic和toc中间部分的代码的运行时间
n=10000000; %生成的随机数组数
x1=unifrnd(22,23,n,1);  % 生成在[22,23]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=x1 - 10;
x3=unifrnd(11,13,n,1);  % 生成在[11,13]之间均匀分布的随机数组成的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中间部分的代码的运行时间

书店 买书 问题 ( 0 - 1 规划)(无脑)

技术图片

技术图片

%% 书店买书问题的蒙特卡罗的模拟
%% (1)预备知识
% (1)unique函数: 剔除一个矩阵或者向量的重复值,并将结果按照从小到大的顺序排列  
% adj.  唯一的; 独一无二的   [ju'ni:k]
unique([1 2 5; 6 8 9;2 4 6])   
unique([5 6 8 8 4 1 6 2 2 4 8 4 5 6])

% (2)randi([a,b],m,n)函数可在指定区间[a,b]内随机取出大小为m*n的整数矩阵
randi([-5,5],2,6)

%% (2)代码求解
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 = k([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

导弹追踪问题

技术图片

技术图片

技术图片

技术图片

技术图片

%%  蒙特卡罗用于模拟导弹追击问题
clear;clc
%% (1)预备知识
% mod(m,n)表示求m/n的余数
mod(8,3)
mod(1000,50)

% 设置横纵坐标的范围并标上字符
x = 1:0.01:3;
y = x .^ 2;
plot(x,y)  % 画出x和y的图形
axis([0 3 0 10])  % 设置横坐标范围为[0, 3] 纵坐标范围为[0, 10]
pause(3)  % 暂停3秒后再继续接下来的命令
text(2,4,'清风')  % 在坐标为(2,4)的点上标上字符串:清风
close % 关闭图形窗口



%% (2) 代码求解
% 1. 不画追击的示意图
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

% 2. 画追击的示意图
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船的距离
for i=1:2
    plot(x(i),y(i),'.k','MarkerSize',1);  % 画出导弹和B船所在的坐标,点的大小为1,颜色为黑色(k),用小点表示
    grid on;  % 打开网格线
    hold on;  % 不关闭图形,继续画图
end
axis([0 30 0 10])  % 固定x轴的范围为0-30  固定y轴的范围为0-10
k = 0;  % 引入一个变量  为了控制画图的速度(因为Matlab中画图的速度超级慢)
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)  计算出cos(α)
    sin_alpha=sqrt(1-cos_alpha^2);  % 利用公式: sin(α)^2 +cos(α)^2 = 1  计算出sin(α)
    x(1)=x(1)+3*v*dt*cos_alpha;   y(1)=y(1)+3*v*dt*sin_alpha;   % 计算新的导弹的位置
    k = k +1 ;  
    if mod(k,500) == 0   % 每刷新500次时间就画出下一个导弹和B船所在的坐标  mod(m,n)表示求m/n的余数
        for i=1:2
            plot(x(i),y(i),'.k','MarkerSize',1);
            hold on; % 不关闭图形,继续画图
        end
        pause(0.001);  % 暂停0.001s后再继续下面的操作
    end
    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

作业

作业1

技术图片

%%  蒙特卡罗求解非线性规划问题
% min f(x) =2*(x1^2)+x2^2-x1*x2-8*x1-3*x2
% s.t.
% (1) 3*x1+x2>9
% (2) x1+2*x2<16
% (3) x1>0 & x2>0

%% (1)初次寻找最小值的代码
clc,clear;
format long g   %可以将Matlab的计算结果显示为一般的长数字格式(默认会保留四位小数,或使用科学计数法)
tic %计算tic和toc中间部分的代码的运行时间
n=10000000; %生成的随机数组数
x1=unifrnd(0,16,n,1);  % 生成在[0,16]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=unifrnd(0,8,n,1);  % 生成在[0,8]之间均匀分布的随机数组成的n行1列的向量构成x2
fmin=+inf; % 初始化函数f的最小值为正无穷(后续只要找到一个比它小的我们就对其更新)
for i=1:n
    x = [x1(i), x2(i)];  %构造x向量, 这里千万别写成了:x =[x1, x2]
    if (3*x(1)+x(2)>9)  &  (x(1)+2*x(2)<16)     % 判断是否满足条件
        result = 2*(x(1)^2)+x(2)^2-x(1)*x(2)-8*x(1)-3*x(2);  % 如果满足条件就计算函数值
        if  result  < fmin  % 如果这个函数值小于我们之前计算出来的最小值
            fmin = result;  % 那么就更新这个函数值为新的最小值
            X = x;  % 并且将此时的x1 x2 保存到相应的变量中
        end
    end
end
disp(strcat('蒙特卡罗模拟得到的最小值为',num2str(fmin)))
disp('最小值处x1 x2的取值为:')
disp(X)
toc %计算tic和toc中间部分的代码的运行时间


%% (2)缩小范围重新模拟得到更加精确的取值
clc,clear;
tic %计算tic和toc中间部分的代码的运行时间
n=10000000; %生成的随机数组数
x1=unifrnd(2,3,n,1);  % 生成在[2,3]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=unifrnd(2,3,n,1);  % 生成在[2,3]之间均匀分布的随机数组成的n行1列的向量构成x2
fmin=+inf; % 初始化函数f的最小值为正无穷(后续只要找到一个比它小的我们就对其更新)
for i=1:n
    x = [x1(i), x2(i)];  %构造x向量, 这里千万别写成了:x =[x1, x2]
    if (3*x(1)+x(2)>9)  &  (x(1)+2*x(2)<16)     % 判断是否满足条件
        result = 2*(x(1)^2)+x(2)^2-x(1)*x(2)-8*x(1)-3*x(2);  % 如果满足条件就计算函数值
        if  result  < fmin  % 如果这个函数值小于我们之前计算出来的最小值
            fmin = result;  % 那么就更新这个函数值为新的最小值
            X = x;  % 并且将此时的x1 x2 保存到相应的变量中
        end
    end
end
disp(strcat('蒙特卡罗模拟得到的最小值为',num2str(fmin)))
disp('最小值处x1 x2的取值为:')
disp(X)
toc %计算tic和toc中间部分的代码的运行时间

作业2

? 某设备上安装有四只型号规格完全相同的电子管,已知电子管寿命(假定为整数)为1000h~2000h之间的均匀分布.当电子管损坏时有两种维修方案,一是每次更换损坏的那一只;二是当其中一只损坏时四只同时更换.已知更换时间为换一只时需1h,4只同时换为2h.更换时机器因停止运转每小时的损失为20元,又每只电子管价格10元,试用模拟方法决定哪一个方案经济合理?

%% 选择决策方案的模拟
% 某设备上安装有四只型号规格完全相同的电子管,已知电子管寿命为1000--2000小时之间的均匀分布(假定为整数)。
% 当电子管损坏时有两种维修方案,一是每次更换损坏的那一只;二是当其中一只损坏时四只同时更换。
% 已知更换时间为换一只时需1小时,4只同时换为2小时。
% 更换时机器因停止运转每小时的损失为20元,又每只电子管价格10元,
% 试用模拟方法决定哪一个方案经济合理?

%% (1)预备知识
% randi([a,b],m,n)  随机生成m*n的矩阵,矩阵中的每个元素都是[a,b]中的随机整数
randi([1, 5],3,2)
randi([1, 5])  % 不写m*n代表只生成1个随机数

% find函数的用法
% find函数的用法在第一期视频:层次分析法那一节讲过,我们当时找最大特征值的位置
a = [2 3 5 1 7 5];
find(a)  % 找到a中所有非0元素的位置
find(a == 5)  % 找到a中等于5的元素的位置
find(a == 5,1)  % 找到a中第一个等于5的元素的位置
find(a == min(a))   % 找到a中最小元素的位置

%% (2)代码部分
clear;clc
T = 100000000;   % T表示模拟的总时间(单位为小时)
t = 0;   % 初始化当前时刻为0小时
c1 = 0; c2 = 0;  % 初始化两种方案的总花费都为0

%%  方案一
life = randi([1000,2000],1,4);  % 随机生成四个电子管的寿命,假设为整数
while t < T  % 只要现在的时刻没有超过总时刻,就不断循环下去
    result = min(life);  % 找出寿命最短的那一个电子管的寿命
    t = t+result+1;  % 现在的时间更改到有电子管损坏的时刻(加上1表示更换电子管需要花费的时间)
    c1 = c1 + 20 * 1 +10;  % 更新方案一的花费 
    k = find(life == result,1);   % 找到哪一个电子管是坏的
    life = life - result -1; % 更新所有电子管的寿命    
    life(k) = randi([1000,2000]);  % 把坏掉的那个电子管的寿命重置
end

%%  方案二
t = 0;   % 初始化当前时刻为0小时
while t < T  % 只要现在的时刻没有超过总时刻,就不断循环下去
    life = randi([1000,2000],1,4); % 随机生成四个电子管的寿命,假设为整数
    result = min(life); % 找出寿命最小的那一个电子管的寿命
    t = t+result+2;  % 现在的时间更改到有电子管损坏的时刻(加上2表示更换所有电子管需要花费的时间)
    c2 =c2 + 20 * 2 +40;  % 更新方案二的花费 
end

%% 两种方案的花费
c1
c2


以上是关于蒙特卡罗(洛)模拟的主要内容,如果未能解决你的问题,请参考以下文章

matlab随机模拟求面积

数据分析之蒙特卡洛模拟

数学模型:4. 蒙特卡罗模拟

蒙特卡罗(Monte Carlo) 模拟

R语言实战案例-蒙特卡罗方法(附实现代码)

R语言实战案例-蒙特卡罗方法(附实现代码)