蒙特卡洛组合法 java 蒙特卡洛

最近疯狂学习数模，接触到了许多数模中常用的算法，今天来过一遍蒙特卡洛算法，算是复习一遍吧。

1.简介

蒙特·卡罗方法（Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。与它对应的是确定性算法。由于是以概率统计理论为基础的方法，故以赌城蒙特卡洛命名。

2.原理介绍

蒙特卡洛算法简单点说可以理解为通过模拟随机过程来估计某个随机事件发生的概率（用频率来估计概率，该算法的迭代次数较多）。实例演示的模式可能更易于理解，接下来就举几个例子来说明吧。

3.实例

例子一：求圆周率
现在设一个圆的半径为R，它的外切正方形的边长就为2*R,圆的面积为
$S_1=\Pi*R^2$,正方形的面积为$S_2=4*R^2$.将二者相除再化简，就可以得到圆周率$\Pi=\frac{4\ast S_1}{S_2}$现在只需要求出圆形和正方形的面积之比就可以求出圆周率了，可以采用随机投点法，在正方形区域内大量投点，然后统计圆形内的点数N₁，总投点数（即正方形内的点数）为N₂，当N₂足够大时可以近似认为:$\frac{S_1}{S_2}=\frac {N_1}{N_2}$现在只需要知道$\frac {N_1}{N_2}$即可，在此借助MATLAB编程来模拟投点过程，代码如下：

clear
%代码中以位于一象限的四分之一单位圆来计算
n=0.0;%用n来统计落入圆中的点
m=10000000;%设置m为投点次数，即为正方形中的点

x=0;y=0;
for i = 1:m
    x=rand();
    y=rand();%每次模拟投点都生成一个满足0<x<1,0<y<1的随机点
    if x*x+y*y<=1.0
        n=n+1;%点如果在四分之一圆内，n就加一
    end   
end
p=4*n/m;%根据公式计算圆周率估计值p
disp('落在四分之一圆内的点的个数为：n=');disp(n)
disp('通过蒙特卡洛模拟得到的圆周率为：p=');disp(p);

运行结果如下图所示：

可以看出p已经很接近理论值了，所以我们就顺利用蒙特卡洛方法求出圆周率啦。

例子二：利用蒙特卡洛算法求定积分
这也算是蒙特卡洛算法较常用的场景了，定积分从几何方面来理解就是被积函数曲线在积分上下限范围内与坐标轴围成图形的面积。设我们所求的定积分表达式为$\int^b_a{f(x)}{\rm d}x$$f(x)$为被积函数，$b,a$分别为积分上下限，蒙特卡洛算法求定积分的主要思想如下: 将区间$[a,b]$分为n个宽相同的小区间，n尽量大，每段区间依次取$x_1,x_2,...,x_n$来计算该区间对应的函数值（由于区间长度很小，可近似认为认为该区间上自变量对应的函数值都相等），函数值分别为$f(x_1),f(x_2),...,f(x_n)$，我们可以生成n个位于区间$[a,b]$内的随机数并求出其对应函数值的平均值$m$,再代入式子：$S=m*n$这就算出了函数与坐标轴围成图形的面积，这个面积就是我们所求的定积分的值。个人感觉蒙特卡洛求定积分就是将被积函数与坐标轴在积分上下限内所围成的不规则图形的面积近似转化成规则矩形的面积，该矩形的长和宽分别为$(b-a),m$。还有要注意不能直接生成随机函数值求平均值来计算!，这样做偏差很大（本人亲测）
下面是代码部分，以求解定积分$\int^1_0{e^x}{\rm d}x$为例子。代码如下：

clear;
%先直接计算定积分的值
syms x1
y1=int(exp(x1),x1,0,1);%matlab中求定积分的函数
disp('该积分的精确值为：');
vpa(y1,5)%将y1保留5位有效数字输出

%画出函数图像
% x1=0:0.1:1;
% y1=exp(x1);
% plot(y1,x1,'r-');

%再用蒙特卡洛模拟来计算该积分的值
n=10000;
x2=unifrnd(0,1,n,1);%生成维度为n*1的位于区间[0,1]的随机数
S=mean(exp(x2));%mean函数用于求平均值
disp('蒙特卡洛计算的定积分值为：');disp(S);

运行结果如下图：

可以看出，计算结果也很接近，这样，我们就算是成功用蒙特卡洛的思想来计算定积分啦。

4.总结

蒙特卡洛算法是数模中特别经典的算法，它的思想很简单，应用也十分广泛，数模必备哦。当然它也有不足之处，那就是它本身是有高计算需求的，即使对于相当简单的问题往往也需要海量的计算，因此，我们要学会在合适的场景下使用该算法。