数学建模十大算法01-蒙特卡洛算法（Monte Carlo）

蒙特卡洛方法也称为 计算机随机模拟方法，它源于世界著名的赌城——摩纳哥的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 引例

为了求得圆周率 π 值，在十九世纪后期，有很多人作了这样的试验：将长为 $l$ 的一根针任意投到地面上，用针与一组相间距离为 $a（l＜a）$的平行线相交的频率代替概率P，再利用准确的关系式：$p=\frac{2l}{\pi a}$ ，求出 π 值。（布丰投针）

注意：当针和平行线相交时有，针的中点x与针与直线的夹角φ满足 $x\le \frac{1}{2}sin\phi$

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次求平均pi
    m = 0;  
    x = rand(1, n) * a / 2 ;
    phi = rand(1, n) * pi;
    for i=1:n
        if x(i) <= l / 2 * sin(phi