【Matlab】蒙特卡洛算法

蒙特卡洛算法 入门教学

假设我们有个y=x^2的表达式，如何用MC方法求得函数在[0,1]区间的定积分呢？

定积分可以用面积来求解，也就是通过求箭头下的面积

x=0:0.01:1;y=x.^2;plot(x,y); % 绘图

staus=10;
for i=1:4  % 4次模拟 % 为观察随数据量的增加，蒙特卡洛求出的正确性是否增加
    point=staus.^i; % 模拟的随机点数
    RandData=rand(2,point); % 根据随机点数，产生随机的(x,y)散点,不明白可以试试   % scatter(RandData(1,:),RandData(2,:))
    Below=find(RandData(1,:).^2>RandData(2,:)); % 寻找位于曲线下的散点
    Outcome(i)=length(Below)/length(RandData); % 最终结果的表示
end

结果如下：

Outcome =

   0.200000000000000   0.280000000000000   0.329000000000000   0.337900000000000

从结果看来，通过不断增加随即点数，结果越与真实值相符。

当数据规模达到1e5时，定积分结果趋于稳定：

point=100000; %模拟的规模
RandData=rand(2,point); 
Below=find(RandData(1,:).^2>RandData(2,:));
Outcome=length(Below)/length(RandData);

% 绘图
BelowData=RandData(:,Below);
hold on
scatter(BelowData(1,:),BelowData(2,:))

例1

题目： y = x 2 y=x^2 y=x2,， y = 12 − x y=12-x y=12−x与4x$轴在第一象限围成一个曲边三角形。设计一个随机实验，求该图形面积的近似值。

解析：在矩形区域[0,12]×[0,9]上产生服从均匀分布的1e7个随机点,统计随机点落在曲边三角形的频数，则曲边三角形的面积近似为上述矩形的面积乘以频率。

x=unifrnd(0,12,[1,10000000]);
y=unifrnd(0,9,[1,10000000]);
pinshu=sum(y