国赛学习——蒙特卡洛模拟

蒙特卡洛方法概述

1. 定义：蒙特卡洛方法又称统计模拟法，是一种随机模拟方法，以概率论和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。
2. 原理：由大数定理可知，当样本容量足够大时，时间的发生频率即为其概率

讨论

1. 蒙特卡洛模拟不是一种算法，是一种方法
2. 蒙特卡洛与计算机仿真的关系：计算机仿真（模拟）早起称为蒙特卡洛方法
3. 蒙特卡洛与枚举法：是枚举法的一种变异

三门问题

你参加⼀档电视节⽬，节⽬组提供了ABC三扇⻔，主持⼈告诉你，其中⼀扇⻔后边有辆汽⻋，其它两扇⻔后是空的。
假如你选择了B⻔，这时，主持⼈打开了C⻔，让你看到C⻔后什么都没有，然后问你要不要改选A⻔？
使用Matlab求解
randi([a,b],m,n)函数可以在指定区间[a,b]内随机取出大小为5*8的整数矩阵
randi([a,b])在区间内随机取出1个整数
字符串的连接方式：

1. strcat('abc','cef')
2. ['字符串1','字符串2']

将数值转换为字符串，num2str函数mystr=num2str(1124)
disp(num2str(1224),'lalala')输出函数
在成功条件下的概率（条件概率）

clear,clc
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	%如果xy相同，那么我们只有不改变主意时才能赢
		a = a + 1;	
		b = b + 0;	
	else 	%如果xy不同，那么我们改变主意时才能赢
		a = a + 0;
		b = b + 1;
	end
end
disp(['蒙特卡洛方法的到的不改变主意的获奖概率为：', num2str(a/n)]);
disp(['蒙特卡洛方法的到的改变主意的获奖概率为：', num2str(b/n)]);

考虑失败情况的代码（无条件概率）

clear,clc
n = 100000; %n代表蒙特卡洛模拟重复次数
a = 0;	%a表示不改变主意时能赢得汽车的次数
b = 0;	%b表示改变主意时能赢得汽车的次数
c = 0;	%表示没有获奖的次数
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	%如果xy相同，那么我们只有不改变主意时才能赢
		if(change == 0)	%不改变主意
			a = a + 1;	
		else		%改变主意
			c = c + 1;
		end
	else 	%如果xy不同，那么我们改变主意时才能赢
		if change == 0	%不改变主意
			c = c + 1;
		else			%改变主意
			b = b + 1;
		end
	end
end
disp(['蒙特卡洛方法的到的不改变主意的获奖概率为：', num2str(a/n)]);
disp(['蒙特卡洛方法的到的改变主意的获奖概率为：', num2str(b/n)]);

模拟排队问题

引入符号
$$\{\begin{array}{c}{c}_i:第i个客户到达时间\\ b_i:第i个客户开始服务时间\\ e_i:第i个客户服务结束时间\end{array}$$
题目中的信息：

1. 第i-1个客户和第i个客户到达的间隔时间xi～E(0,1) 参数入=0.1的指数分布（均值为10）
2. 第i个客户的服务持续时间yi～N(10,4)，若yi<1则令yi=1

$$\{\begin{array}{c}{c}_i=c_{i-1}+x_i（初始值c_0=0）\\ e_i=b_i+y_i（初始值e_0=0）\\ b_i=max(c_i,e_{i-1})（初始值b_1=c_1）\end{array}$$
第i个客户开始服务的时间取决于该客户的到达时间和上一个客户的服务结束时间
使用Matlab求解：
normrnd(MU,SIGMA):生成一个服从正态分布(MU参数代表均值,SIGMA参数代表标准差,方差开根号是标准差)的随机数
exprnd(M)表示生成一个均值为M的指数分布随机数(其对应的参数为1/M)
mean函数是用来求解均值的函数mean([1,2,3])
tic函数和toc函数可以用来返回代码运行的时间，例如我们要计算一段代码的运行时间，就可以在这段代码前加上tic，在这段代码后加上toc

tic
a = 2^100
toc

有约束的非线性规划问题

线性规划问题详见国赛学习——5种数学规划模型.
用蒙特卡洛模拟的思路：

1. 从约束条件中求出每个变量受限的大致范围（ai <= xi <= bi）
2. 在上述范围中用随机数生成若干组试验点，并验证它们是否满足所有的约束条件，若满足，则将其划分到可行组，再从可行组中找到函数的最大值或最小值

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

format long g
5/7
5895*514100

unifrnd(a,b,m,n)可以输出在[a,b]之间均匀分布的随机数组成的m行n列的矩阵。(等价于 a + rand(m,n)*(b-a)),e.g.unifrnd(0,5,4,3)

% max f(x) = x1x2x3
% s.t.
(1) -x1+2x2+2x3>=0
(2) x1+2x2+2x3<=72
(3) x2<=20 & x2>=10
(4) x1-x2 == 10

clc,clear;
tic %计算tic和toc中间部分的代码的运行时间
n=10000000;
x1=unifrnd(20,30,n,1); 
x2=x1 - 10;
x3=unifrnd(-10,16,n,1); 
fmax=-inf;
for i=1:n
    x = [x1(i), x2(i), x3(i)]; 
    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;  
        end
    end
end
disp(strcat('蒙特卡罗模拟得到的最大值为',num2str(fmax)))
disp('最大值处x1 x2 x3的取值为：')
disp(X)
toc %计算tic和toc中间部分的代码的运行时间