这部分主要得会写代码

引例

蒲丰投针问题：

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

1） 取一张白纸，在上面画上许多条间距为a的平行线。
2） 取一根长度为l（l≤a） 的针，随机地向画有平行直线的纸上掷n次，观察针与直线相交的次数，记为m。
3）计算针与直线相交的概率。

布丰本人证明了，这个概率是：P=2lπaP=\tfrac{2l}{πa}(其中π为圆周率）

证明：当x≤l2sinϕ,0≤x≤a2,0≤ϕ≤πx\leq \tfrac{l}{2} sin\phi,0\le x\le \tfrac{a}{2},0\le \phi \le \pi时，才符合题意，其中x为直线的中点到最近的平行线的距离 ϕ\phi为直线与最近平行线的夹角，由几何概率的定义可以求得P，由于大数定理，当n足够大的时候，频率趋近于概率，以此来得到π\pi的值。

用matlab模拟：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
% rand(m,n)可以生成一个m行n列的矩阵，每个元素值在[0,1]之间
% rand(n)只有一个值就是n*n的方阵
% a+rand(m,n)*(b-a) 可以生成一个m行n列的矩阵，每个元素值在[a,b]之间 等价于unifrnd(a,b,m,n)

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)])

可以重复多次求平均值

概述

定义

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

原理

大数定理 当样本容量足够大时，事件的发生频率即为其概率。

蒙特卡罗模拟不是一个通用的代码，而是一种思想。

问题

问题一：三门问题

问题的提出：你参加一档节目，节目有ABC三扇门，其中一扇门后面有汽车，另外两扇门后是空的。假设你选择了B门，这时C门被打开了，C门后面什么都没有，那么要不要改选A门？
结论：如果换了A得奖的概率是三分之二，不换门概率是三分之一。（好神奇啊）
蒙特卡罗模拟：
在赢的条件下，没考虑输的，条件概率：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
% randi([a,b],m,n)函数可在指定区间[a,b]内随机取出大小为m*n的整数矩阵
n = 100000;  % n代表蒙特卡罗模拟重复次数
a = 0;  % a表示不改变主意时能赢得汽车的次数
b = 0;  % b表示改变主意时能赢得汽车的次数
for i= 1 : n  % 开始模拟n次
    x = randi([1,3]);  % 随机生成一个1-3之间的整数x表示汽车出现在第x扇门后
    y = randi([1,3]);  % 随机生成一个1-3之间的整数y表示自己选的门
    % 下面分为两种情况讨论：x=y和x~=y
    if x == y   % 如果x和y相同，那么我们只有不改变主意时才能赢
        a = a + 1;     b = b + 0;
    else  % x ~= y ，如果x和y不同，那么我们只有改变主意时才能赢
        a = a + 0;     b = b +1;
    end
end
disp(['蒙特卡罗方法得到的不改变主意时的获奖概率为：', num2str(a/n)]);
disp(['蒙特卡罗方法得到的改变主意时的获奖概率为：', num2str(b/n)]);

有输的情况，无条件概率：

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
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. 顾客到来的间隔时间服从参数为0.1的指数分布
  2. 每个顾客的服务时间服从均值为10，⽅差为4的正态分布(单位为分钟，若服务时间⼩于1分钟，则按1分钟计算)
  3. 排队按先到先服务的规则，且不限制队伍的⻓度，每天⼯作时⻓为8⼩时。

试回答下⾯的问题：

  1. 模拟⼀个⼯作⽇，在这个⼯作⽇共接待了多少客户，客户平均等待的时间为多少?
  2. 模拟100个⼯作⽇，计算出平均每⽇接待客户的个数以及每⽇客户的平均等待时⻓。

