遇见蒙特卡洛

　　昨天，呃，不，是今天凌晨，我问一位大牛有没有什么既简单又强势的算法，他说了蒙特卡洛。今天查了些资料，见识到它的强大与及简洁。参考了这篇文章，和维基百科。

现在用MATLAB实现蒙特卡洛方法的几个应用。

　　1.计算圆周率

status = 10;

for i=1:6                                   %6次模拟
    point=status.^i;                    %模拟的随机点数
    RandData=rand(2,point);       % 根据随机点数，产生随机的(x,y)散点
    x = RandData(1,:);                  
    y = RandData(2,:);
    inner=find(x.^2+y.^2 <= 1);  % find返回符合条件（在圆内）的点的索引值，构成数组
    Outcome(i)=length(inner)/length(RandData);    % 计算圆内的点占所有点的比例，即为pi/4
    my_pi(i) = Outcome(i)*4;
end
    my_pi
   
    

　　输出如下：

　　可以看到，随着随机点数的增加，算得的pi精度越来越高。（注意这里，我将脚本文件命名为pi.m，系统会有冲突警告，所以最好改下文件名。）

　　2.计算积分

　　这里以计算y=x^2在[0,1]上的定积分为例。

status = 10;

for i=1:4                         %4次模拟
    point=status.^i;          %模拟的随机点数
    RandData=rand(2,point); % 根据随机点数，产生随机的(x,y)散点
    x = RandData(1,:);
    y = RandData(2,:);
    %scatter(x,y)
    Below=find(x.^2>y);      % y=x^2曲线下的点（的下标）
    Outcome(i)=length(Below)/length(RandData);
    subplot(2,2,i)
    scatter(x(Below),y(Below),'r')
end
    axis([0,1,0,1])
    Outcome
    

　　输出如下：

 

　　3.计算测度

　　这里参考一篇论文，讲解关于测度的MATLAB实现。

 　　相关代码如下：

clear all
n = input('please input the dimensionn n of variable=:');
m=0;
for i=1:n
    x=rand(1);
    y=rand(1);
    z=rand(1);
    if x^2+sin(y)<=z
        if x-z+exp(y)<=1
            m=m+1;
        end
    end
end
ration=m/n

　　运行结果如下：

 　　可见结果运行还是很稳定的。

　　4.计算非线性规划的优化解

　　代码如下：

clear all
L=10000;    %要求最小值，设置足够大的上限
M = input('set the number of dots M=')
for N=1:M
    % 将不等关系巧妙地转化为乘积关系
    x2=-9+13*rand(1);
    x1=((x2)^2+9)*rand(1);
    x3=3*(x1)+(x2)-5+(12-(x1)*(x2)-3*(x1)-(x2))*rand(1);
    if 12-(x1)*(x2)-3*x1-x2>=0
        f=3*(x1)^2-4*(x2)^2+5*(x1)*(x3);
        if f<L
            L=f;
        end
    end
end
L

　　运行输出结果：

　　　并且可行域和目标函数的复杂性不对MC方法造成影响。

　　　最后，看到知乎上有句话总结了蒙特卡洛方法——尽量找好的，但不保证是最好的。多么值得思考的一句话啊。