蒙特卡羅實驗報告

蒙特卡羅方法

實驗一實驗報告蒙特卡羅方法實驗一實驗報告一、實驗目的1、了解蒙特卡羅方法方法的基本思想；2、掌握蒙特卡羅方法計算面積、體積的方法；3、掌握由已知分布的隨機抽樣方法。二、實驗原理MonteCarlo方法，又稱統(tǒng)計模擬方法或計算機隨機模擬方法，是一種基于“隨機數(shù)”進行數(shù)值模擬的方法，一種采用統(tǒng)計抽樣理論近似求解物理或數(shù)學問題的方法。倘若待求量可以表述成某些特征量的期望值、某些事件出現(xiàn)的概率或兩者的函數(shù)形式，那么可采用蒙特卡羅方法求解。在求解某些特征量的期望值或某些事件出現(xiàn)的概率時，必須構建合符實際的數(shù)學模型。例如采用蒙特卡羅方法計算某函數(shù)所圍面積時，構建的數(shù)學模型是構造一已知面積的可均勻抽樣區(qū)域，在該區(qū)域投點，由伯努利定理大數(shù)定理可知，進入待求區(qū)域投點的頻率依概率1收斂于該事件出現(xiàn)的概率(面積之比)。由已知分布的隨機抽樣方法指的是由已知分布的總體中抽取簡單子樣。抽樣方法有：直接抽樣方法：離散型分布隨機抽樣方法、連續(xù)型分布直接抽樣方法；挑選抽樣方法；復合抽樣方法；隨機抽樣一般方法：加抽樣方法、減抽樣方法、乘抽樣方法、乘加抽樣方法、乘減抽樣方法、對稱抽樣方法、替換抽樣方法、多為分布抽樣方法、積分抽樣方法；隨機抽樣其他方法：偏倚抽樣方法、近似分布抽樣方法、近似-修正抽樣方法。三、實驗內容1、安裝所需計算工具(MATLAB、fortran、C++等)；2、編寫一偽隨機數(shù)發(fā)生器；(如乘加同余a=1366,c=150889,M=714025、a=9301,c=49297,M=233280；乘同余a=16807,M=232-1；或采用其它方法)以下內容選取一個采用自編偽隨機數(shù)發(fā)生器進行計算，其余采用工具軟件中自帶偽隨機數(shù)發(fā)生器進行計算。3、求解以下區(qū)域的面積、體積：3.1、給定曲線y=2-x2和曲線y3=x2，曲線的交點為：P1(-1，1)、P2(1，1)。曲線圍成平面有限區(qū)域，用蒙特卡羅方法計算區(qū)域面積；z>.Jx2+y23.2、計算,所圍體積z<1+J1-x2-y2其中。={(x,y,z)|-1<x<1,-1<y<1,0<z<2}。4、對以下已知分布進行隨機抽樣：4.1、f(x)=2x3+-(1-x)2,xe[0,1]；24.2、f(x)=二-f(x),xeh,2E+1]Ek(E)1〃)_，E+1—x)2111其中(x}=r~E^rJ+x—X2+X3k(E)=[1—2(E+1)]in(2E+1)+1+-—1E22E2(2E+1「四、實驗程序及其相關情況第2題function[SJS]=suiji(ZHONG)%SJS產生的隨機數(shù)%ZHONG輸入的隨機數(shù)種子a=1366;c=150889;M=714025;SJS=mod(ZHONG*a+c,M)/M;第3.1題clear;clc;M=0;%記錄投點在所圍圖形中的個數(shù)N=input('請輸入總投點個數(shù):\n');ksi=0.89656;%用輸入的方式tic;fori=1:N其中ksi=suiji(ksi);x=2*ksi-1;ksi=suiji(ksi);y=2*ksi;ify<=2-xA2ifyA3>=xA2M=M+1;endendendtocMIANJI=M/N*4clearMNixy;計算結果：N=50000時面積為2.1431，計算時間約0.688s。第3.2題clear;clc;M=0;%記錄投點在所圍圖形中的個數(shù)N=input('請輸入總投點個數(shù):\n');tic;fori=1:Nx=2*rand()-1;y=2*rand()-1;z=2*rand();t=xA2+yA2;s=zA2;ifs>=tift<=-s+2*zM=M+1;endendendtocMIANJI=M/N*8clearMNixy;計算結果：N=50000時面積為3.1350，計算時間約0.282s。第4.1題clear;clc;M=input('輸入所需產生隨機變量的個數(shù):\n');x=zeros(M,1);tic;fori=1: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());endendtocclearM;第4.2題clear;clc;E=input('輸入入射光子的能量(單位keV):\n');E=E/511;%計算系數(shù)KE=(1-2*(E+1)/(EA2))*log(2*E+1)+0.5+4/E-1/(2*(2*E+1)A2);p1=(E/(2*E+1))A2/KE;p2=p1+(2/E+2*E/(2*E+1)A2)/KE;p3=p1+p2+(1-2*(E+1)/(EA2))*log(2*E+1)/KE;

tic;a=rand();ifa<=p1Cshe=(2*E+1)/(2*E*max(rand(),rand())-1);elseifa<=p2Cshe=(2*E+1)/(2*E*rand()+1);elseifa<=p3Cshe=(2*E+1)Arand();%可采用倒數(shù)抽樣方法修正elseCshe=2*E*rand()+1;endtocclearE;Cshe=Cshe*0.511Z1，t22dt<xZ1，t22dt<x序號12314567試驗次數(shù)1035x1045x1041.2x1051.5x1051.8x1062.0x107試驗時間計算結果實驗誤差1、3.1和3.2抽樣結果誤差隨抽樣次數(shù)的關系圖，及相關原因;表1實驗記錄表蒙特卡羅方法的近似值與真實值得誤差問題，概率論中心極限定理指出，如果隨機變量序列Z2，。。。，Zn獨立同分布，且具有有限異于零的方差，其中b是隨機變量Z的均方差b2=E（Z-E（Z））2=j（t-E（Z）》f（t）dtF（x）是Z的分布密度函數(shù)。當N充分大時，有如下近似式PZ-E（Z）vk1「t2jxe2dt=1-a-x<2k,其中a稱為置信度，1-a成為置信水平。該式的含義是：在1-F（x）是Z的分布密度函數(shù)。當N充