Matlab学习系列18-蒙特卡罗方法

18.蒙特卡罗方法（一）概述一、原理蒙特卡罗（MonteCarlo）方法，是一种基于“随机数”的计算机随机模拟方法，通过大量随机试验，利用概率统计理论解决问题的一种数值方法。其理论依据是：大数定律、中心极限定理（估计误差）。常用来解决如下问题：1.求某个事件的概率，基于“频率的极限是概率”；2.可以描述为“随机变量的函数的数学期望”的问题，用随机变量若干个具体观察值的函数的算术平均值代替。一般形式为求积分：()()()dbaEfXfxxxX为自变量（随机变量），定义域为[a,b],f(x)为被积函数，()x为概率密度函数。蒙特卡罗法步骤为：(1)依据概率分布()x不断生成N个随机数x,依次记为x1,…,xN,并计算f(xi);(2)用f(xi)的算术平均值1()NiifxN近似估计上述积分类比：“积分”同“求和”，“dx”同“1/N”，“()x”同“服从()x分布的随机数”；(3)停止条件：至足够大的N停止；或者误差小于某值停止。3.利用随机数模拟各种分布的随机现象，进而解决实际问题。二、优缺点优点：能够比较逼真地描述具有随机性质的事物的特点及物理实验过程；受几何条件限制小；收敛速度与问题的维数无关；误差容易确定。缺点：收敛速度慢；误差具有概率性；进行模拟的前提是各输入变量是相互独立的。三、应用随机模拟实验，随机最优化问题，含有大量不确定因素的复杂决策系统进行风险模拟分析（金融产品定价、期权）。（二）用蒙特卡罗法求事件概率一、著名的“三门问题”源自博弈论的一个数学游戏：参赛者面前有三扇关闭的门，其中一扇门的后面藏有一辆汽车，而两扇门的后面各藏有一只山羊。参赛者从三扇门中随机选取一扇，若选中后面有车的那扇门就可以赢得该汽车。当参赛者选定了一扇门，但尚未开启它的时候，节目主持人会从剩下的两扇门中打开一扇藏有山羊的门，然后问参赛者要不要更换自己的选择，选取另一扇仍然关着的门。该问题涉及到的问题是：参赛者更换自己的选择是否会增加赢得汽车的概率？解：这是全概率问题。记结果事件B=“赢得汽车”，造成结果的原因事件有两个：A1=“第一次选到汽车”，A2=“第一次选到山羊”则P(A1)=1/3,P(A2)=2/3;P(B|A1)=0,P(B|A2)=1.利用全概率公式，P(B)=P(A1)×P(B|A1)+P(A2)×P(B|A2)=2/3而参赛者不更换选择，抽中汽车的概率为1/3.可见，参赛者更换选择可以增加一倍的获奖几率。下面用蒙特卡罗方法来模拟验证上面的理论结果：将问题“随机数”化，对羊编号为1,2；对汽车编号为3.先从1,2,3中随机选取一个数，若开始选中1或2，则更换选择后选中3,即赢得汽车；若开始选中3,则更换选择后选中1或2，即得不到汽车。这样的试验重复n次，记录开始选中1或2的次数m（即更换选择后赢得汽车的次数），从而可以确定更换选择后赢得汽车的频率m/n.由伯努利大数定律，当试验次数n增大时，频率m/n将趋于更换选择后赢得汽车的概率。代码：先编写SheepAndCar.m函数functionp=SheepAndCar(n)%n可以是正整数，也可以是正整数的向量fori=1:length(n)x=randsample(3,n(i),'true');%从1,2,3中随机抽取n(i)个数,'true'表示可以重复p(i)=sum(x~=3)/n(i);end再执行代码：（分别模拟10次、100次、……、1000000次）p=SheepAndCar([10,100,1000,10000,100000,1000000])运行结果：p=0.60000.72000.69100.66410.66710.6666二、用蒲丰投针法求圆周率π【蒲丰投针问题】：平面上画有间隔为d的等距平行线，向平面内任意投掷一枚长为h（hd）的针，求针与任一平行线相交的概率。用Y表示针的中点与最近一条线的距离，用X表示针与此直线间的夹角，则(X,Y)为二维随机变量，其样本空间为平面上的矩形区域：(,):0,02dxyxy由于是任意投掷，(X,Y)在Ω上服从均匀分布。记事件A=“针与平行线相交”，则A所对应的区域为(,):sin2hAxyyx由几何概率知，事件A的概率为0(/2)sind)2()()/2hxxSAhPASdd（若h,d已知，代入上式就可求出P(A).反过来，由P(A)的值也可以求2()hdPA.P(A)可以通过蒙特卡罗方法随机模拟获得：在区域Ω随机均匀投点，落在区域A中的点的频率fA，随着投点次数的增大，就趋于概率P(A).从而2Ahdf代码：先编写BuffonMonteCarlo.m函数function[p0,pm,piv]=BuffonMonteCarlo(d,h,N)%返回p0为理论概率,pm为基于蒙特卡罗法的模拟概率,piv为pi的模拟值.%输入参数d为相邻两条平行线的间距，h为针的长度，N为模拟投针的次数.ifh=derror('针的长度应小于相邻平行线的间距')endp0=2*h/(d*pi);%计算针与任一平行线相交的理论概率x=0;y=0;m=length(N);pm=zeros(1,m);pival=pm;%赋变量初值fori=1:mx=unifrnd(0,pi,1,N(i));%生成[0,pi]上服从均匀分布的1×N(i)个随机数y=unifrnd(0,d/2,1,N(i));yb=h*sin(x)/2;pm(i)=sum(y=yb)/N(i);%频率pival(i)=2*h/(d*pm(i));%求圆周率的模拟值end再执行代码：[p0,pm,pival]=BuffonMonteCarlo(10,8,[10,100,1000,10000,…100000,1000000])运行结果：p0=0.5093pm=0.50000.53000.50200.51180.51120.5093pival=3.20003.01893.18733.12623.13003.1416（三）用蒙特卡罗法求积分蒙特卡罗法计算四重以下的积分，效率可能不如一般数值积分方法，但蒙特卡罗法与积分重数无关，所以在计算高重积分特别有效。蒙特卡罗法求n重积分的原理：设D为Rn中的区域，f：D→R，其n重积分表示为()dDIfxx先找一个包含区域D的矩体C，在C中随机均匀地投n个点，统计落入D中的点，记为k个，则区域D的测度()()kmCmDn函数f在D上的均值为1()DiiDffxk由于是均匀投点，故有()()()DiiDmCImDffxn注：C为矩体，m(C)即为其各边长的乘积。例1求曲线yx与直线y=x所围成的阴影部分的面积。其理论值为101()d6Sxxx取D=C=[0,1].代码：f=@(x)sqrt(x)-x;I0=quad(f,0,1)%理论值(数值解)n=10000;X=unifrnd(0,1,1,n);formatlongI=sum(f(X))/n运行结果：I0=0.166659569272020I=0.166386907074713例2求球体2224xyz被圆柱面222xyx所截得的立体（含在圆柱面内的部分）的体积。如图为立体在第一卦限的部分，在xoy面的投影区域为D，转化为极坐标2(,):02,02Dxyxyxx2244dd9.644DVxyxy取C=[0,2]×[0,1],代码：I0=4*quad2d(@(x,y)sqrt(4-x.^2-y.^2),0,2,0,…@(x)sqrt(1-(1-x).^2))%理论值(数值解)n=100000;x1=unifrnd(0,2,1,n);x2=unifrnd(0,1,1,n);ind=(x2=sqrt(2*x1-x1.^2));formatlongI=4*2*sum(sqrt(4-x1(ind).^2-x2(ind).^2))/n运行结果：I0=9.644049708032036I=9.644614410855425例3计算积分1121122221233211dddxxxxxxxxxxxx已知数值结果为179.2969.取C=[1,2]×[1,4]×[1,16].代码：n=100000;x1=unifrnd(1,2,1,n);x2=unifrnd(1,4,1,n);x3=unifrnd(1,16,1,n);ind=(x2=x1)&(x2=2*x1)&(x3x1.*x2)&(x32*x1.*x2);formatlongI=(4-1)*(16-1)*sum(x1(ind).*x2(ind).*x3(ind))/n运行结果：I=1.790212760006343e+002