建立模型：
x(i)表示第i-1个客户和第i个客户到达的间隔时间，服从参数为0.1的指数分布
y(i)表示第i个客户的服务持续时间，服从均值为10方差为4(标准差为2)的正态分布 (若小于1则按1计算)
c(i)表示第i个客户的到达时间，那么c(i) = c(i-1) + x(i)，初始值c0=0
b(i)表示第i个客户开始服务的时间
e(i)表示第i个客户结束服务的时间，初始值e0=0
第i个客户结束服务的时间 = 第i个客户开始服务的时间 + 第i个客户的服务持续时间
即：e(i) = b(i) + y(i）
第i个客户开始服务的时间取决于该客户的到达时间和上一个客户结束服务的时间
即：b(i) = max(c(i),e(i-1))，初始值b1=c1;
第i个客户等待的时间 = 第i个客户开始服务的时间 - 第i个客户到达银行的时间
即：wait(i) = b(i) - c(i)
w表示所有客户等待时间的总和
假设一天内银行最终服务了n个顾客，那么客户的平均等待时间t = w/n

预备知识

1
2
3
4
5
6
7
8
9
10
% normrnd(MU,SIGMA):生成一个服从正态分布(MU参数代表均值,SIGMA参数代表标准差,方差开根号是标准差)的随机数
normrnd(10,2)  % 均值为10 标准差为2（方差为4）的正态分布随机数
% exprnd(M)表示生成一个均值为M的指数分布随机数(其对应的参数为1/M)
exprnd(5)  % 均值为5的指数分布随机数（对应的参数为0.2）
% mean函数是用来求解均值的函数（第一期视频第五讲）
mean([1,2,3])
% tic函数和toc函数可以用来返回代码运行的时间，例如我们要计算一段代码的运行时间，就可以在这段代码前加上tic，在这段代码后加上toc (我的微信公众号"数学建模学习交流"中有一篇推送《为什么要对代码初始化》中使用过这对函数)
tic
a = 2^100
toc

问题1的代码

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

问题2的代码

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
tic  %计算tic和toc中间部分的代码的运行时间
day = 100;  % 假设模拟100天
n = zeros(day,1); % 初始化用来保存每日接待客户数结果的矩阵
t = zeros(day,1); % 初始化用来保存每日客户平均等待时长的矩阵
for k = 1:day
    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(k) = i-1; % n(k)表示银行第k天服务的客户人数
    t(k) = w/n(k); % t(k)表示该银行第k天客户的平均等待时间
end
disp([num2str(day),'个工作日中，银行每日平均服务的客户人数为: ',num2str(mean(n))])
disp([num2str(day),'个工作日中，银行每日客户的平均等待时间为: ',num2str(mean(t))])
toc  %计算tic和toc中间部分的代码的运行时间

matlab代码汇总

还有一些经典的问题这里不再赘述，题目很多，蒙特卡洛模拟主要是用其随机生成的思想来模拟问题的情况，用频率估计概率。遇到具体问题还要具体分析，这里汇总一些常用的matlab语句。

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
rand(m,n) %可以生成一个m行n列的矩阵，每个元素值在[0,1]之间

rand(n) %只有一个值就是n*n的方阵

a+rand(m,n)*(b-a) %可以生成一个m行n列的矩阵，每个元素值在[a,b]之间 等价于
unifrnd(a,b,m,n)

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

normrnd(MU,SIGMA) %生成一个服从正态分布(MU参数代表均值,SIGMA参数代表标准差,方差开根号是标准差)的随机数

exprnd(M) %表示生成一个均值为M的指数分布随机数(其对应的参数为1/M)

mean[a,b,c] %用来求解均值的函数

tic% tic函数和toc函数可以用来返回代码运行的时间
a = 2^100
toc

format long g  %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）

unique([5 6 8 8 4 1 6 2 2 4 8 4 5 6]) %剔除一个矩阵或者向量的重复值，并将结果按照从小到大的顺序排列 

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

find %假设a是一个向量，那么find(a)可以用来返回这个向量中非零元素的下标，如果a中所有元素都为0，则返回空值
find(a)  % 找到a中所有非0元素的位置
find(a == 5)  % 找到a中等于5的元素的位置
find(a == 5,1)  % 找到a中第一个等于5的元素的位置
find(a == min(a))   % 找到a中最小元素的位置

isempty(A) %函数可以用来判断A是否为空, 如果A为空， isempty(A) 返回逻辑值1(true)，否则返回逻辑值0(false)。

randsrc(m,n,[alphabet; prob]) % 以一定的概率产生随机数
% m和n表示生成的随机数矩阵的行数和列数
% alphabet表示需要产生的随机数的数字，用一个行向量表示
% prob表示这些数字出现的概率大小，用一个行向量表示，向量长度和alphabet向量要完全相同, 且这些概率的和要为1
% 比如：要产生1、4、 6这三个数。它们分别出现的概率为 0.1、0.2、0.7，如何设计程序使得按照这个概率产生10个随机数呢？
alphabet = [1 4 6]; prob = [0.1 0.2 0.7];
randsrc(10,1,[alphabet; prob])