蒙特卡罗实验报告.docx
蒙特卡罗方法试验一实验报告蒙特卡罗方法试验一试验报告一、试验目的1、了解蒙特卡罗方法方法的基本思想；2、把握蒙特卡罗方法计算面积、体积的方法;3、把握由已知分布的随机抽样方法。二、试验原理Monte Carlo方法，又称统计模拟方法或计算机随机模拟方法，是一种基于“随机数”进行数值模拟的方法，一种采纳统计抽样理论近似求解物理或数学问题的方法。如果待求量可以表述成某些特征量的期望值、某些大事消失的概率或两者的函数形式, 那么可采纳蒙特卡罗方法求解。在求解某些特征量的期望值或某些大事消失的概率时，必需 构建合符实际的数学模型。例如采纳蒙特卡罗方法计算某函数所围面积时，构建的数学模型 是构造一已知面积的可匀称抽样区域，在该区域投点，由伯努利定理大数定理可知，进入待 求区域投点的频率依概率1收敛于该大事消失的概率(面积之比)。由已知分布的随机抽样方法指的是由已知分布的总体中抽取简洁子样。抽样方法有：直接抽样方法：离散型分布随机抽样方法、连续型分布直接抽样方法；选择抽样方法; 复合抽样方法；随机抽样一般方法：加抽样方法、减抽样方法、乘抽样方法、乘加抽样方法、 乘减抽样方法、对称抽样方法、替换抽样方法、多为分布抽样方法、积分抽样方法；随机抽 样其他方法：偏倚抽样方法、近似分布抽样方法、近似-修正抽样方法。三、试验内容1、安装所需计算工具(MATLAB、fortran> C+等)；2、编写一伪随机数发生器；(如乘加同余a=1366,c=150889,M=714025、 a=9301 ,c=49297,M=233280；乘同余 a= 16807,M=232-l ；或采纳其它方法)以下内容选取一个采纳自编伪随机数发生器进行计算，其余采纳工具软件中自带伪随机 数发生器进行计算。3、求解以下区域的面积、体积：3.1、 给定曲线y=2-x2和曲线y3=x2,曲线的交点为：Pi(-1, 1 )> P2( 1, 1 )o曲线 围成平面有限区域，用蒙特卡罗方法计算区域面积；z > Jx2 + y23.2、 计算17 二所围体积Z < + J _%2 _ y2其中。=(x, y, 2) | 1 < x < 1, 1 < y < L 0 < z < 2。4、对以下已知分布进行随机抽样：394.1 > /(%) = 2%3 + (1 -x)0,1 ；4.2、人（力看"G）”L2e+i/(%) =其中<k(E)=E + l-x? 1Ex1 1 1-I123XXX2(£ + l)E2-ln(2E + l)1 4+ +2 E2(2E + 1)2四、试验程序及其相关状况第2题function SJS=suiji(ZHONG)%SJS产生的随机数% ZHONG输入的随机数种子a= 1366;c=l 50889;M=714025;SJS=mod(ZHONG*a+c,M)/M;第31题clear;clc;M=0;%纪录投点在所围图形中的个数N= input，请输入总投点个数:n);ksi=0.89656; %用输入的方式tic;for i=l:Nksi=suiji(ksi);x=2*ksi-l;ksi=suiji(ksi);y=2*ksi;if y<=2-xA2ifyA3>=xA2M=M+1;endendendtoeMIANJI=M/N*4clear M N i x y;计算结果：N=50000时面积为2.1431,计算时间约0.688s。第3.2题clear;clc;M=0;%纪录投点在所围图形中的个数N= input，请输入总投点个数:n);tic;for i=l:Nx=2*rand()-l;y=2*rand()-l;z=2*rand();t=xA2+yA2;s=zA2;if s>=tif t<=-s+2*zM=M+1;endendendtoeMIANJI=M/N*8clear M N i x y;计算结果：250000时面积为3.1350,计算时间约0.282s。第4.1题clear;clc;M = input。输入所需产生随机变量的个数:n);x = zeros(M,l);tic;for i=l:Mif(rand()<=0.5)x(i) = max(rand(),rand()x(i) = max(x(i),rand();x(i) = max(x(i),rand();elsex(i) = min(rand(),rand();x(i) = min(x(i),rand();endendtoeclear M;第4.2题clear;clc;E = input。输入入射光子的能量(单位keV):nf);E=E/511;%计算系数KE=( 1 -2*(E+1 )/(E 八 2)*log(2*E+1 )+0.5+4/E- 1/(2*(2*E+ 1)A2);pl=(E/(2*E+l)A2/KE;p2=pl +(2/E+2*E/(2*E+ 1)A2)/KE;p3=p 1 +p2+( 1 -2*(E+1 )/(E 八 2)*log(2*E+1 )/KE;tic; a=rand(); if a<=plCshe = (2*E+1 )/(2*E*max(rand(),rand()-1);elseif a<=p2Cshe = (2*E+l)/(2*E*rand()+l);elseif a<=p3Cshe = (2*E+1)八rand();%可采纳倒数抽样方法修正elseCshe = 2*E*rand()+l; end clear E;Cshe = Cshe*0.511五、试验结果与分析1、3.1和3.2抽样结果误差随抽样次数的关系图，及相关缘由;表1试验纪录表蒙特卡罗方法的近似值与真实值得误差问题，概率论中心极限定理指出，假如随机变量序列序号12314567试验次数1035xl045xl041.2X1051.5X1051.8X1062.0x107试验时间计算结果试验误差Zh Z2,ooo, Zn独立同分布，且具有有限异于零的方差，则t2：Qdt其中。是随机变量Z的均方差<72=E(Z-E(Z)2=J(t-E(Z)2f (t) dtF (x)是Z的分布密度函数。当N充分大时，有如下近似式P ZN-E(Z)P ZN-E(Z)1-a其中a称为置信度，1-二成为置信水平。该式的含义是：在1-。置信水平下，近似值与真值的误差为值的误差为XCTO依据问题的要求确定出置信水平后，查标准正态分布表，确定出X,就AY ZT可以得到近似值和真值的误差Zn-石(Z) <源这样我们很简洁发觉误差的大小在ox确定的状况下，N越大误差越小2、给出4.1和4.2的抽样框图、试验累积频率与理论累积频率关系图，并给出抽样次 数(>1()6)与抽样时间。
