蒙特卡洛算法代码实现

 01 算法概述 蒙特卡罗算法是以概率和统计理论方法为基础的一种计算方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。又称统计模拟法、随机抽样技术。由S.M.乌拉姆和J.冯·诺伊曼在20世纪40年代为研制核武器而首先提出。
它的基本思想是，为了求解数学、物理、工程技术以及管理等方面的问题 ，首先建立一个概率模型或随机过程，使它们的参数，如概率分布或数学期望等问题的解；然后通过对模型或过程的观察或抽样试验来计算所求参数的统计特征，并用算术平均值作为所求解的近似值。对于随机性问题，有时还可以根据实际物理背景的概率法则，用电子计算机直接进行抽样试验，从而求得问题的解答。简单来讲即采样越多越接近最优解。

02 算法特点 直接追踪粒子，物理思路清晰，易于理解。

采用随机抽样的方法，较真切的模拟粒子输运的过程，反映了统计涨落的规律。
不受系统多维、多因素等复杂性的限制，是解决复杂系统粒子输运问题的好方法。
MC程序结构清晰简单。
研究人员采用MC方法编写程序来解决粒子输运问题，比较容易得到自己想得到的任意中间结果，应用灵活性强。

%画函数图像
x=0:0.25:10;
y1=3*x;
y2=8-x;
plot(x,y1,x,y2)
 
axis([0 10 0 10]);
legend('y1=3x','y2=8-x');
title('蒙特卡洛算法');
text(2,6,'交点');
grid on
%产生随机点
x=unifrnd(0,8,[1,10000000]);
y=unifrnd(0,6,[1,10000000]);
%统计所在所求图形中的点
frequency=sum(y<3*x & x<=2)+sum(y<8-x & x>2);
%计算面积
area=6*8*frequency/10^7;

这段 MATLAB 代码用于演示蒙特卡洛算法估计平面图形的面积。以下是对每一句代码的解释：

1. matlabCopy code

   x=0:0.25:10; y1=3*x; y2=8-x; plot(x,y1,x,y2)

   • 这部分代码生成 x 范围在 0 到 10 之间、间隔为 0.25 的数据点，然后计算两个函数 y1=3*x 和 y2=8-x 的对应值。
   • 使用 plot 函数绘制这两个函数的图像。
   • axis([0 10 0 10]) 设置图形坐标轴范围，确保在给定范围内进行绘制。
   • legend('y1=3x','y2=8-x') 添加图例标签。
   • title('蒙特卡洛算法') 添加图形标题。
   • text(2,6,'交点') 在图中添加文本，表示两个函数的交点位置。
   • grid on 打开坐标网格。

2. matlabCopy code

   % 产生随机点 x=unifrnd(0,8,[1,10000000]); y=unifrnd(0,6,[1,10000000]);

   • 这部分代码使用 unifrnd 函数生成 1000 万个在指定区间内的均匀分布的随机点。x 取值在 [0, 8]，y 取值在 [0, 6]。

3. matlabCopy code

   % 统计所在所求图形中的点 frequency=sum(y<3*x & x<=2)+sum(y<8-x & x>2);

   • 统计这些随机点中有多少点落在给定的图形区域内。使用条件 y<3*x & x<=2 和 y<8-x & x>2 分别计算两个函数所对应的区域内的点数。

4. matlabCopy code

   % 计算面积 area=6*8*frequency/10^7;

   • 根据蒙特卡洛算法，通过统计区域内的点数来估计图形的面积。
   • 面积的计算公式为：图形的高（6）乘以宽（8），再乘以统计的点数与总点数的比例。

这个示例演示了如何使用蒙特卡洛方法来估计给定图形的面积。

N=1000000;    %随机点的数目
x=rand(N,1);  %rand 生成均匀分布的伪随机数。分布在（0~1）之间
y=rand(N,1);  %矩阵的维数为N×1
count=0;
for i=1:N
   if (x(i)^2+y(i)^2<=1)
     count=count+1;
    end
end
PI=4*count/N

这段 MATLAB 代码的目的是通过 Monte Carlo 方法估算 π（圆周率）。以下是对每一句代码的解释：

1. matlabCopy code

   N=1000000; %随机点的数目

   • 设置要生成的随机点的数量 N。

2. matlabCopy code

   x=rand(N,1); %rand 生成均匀分布的伪随机数。分布在（0~1）之间 y=rand(N,1); %矩阵的维数为N×1

   • 生成两个 N×1 的列向量 x 和 y，这两个向量包含了均匀分布在 [0, 1) 区间内的随机数。

3. matlabCopy code

   count=0;

   • 初始化一个计数器 count，用于记录落在单位圆内的点的数量。

4. matlabCopy code

   for i=1:N if (x(i)^2+y(i)^2<=1) count=count+1; end end

   • 使用循环遍历所有的随机点 (x(i), y(i))，并检查每个点是否落在单位圆内。如果是，就将计数器 count 增加。

5. matlabCopy code

   PI=4*count/N

   • 最后，通过 Monte Carlo 估算 π 的值。具体计算是通过单位圆的面积占整个正方形的面积的比例，然后乘以 4。这个比例等于单位圆的面积除以正方形的面积，而单位圆的面积为 π。因此，通过 Monte Carlo 方法估算的 π 等于 4 乘以单位圆内的点数除以总的点数 N。

clear
clc
N=10000;  
x=rand(N,1); 
y=rand(N,1);
count=0;
for i=1:N
   if (y(i)<=x(i)^2)
     count=count+1;
   end
end
result=count/N

这段 MATLAB 代码的目的是通过 Monte Carlo 方法估算概率值。以下是对每一句代码的解释：

1. matlabCopy code

   clear clc

   • clear 用于清除当前工作区的所有变量，clc 用于清除命令窗口的内容，以确保开始时工作区变量清空，命令窗口清屏。

2. matlabCopy code

   N=10000;

   • 设置要生成的随机点的数量 N。

3. matlabCopy code

   x=rand(N,1); y=rand(N,1);

   • 生成两个 N×1 的列向量 x 和 y，这两个向量包含了均匀分布在 [0, 1) 区间内的随机数。

4. matlabCopy code

   count=0;

   • 初始化一个计数器 count，用于记录满足某个条件的点的数量。

5. matlabCopy code

   for i=1:N if (y(i)<=x(i)^2) count=count+1; end end

   • 使用循环遍历所有的随机点 (x(i), y(i))，并检查每个点是否满足条件 y <= x^2。如果是，就将计数器 count 增加。

6. matlabCopy code

   result=count/N

   • 最后，通过 Monte Carlo 估算概率值。具体计算是通过满足条件 y <= x^2 的点数除以总的点数 N。这个例子中，估算的概率值即为 result。