Monte Carlo方法的方法的的基本思想是首先建立一个概率模型，使所求问题的解正好是该模型参数或其他有关特征量，然后通过模拟的基本思想是首先建立一个概率模型，使所求问题的解正好是该模型参数或其他有关特征量，然后通过模拟 统计试验，即多次随机抽样试验（确定m和n）,统计出某事件发生的百分比，只要试验次数很大，该百分比就近拟于事件发生的概率。这实际上就是事件发生概率的统计定义。利用建立的概率模型，求出要估计的参数。蒙特卡洛属于试验数学的一个分支。统计试验，即多次随机抽样试验（确定m和n）,统计出某事件发生的百分比，只要试验次数很大，该百分比就近拟于事件发生的概率。这实际上就是事件发生概率的统计定义。利用建立的概率模型，求出要估计的参数。蒙特卡洛属于试验数学的一个分支。Monte Carlo仿真原理Monte Carlo仿真原理MATLAB实现Buffon问题仿真求解程序MATLAB实现Buffon问题仿真求解程序clear all; L=1; %针的长度；针的长度； d=2; %平行线间的距离（平行线间的距离（dL）；）； m=0; %统计满足针与线相交条件的次数并赋初值；统计满足针与线相交条件的次数并赋初值； n=10000; %投针试验次数投针试验次数 for k=1:n %迭代次数迭代次数 x=unifrnd(0,d/2); %随机产生数的长度，即投针之后针中点与平行线的距离随机产生数的长度，即投针之后针中点与平行线的距离 p=unifrnd(0,pi); %随机产生的针与线相交的角度随机产生的针与线相交的角度 if xL） m=0; %统计满足针与线相交条件的次数并赋初值统计满足针与线相交条件的次数并赋初值 n=10000; %投针试验次数投针试验次数 for k=1:n %迭代次数迭代次数 x=unifrnd(0,d/2); %随机产生数的长度，即投针之后针中点与平行线的距离随机产生数的长度，即投针之后针中点与平行线的距离 p=unifrnd(0,pi); %随机产生的针与线相交的角度随机产生的针与线相交的角度 if x0.5%模拟整个赌博过程模拟整个赌博过程 r=(rand=y %落在阴影区域内点的频数落在阴影区域内点的频数 Sn(i)=n/N(i); %落到阴影区域内点的频率，即概率的模拟值落到阴影区域内点的频率，即概率的模拟值 end S0 S0 =0.1667 Sn =0 0.1600 0.1620 0.1742 0.1676 0.1667 运行结果运行结果教材 P214 例：计算定积分教材 P214 例：计算定积分2 12 01 2x Jedx=源程序源程序1如下：运行结果如下：运行结果clear all;close all; N=10,100,1000,10000,100000,1000000,10000000; %输入参数 N 是随机投点的个数 syms x;S0=int(exp(-x2/2)/sqrt(2*pi),0,1); %面积的理论值（解析解） S0=vpa(S0,6);%vpa 函数得出确定小数点位数后的确切解 %S0=quad(x)exp(-x2/2)/sqrt(2*pi),0,1); %面积的理论值（数值解） %计算阴影区域的面积的 Monte Carlo 模拟值 for i=1:length(N) x=unifrnd(0,1,N(i),1); %点的横坐标 y=unifrnd(0,1,N(i),1); %点的纵坐标 n=sum(exp(-x(i)2/2)/sqrt(2*pi)=y); %落在阴影区域内点的频数 Sn(i)=vpa(n/N(i),6); %落到阴影区域内点的频率，即概率的模拟值 end S0 clear all;close all; N=10,100,1000,10000,100000,1000000,10000000; %输入参数 N 是随机投点的个数 syms x;S0=int(exp(-x2/2)/sqrt(2*pi),0,1); %面积的理论值（解析解） S0=vpa(S0,6);%vpa 函数得出确定小数点位数后的确切解 %S0=quad(x)exp(-x2/2)/sqrt(2*pi),0,1); %面积的理论值（数值解） %计算阴影区域的面积的 Monte Carlo 模拟值 for i=1:length(N) x=unifrnd(0,1,N(i),1); %点的横坐标 y=unifrnd(0,1,N(i),1); %点的纵坐标 n=sum(exp(-x(i)2/2)/sqrt(2*pi)=y); %落在阴影区域内点的频数 Sn(i)=vpa(n/N(i),6); %落到阴影区域内点的频率，即概率的模拟值 end S0 SnSn S0 =.341344 Sn =.400000 .330000 .354000 .389800 .398080 .313488 .342588 源程序源程序2如下：如下：clear all;close all; 运行结果运行结果N=10,100,1000,10000,100000,1000000,10000000; %输入参数输入参数 N 是随机投点的个数是随机投点的个数 syms x;S0=int(exp(-x2/2)/sqrt(2*pi),0,1); %面积的理论值（解析解）面积的理论值（解析解） S0=vpa(S0,6);%vpa 函数得出确定小数点位数后的确切解函数得出确定小数点位数后的确切解 %S0=quad(x)exp(-x2/2)/sqrt(2*pi),0,1); %面积的理论值（数值解）面积的理论值（数值解） %计算阴影区域的面积的计算阴影区域的面积的 Monte Carlo 模拟值模拟值 for i=1:length(N) x=unifrnd(0,1,N(i),1); %点的横坐标点的横坐标 y=unifrnd(0,1,N(i),1); %点的纵坐标点的纵坐标 n=sum(exp(-x.2/2)/sqrt(2*pi)=y); %落在阴影区域内点的频数落在阴影区域内点的频数 Sn(i)=vpa(n/N(i),6); %落到阴影区域内点的频率，即概率的模拟值落到阴影区域内点的频率，即概率的模拟值 end S0 Sn S0 =.341344 Sn = .300000, .340000, .334000, .346200, .343860, .339903, .341463 对于一般区间对于一般区间 , a b上的定积分 上的定积分 ( )baJg x dx =作线性变换作线性变换xayba=，即可化成，即可化成0,1区间上的积分. 进一步若区间上的积分. 进一步若( )cg xd,可令 ,可令 1( ) () f yg aba ycdc=+则则0( )1f y. 此时有 . 此时有 10( )()()( )()baJg x dxba dcf y dyc ba =+用上述 Monte Carlo 方法计算上式中用上述 Monte Carlo 方法计算上式中10( )f y dy，代入上式即可求出 ，代入上式即可求出 ( )baJg x dx =例教材P218 习题16 计算定积分例教材P218 习题16 计算定积分121xJe dx =作线性变换作线性变换1 2xaxyba+=，即化成，即化成0,1区间上 的积分，但区间上 的积分，但11( )xceg xeed=，故令 ，故令 1( ) () f yg aba ycdc=+1 11exp(21)yeee =即即0( )1f y，此时有 ，此时有 11111 211012()exp(21)2xJe dxeeyedyeee =+按上述Monte carlo方示计算出上述积分项，再代入即可按上述Monte carlo方示计算出上述积分项，再代入即可编写程序实现编写程序实现%J2 的随机投点法计算的随机投点法计算 clear all;close all; N=10,100,1000,10000,100000,1000000,10000000; %输入参数输入参数 N 是随机投点的个数是随机投点的个数 %syms x;S0=int(exp(x),-1,1); %面积的理论值（解析解）面积的理论值（解析解） %S0=vpa(S0,6);%vpa 函数得出确定小数点位数后的确切解函数得出确定小数点位数后的确切解 syms x;S0=quad(x)(exp(x),-1,1); %面积的理论值（数值解）面积的理论值（数值解） S0=vpa(S0,6);%vpa 函数得出确定小数点位数后的确切解函数得出确定小数点位数后的确切解 %计算积分项阴影区域的面积的计算积分项阴影区域的面积的 Monte Carlo 模拟值模拟值 for i=1:length(N) x=unifrnd(0,1,N(i),1); %点的横坐标点的横坐标 y=unifrnd(0,1,N(i),1); %点的纵坐标点的纵坐标 n=sum(exp(-1+2*x.1)-exp(-1)/(exp(1)-exp(-1)=y); %计算上述式子积分项中被积函数落在阴影区域内点的频数计算上述式子积分项中被积函数落在阴影区域内点的频数 Sn(i)=vpa(n/N(i),6); %上述积分项落到阴影区域内点的频率，即概率的模拟值上述积分项落到阴影区域内点的频率，即概率的模拟值 end S0 Sn %输出积分项值输出积分项值 S=vpa(Sn.1*2*(exp(1)-exp(-1)+2*exp(-1),6) %将上述式子中的积分项代回将上述式子中的积分项代回 J2中，并取小数点后中，并取小数点后 6 位值位值 S0 =2.35040 Sn = .500000, .340000, .312000, .338300, .344010, .343157, .343434S = 3.08616, 2.33403, 2.20241, 2.32604, 2.35288, 2.34887, 2.35018 运行结果运行结果2. 用Monte Carlo方法计算重积分2. 用Monte Carlo方法计算重积分见文献 ： 谢中华MATLAB统计分析与应用40个案例分析北京航空航天大学出版社2010.6 P145页见文献 ： 谢中华MATLAB统计分析与应用40个案例分析北京航空航天大学出版社2010.6 P145页用Monte Carlo方法计算积分(用Monte Carlo方法计算积分(平均值法平均值法) ) 2. 用Monte Carlo方法计算定积分2. 用Monte Carlo方法计算定积分( (P216 教材 例教材 例4.2.4 ) ) 用用Monte Carlo方法计算定积分（方法计算定积分（平均值法平均值法）为计算定积分 为计算定积分 10( )Jf x dx=设设(0,1)XUniform，则，则()Yf X=的数学期望 的数学期望 10 ()( )E f Xf x dxJ=所以估计所以估计J的值就是估计的值就