蒙特卡洛

简述

　　蒙特卡罗方法又称随机抽样技巧或统计试验方法。半个多世纪以来，由于科学技术的发展和电子计算机的发明 ，这种方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。蒙特卡罗方法是一种计算方法，但与一般数值计算方法有很大区别。它是以概率统计理论为基础的一种方法。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。

基本思想

　　当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。这就是蒙特卡罗方法的基本思想。

因此，可以通俗地说，蒙特卡罗方法是用随机试验的方法计算积分，即将所要计算的积分看作服从某种分布密度函数f(r)的随机变量ｇ(r)的数学期望

                       

通过某种试验，得到Ｎ个观察值r1，r2，…，rN（用概率语言来说，从分布密度函数f(r)中抽取Ｎ个子样r1，r2，…，rN，），将相应的Ｎ个随机变量的值g(r1)，g(r2)，…，g(rN)的算术平均值 作为积分的估计值（近似值）。

                        

优点 能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。 受几何条件限制小。 收敛速度与问题的维数无关。 具有同时计算多个方案与多个未知量的能力。 误差容易确定。 程序结构简单，易于实现。

缺点 收敛速度慢。 误差具有概率性。 在粒子输运问题中，计算结果与系统大小有关。

实例

蒲丰投针

　　设针投到地面上的位置可以用一组参数（x,θ）来描述，x为针中心的坐标，θ为针与平行线的夹角，如图所示。 任意投针，就是意味着x与θ都是任意取的，但x的范围限于［0，a］，夹角θ的范围限于［0，π］。在此情况下，针与平行线相交的数学条件是

       

 

x在［0，a］上任意取值，表示x在［0，a］上是均匀分布的，其分布密度函数为：

 

类似地，θ的分布密度函数为：

 eg：用蒲丰投针法在计算机上计算π值，取a=4、l=3。

%clear;clc;
N=10000;a=4;l=3;counter=0;
Rand1=unifrnd(0,a,1,N);%产生n个(0,a)之间均匀分布bai的随机数，du这里a/2是投针的中点到最近的平行线的距离
Rand2=unifrnd(0,pi,1,N);% 产生n个(0,pi)之间均匀分布的随机数，这里pi是投针到最近的平行线的角度
for i=1:N
    if Rand1(i)<l*sin(Rand2(i))
        counter=counter+1;
    end
end
frequency=counter/N;
PI=2*l/(a*frequency)

 

 eg：求连续掷两颗骰子，点数之和大于6且第一次掷出的点数大于第二次掷出点数的概率。

理论计算为9/36。

clear;clc;
N=1000000;counter=0;
Rand1=unidrnd(6,1,N);
Rand2=unidrnd(6,1,N);
for i=1:N
    if (Rand1(i)+Rand2(i))>6 && Rand1(i)>Rand2(i)
        counter=counter+1;
    end
end
frequency=counter/N

结果为0.2508。与理论非常相近。

eg：

在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏．为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点．经过长期观察发现，我方指挥所对敌方目标的指示有50％是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部毁伤敌人火炮．现在希望能用某种方式把我方将要对敌人实施的次打击结果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值。使用蒙特卡洛方法模拟次打击结果：

function [out1,out2,out3,out4] = Msc(N)
%UNTITLED8 此处显示有关此函数的摘要
%   此处显示详细说明
k1=0;k2=0;k3=0;
for i=1:N
    x0=randperm(2)-1;
    y0=x0(1);
    if y0==1
        fprintf('第%d次：指示正确||',i);
        x1=randperm(6);
        y1=x1(1);
        if y1==1||y1==2||y1==3
            fprintf('第%d次：击中0炮||',i);
            k1=k1+1;
        else
            if y1==4||y1==5
                fprintf('第%d次：击中1炮||',i);
                k2=k2+1;
            else
                fprintf('第%d次：击中2炮||',i);
                k3=k3+1;
            end
        end
    else
        fprintf('第%d次：指示错误，击中0炮||',i);
    end
    fprintf('\n');
    out1=(k2+k3)/N;  %有效打击的概率
    out2=k2/N;   %每次击中摧毁一门火炮的概率
    out3=k3/N;     %每次击中摧毁两门火炮的概率
    out4=(k2+k3*2)/N;   %平均每次击中的火炮个数
end

 

 

 

 a:有效打击的概率  b:每次击中摧毁一门火炮的概率   c:每次击中摧毁两门火炮的概率  d:平均每次击中的火炮个数