蒙特卡羅Monte Carlo模擬誤差分析課程設(shè)計哈工大誤差原理課程設(shè)計

1、Monte Carlo模擬誤差分析課程設(shè)計Monte Carlo模擬誤差分析課程設(shè)計1. 實驗?zāi)康?.1 學(xué)習(xí)并掌握MATLAB軟件的基本功能和使用。1.2 學(xué)習(xí)并掌握基于Monte Carlo Method(MCM)分析的不確定度計算方法。1.3 研究Guide to the expression of Uncertainty in Measurement(GUM)法與MCM法的區(qū)別與聯(lián)系和影響因素，自適應(yīng)MCM方法，基于最短包含區(qū)間的MCM法。2. MATLAB軟件介紹實驗內(nèi)容2.1 介紹MATLAB軟件的基本知識MATLAB名字由MATrix和LABoratory 兩詞的前三個字母組合而

2、成。20世紀(jì)七十年代，時任美國新墨西哥大學(xué)計算機(jī)科學(xué)系主任的Cleve Moller出于減輕學(xué)生編程負(fù)擔(dān)的動機(jī)，為學(xué)生設(shè)計了一組調(diào)用LINPACK和EISPACK矩陣軟件工具包庫程序的“通俗易用”的接口，此即用FORTRAN編寫的萌芽狀態(tài)的MATLABMATLAB語言的主要特點(diǎn)(1). 具有豐富的數(shù)學(xué)功能(2). 具有很好的圖視系統(tǒng)(3). 可以直接處理聲言和圖形文件。(4). 具有若干功能強(qiáng)大的應(yīng)用工具箱。(5). 使用方便，具有很好的擴(kuò)張功能。(6). 具有很好的幫助功能演示內(nèi)容：(1). MATLAB的數(shù)值計算功能在“命令行”Command提示窗口中鍵入：“A=eye(5,5);A=ze

3、ros(5,5);A=ones(5,5)”等命令生成各類矩陣；在“命令行”Command提示窗口中鍵入：“v,d=eig (A)”生成特征矩陣和特征向量；在“命令行”Command提示窗口中鍵入：“expm(A)”對矩陣A求冪；在“命令行”Command提示窗口中鍵入：x=1 3 5;y=2 4 6;z=conv(x,y)；顯示結(jié)果：z = 2 10 28 38 30(2). MATLAB的符號計算功能在“命令行”Command提示窗口中鍵入：syms a x；f=sin(a*x); df=diff(f,x); dfa=diff(f,a)；Command提示窗口顯示結(jié)果： df =cos(a*

4、x)*a； dfa =cos(a*x)*x；2.2 MATLAB軟件畫圖特性(1). MATLAB二維繪圖命令函數(shù)：plot 參數(shù)：線型、顏色、多重線、網(wǎng)格和標(biāo)記、畫面窗口分割、其他方式、隱函數(shù)的描繪)(2). MATLAB三維畫圖曲面與網(wǎng)格圖命令函數(shù)：mesh三維帶陰影曲面圖：surf三維曲線命令：plot3演示內(nèi)容：(1). MATLAB的二維繪圖功能在命令行Command提示窗口中鍵入：close all; x=linspace(0, 2*pi, 100); % 100個點(diǎn)的x座標(biāo) y=sin(x); % 對應(yīng)的y座標(biāo) plot(x,y); 得到如下的結(jié)果：圖1在命令行Command提示

5、窗口中鍵入：“plot(x, sin(x), x, cos(x);”得到如下的結(jié)果：圖2在命令行Command提示窗口中鍵入：plot(x, sin(x), 'co', x, cos(x), 'g*'); 得到如下的結(jié)果：圖3在命令行Command提示窗口中鍵入：xlabel('Input Value'); % x軸注解 ylabel('Function Value'); % y軸注解 title('Two Trigonometric Functions'); % 圖形標(biāo)題 legend('y = sin(

6、x)','y = cos(x)'); % 圖形注解 grid on; % 顯示格線 得到如下的結(jié)果：圖4(2). MATLAB的多維繪圖功能在命令行Command提示窗口中鍵入：X,Y = meshgrid(-3:0.125:3); % 生成二維網(wǎng)格點(diǎn)Z = peaks(X,Y); % 生成某種內(nèi)置函數(shù)mesh(X,Y,Z);得到如下的結(jié)果圖5其他的演示功能詳見“MATLAB畫圖文檔”3. Monte Carlo模擬誤差分析的實驗原理在誤差分析的過程中，常用的方法是通過測量方程推導(dǎo)出誤差傳遞方程，再通過不確定度的合成公式獲得間接測量量的標(biāo)準(zhǔn)不確定度和擴(kuò)展不確定度(GUM

7、)。在有些場合下，測量方程較難獲得，在這種情況下研究誤差的特性就需要借助于模擬統(tǒng)計的方式進(jìn)行計算。Monte Carlo(MCM)法就是較為常用的數(shù)學(xué)工具，具體原理相見相關(guān)資料。此次課程設(shè)計中按照實驗要求產(chǎn)生的隨機(jī)數(shù)可以模擬測量誤差，通過對這些隨機(jī)數(shù)的概率密度分布函數(shù)的面積、包絡(luò)線和概率特征點(diǎn)的求取，可以獲得隨機(jī)誤差的標(biāo)準(zhǔn)不確定度(MCM)，并與理論上估計標(biāo)準(zhǔn)不確定度的Bessel公式、極差法作(GUM)比較，完成實驗內(nèi)容。并以此作為基礎(chǔ)，分析GUM法與MCM法的區(qū)別與聯(lián)系，影響MCM法的參數(shù)，自適應(yīng)MCM法和基于最短包含區(qū)間的MCM法。已知兩項誤差分量服從正態(tài)分布，標(biāo)準(zhǔn)不確定度分別為mV，

8、mV，用統(tǒng)計模擬分析法給出兩項誤差和的分布(誤差分布的統(tǒng)計直方圖，合成的標(biāo)準(zhǔn)差，合成的置信概率 P為99.73%的擴(kuò)展不確定度)。4. Monte Carlo模擬誤差分析的實驗內(nèi)容4.1 MCM法與GUM法進(jìn)行模擬誤差分析和不確定度計算(1). 利用MATLAB軟件生成0，1區(qū)間的均勻分布的隨機(jī)數(shù)；(2). 給出誤差分量的隨機(jī)值：利用MATLAB，由均勻分布隨機(jī)數(shù)生成標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)，誤差分量隨機(jī)數(shù)可表示為mV；同理得 mV(3). 求和的隨機(jī)數(shù)：誤差和的隨機(jī)數(shù)；(4). 重復(fù)以上步驟，得誤差和的隨機(jī)數(shù)系列： ；(5). 作誤差和的統(tǒng)計直方圖：以誤差數(shù)值為橫坐標(biāo)，以頻率為縱坐標(biāo)作圖。作圖區(qū)間

9、應(yīng)包含所有數(shù)據(jù)，按數(shù)值將區(qū)間等分為組（盡可能大），每組間隔為，記數(shù)各區(qū)間的隨機(jī)數(shù)的數(shù)目，以為底，以為高作第()區(qū)間的矩形，最終的組矩形構(gòu)成誤差和的分布直方圖，該圖包絡(luò)線線即為實驗的誤差分布曲線。(6). 以頻率為界劃定區(qū)間，該區(qū)間半寬即為測量總誤差的置信概率為99.73%的擴(kuò)展不確定度。(7). 合成的標(biāo)準(zhǔn)不確定度：A. 總誤差隨機(jī)數(shù)平均值B. 各誤差隨機(jī)數(shù)的殘差C. 按照Bessel公式估計標(biāo)準(zhǔn)不確定度 實驗流程圖：圖6實驗1本實驗中隨機(jī)數(shù)種子為28。以下為N分別為100000點(diǎn)和500000點(diǎn)兩種情況下，M分別等于N/10、N/5、N/2、N、2N、5N六種情況下的模擬圖像。1）程序tic

10、;clear;clc;close all;%設(shè)定參數(shù)值%隨機(jī)信號點(diǎn)數(shù)N，均值為1，標(biāo)準(zhǔn)差u1,u2%N=100000;M=N/10;x=0:1:M;x_=1:M;u1=0.005;u2=0.007;%產(chǎn)生兩個在(0,1)上服從均勻分布的,種子為28,每一次都相同的隨機(jī)數(shù)X1和X2%rand('state',28);X1=rand(1,N);X2=rand(1,N);%按照Box-Mueller變換方法產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布Y1和Y2%Y1=sqrt(-2*log(X1).*cos(2*pi*X2);Y2=sqrt(-2*log(X1).*sin(2*pi*X2);% 為做直方圖先定義

11、好X軸的坐標(biāo)數(shù)據(jù)%delta1=u1*Y1;delta2=u2*Y2;delta=delta1+delta2;d_delta=(max(delta)-min(delta)/(M-1); %d_delta為誤差分布的間距delta_n=min(delta):d_delta:max(delta); %delta_n為誤差分布序列%作圖%高斯隨機(jī)信號%figure(1),axis(0,N,-max(5*Y1),max(5*Y1)plot(Y1);grid on;title('學(xué)號：13S101028 姓名：李明');figure(2),axis(0,N,-max(5*Y2),max(

12、5*Y2)plot(Y2);grid on;title('學(xué)號：13S101028 姓名：李明');% hold on% plot(x,0,'k');grid on;% plot(x,1,'r-');grid on;% plot(x,-1,'r-');grid on;% hold on%變換為任意均值和方差的正態(tài)分布%Z1=Sigma*Y1+Mu;%作圖%高斯隨機(jī)信號% subplot(2,2,2)% axis(0,N,-6,6)% plot(Z1);grid on;% hold on% plot(x,Mu,'k'

13、);% plot(x,Mu+Sigma,'r-');grid on;% plot(x,Mu-Sigma,'r-');grid on;% hold on%正態(tài)分布誤差1幅度直方圖%figure(3)axis(-1,1,0,N)hist(delta1,M);grid on;title('學(xué)號：13S101028 姓名：李明');%正態(tài)分布誤差2幅度直方圖%figure(4)axis(-1,1,0,N)hist(delta2,M);grid on;title('學(xué)號：13S101028 姓名：李明');%合成誤差幅度直方圖%figure

14、(5)axis(-1,1,0,N)H=hist(delta,M);hist(delta,M);grid on;title('學(xué)號：13S101028 姓名：李明');%畫包絡(luò)線%figure(6)HH=envelope(x_,H);plot(delta_n,HH,'b:');grid on;title('學(xué)號：13S101028 姓名：李明');hold on;%計算直方圖的面積%S=sum(d_delta*abs(H);% 計算直方圖的面積%s_1表示正向直方圖的每一個單元的面積%s_2表示反向直方圖的每一個單元的面積%s_表示正反向兩兩對稱每

15、一對單元的面積%area表示以中心為對稱軸的累加面積i=1:1:M/2;s_1(i)=d_delta*abs(floor(H(floor(M/2+i);s_2(i)=d_delta*abs(floor(H(floor(M/2-i+1);s_(i)=s_1(i)+s_2(i);area(1)=s_(1);for ii=1:M/2-1area(ii+1)=area(ii)+s_(ii);end% 計算99.73%的直方圖面積for iii=1:M/2; area(iii);if (area(iii)-0.9973*S)>=0; breakendendplot(delta_n(M/2-iii+

16、1),delta_n(M/2+iii),H(M/2-iii+1),H(M/2+iii),'ro');grid on; title('學(xué)號：13S101028 姓名：李明'); delta_n_u=(delta_n(floor(M/2+iii)-delta_n(floor(M/2-iii+1)/6 %理論計算標(biāo)準(zhǔn)不確定度%delta_mean=mean(delta);delta_cancha=delta-delta_mean;s=sqrt(sum(delta_cancha.2)/(N-1) toc;2）程序運(yùn)行結(jié)果圖（1）當(dāng)N=100000，M=N/10時：圖 1

17、圖 2圖 3圖 4圖 5圖 6計算結(jié)果：s=0.0086，delta_n_u= 0.0085。（2）當(dāng)更改N與M不同的倍數(shù)關(guān)系時，可得到不同的計算結(jié)果，如下各圖所示：圖7 N=100000, M=N/5; s=0.0086, delta_n_u= 0.0086圖8 N=100000, M=N/2; s=0.0086, delta_n_u= 0.0086圖9 N=100000, M=N; s=0.0086, delta_n_u= 0.0086圖10 N=100000, M=2N; s=0.0086, delta_n_u= 0.0086圖11 N=100000, M=5N; s=0.0086; d

18、elta_n_u= 0.0086（3）當(dāng)N=500000時，計算結(jié)果如下所示：圖12圖13圖14圖15圖16圖17計算結(jié)果：s=0.0086，delta_n_u= 0.0088。（4）當(dāng)更改N與M不同的倍數(shù)關(guān)系時，可得到不同的計算結(jié)果，如下各圖所示：圖18 N=500000, M=N/5; s=0.0086, delta_n_u= 0.0088圖19 N=500000, M=N/2; s=0.0086, delta_n_u= 0.0088圖20 N=500000, M=N; s=0.0086, delta_n_u= 0.0088圖21 N=500000,M=2N; s=0.0086, delt

19、a_n_u= 0.0088圖22 N=500000, M=5N; s=0.0086, delta_n_u=0.0088表1 N=100000時，N與M成不同倍數(shù)k時，直方圖計算結(jié)果與理論計算結(jié)果的差異k=N/M105211/21/5s0.00860.00860.00860.00860.00860.0086delta_n_u0.00850.00860.00860.00860.00860.0086|delta_n_us|0.000100000表2 N=500000時，N與M成不同倍數(shù)k時，直方圖計算結(jié)果與理論計算結(jié)果的差異k=N/M105211/21/5s0.00860.00860.00860.0

20、0860.00860.0086delta_n_u0.00880.00880.00880.00880.00880.0088|delta_n_us|0.00020.00020.00020.00020.00020.00023）實驗需要討論的問題（1）N(采樣點(diǎn)數(shù))，M(劃分的區(qū)間數(shù))與直方圖的關(guān)系(平滑，Y軸的高度)。根據(jù)以上各圖分析知：當(dāng)N固定的情況下，隨著M值得增大直方圖的平滑性變差，Y軸高度下降。其中，M<N時，Bessel公式計算的標(biāo)準(zhǔn)不確定度與99.73%直方圖面積的擴(kuò)展不確定度兩者之間的誤差隨著M的增大而逐漸減小。當(dāng)N改變時，即當(dāng)N增大時可使直方圖更為精細(xì)，且此時不改變直方圖包絡(luò)的

21、基本形狀。 （2） Bessel公式計算的標(biāo)準(zhǔn)不確定度與99.73%直方圖面積的擴(kuò)展不確定度兩者之間會存在誤差，這個誤差與哪些因素有關(guān)(N,M,N>=M)此誤差的大小和M、N的相對大小值有關(guān)。當(dāng)N>=M時，由于對離散的誤差值統(tǒng)計運(yùn)算存在舍入誤差導(dǎo)致誤差，此誤差隨著M的增大可消除此項舍入誤差。當(dāng)M>N時，增大M值，誤差值穩(wěn)定，但不能改善誤差值。實驗2自適應(yīng)MCM法在執(zhí)行自適應(yīng)蒙特卡洛方法的基本過程中，蒙特卡洛試驗次數(shù)不斷增加，直至所需要的各種結(jié)果達(dá)到統(tǒng)計意義上的穩(wěn)定。如果某結(jié)果的兩倍標(biāo)準(zhǔn)偏差小于標(biāo)準(zhǔn)不確定度的數(shù)值容差時，則認(rèn)定該數(shù)值結(jié)果穩(wěn)定。(1). 基于前一個實驗，構(gòu)建衡量M

22、onte Carlo法和GUM法計算得到標(biāo)準(zhǔn)不確定度差值的函數(shù)。(2). 將隨機(jī)數(shù)個數(shù)N，分割區(qū)間數(shù)M分別作為該函數(shù)的自變量，定義自變量的取值范圍，從而獲得相應(yīng)的函數(shù)值。(3). 分別進(jìn)行三維網(wǎng)格作圖和三維曲線作圖，通過觀察曲線獲得最佳的N，M組合。1）程序tic;warning off;a,b=meshgrid(logspace(1,6);for j=1:max(size(a); for jj=1:max(size(b); Result1(j,jj)=shiyan(a(j),b(jj); endendfigure(1),surfc(a,b,Result1);title('學(xué)號：13S

23、101028 姓名：李明');c=logspace(1,6);d=logspace(1,6);for jjj=1:max(size(c); Result2(jjj)=shiyan(c(jjj),d(jjj);endfigure(2),plot3(c,d,Result2);grid on;title('學(xué)號：13S101028 姓名：李明');toc;其中shiyan（）子程序為：function y=shiyan(N,M)%clear;clc;close all;%bdclose all;%設(shè)定參數(shù)值%隨機(jī)信號點(diǎn)數(shù)N，均值為1，標(biāo)準(zhǔn)差u1,u2%N=105;Mu=1;%

24、M=N/10;x=0:1:floor(M);x_=1:floor(M);u1=0.005;u2=0.007;%產(chǎn)生兩個在(0,1)上服從均勻分布的,種子為28,每一次都相同的隨機(jī)數(shù)X1和X2%rand('state',28); X1=rand(1,floor(N);X2=rand(1,floor(N);%按照Box-Mueller變換方法產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布Y1和Y2%Y1=sqrt(-2*log(X1).*cos(2*pi*X2);Y2=sqrt(-2*log(X1).*sin(2*pi*X2);% 為做直方圖先定義好X軸的坐標(biāo)數(shù)據(jù)%delta1=u1*Y1;delta2=u2*

25、Y2;delta=delta1+delta2;d_delta=(max(delta)-min(delta)./(floor(M)-1); %d_delta為誤差分布的間距delta_n=min(delta):d_delta:max(delta); %delta_n為誤差分布序列%作圖%高斯隨機(jī)信號% figure(1),% axis(0,N,-max(5*Y1),max(5*Y1)% plot(Y1);grid on;% % figure(2),% axis(0,N,-max(5*Y2),max(5*Y2)% plot(Y2);grid on;% % hold on% % plot(x,0,&

26、#39;k');grid on;% % plot(x,1,'r-');grid on;% % plot(x,-1,'r-');grid on;% % hold on% % %變換為任意均值和方差的正態(tài)分布% %Z1=Sigma*Y1+Mu;% % %作圖% %高斯隨機(jī)信號% % subplot(2,2,2)% % axis(0,N,-6,6)% % plot(Z1);grid on;% % hold on% % plot(x,Mu,'k');% % plot(x,Mu+Sigma,'r-');grid on;% % plo

27、t(x,Mu-Sigma,'r-');grid on;% % hold on% %正態(tài)分布誤差1幅度直方圖% figure(3)% axis(-1,1,0,N)% hist(delta1,M);grid on;% %正態(tài)分布誤差2幅度直方圖% figure(4)% axis(-1,1,0,N)% hist(delta2,M);grid on;% %合成誤差幅度直方圖% figure(5)% axis(-1,1,0,N) H=hist(delta,floor(M);% hist(delta,M);grid on;% %畫包絡(luò)線% figure(6) % HH=envelope(x

28、_,H);% plot(delta_n,HH,'b:');grid on;% hold on;%計算直方圖的面積% S=sum(d_delta.*abs(H);% 計算直方圖的面積%s_1表示正向直方圖的每一個單元的面積%s_2表示反向直方圖的每一個單元的面積%s_表示正反向兩兩對稱每一對單元的面積%area表示以中心為對稱軸的累加面積i=1:1:floor(M./2);s_1(i)=d_delta.*abs(floor(H(floor(M./2+i);s_2(i)=d_delta.*abs(floor(H(floor(M./2-i+1);s_(i)=s_1(i)+s_2(i)

29、;area(1)=s_(1);for ii=1:floor(M./2)-1area(ii+1)=area(ii)+s_(ii);end% 計算99.73%的直方圖面積for iii=1:floor(M./2); area(iii);if (area(iii)-0.9973*S)>=0; breakendend%plot(delta_n(M/2-iii+1),delta_n(M/2+iii),H(M/2-iii),H(M/2+iii),'ro');grid on; delta_n_u=(delta_n(floor(M./2+iii)-delta_n(floor(M./2-i

30、ii+1)./6;%理論計算標(biāo)準(zhǔn)不確定度%delta_mean=mean(delta);delta_cancha=delta-delta_mean;s=sqrt(sum(delta_cancha.2)./(floor(N)-1);y=abs(delta_n_u-s);2）程序運(yùn)行結(jié)果圖23 圖24 3）實驗需要討論的問題：如何根據(jù)三維網(wǎng)格曲線和三維曲線獲得最佳的N，M組合。根據(jù)shiyan（）子函數(shù)知：程序返回值為y=abs(delta_n_u-s)；顯然，當(dāng)y=0時即可獲得N，M的最佳組合，即三維網(wǎng)格曲線和三維曲線的Z坐標(biāo)為0時的N，M。實驗3基于最短包含區(qū)間的MCM法如果PDF不對稱，可采

31、用最短包含區(qū)間。確定，使得，可得最短包含區(qū)間。(1). 先確定q(P=0.9973,M=1000000,q=PM=10000)(2). 重新計算包絡(luò)線下的面積(不是對稱的時候)(3). 根據(jù)算法：，計算r*(4). r*對應(yīng)的區(qū)間長度極為最短包含區(qū)間1）程序%x為均勻分布，正態(tài)分布，反正弦分布%y=sin(x)為何種分布tic;clear;clc;close all;%設(shè)定參數(shù)值%隨機(jī)信號點(diǎn)數(shù)N，均值1，標(biāo)準(zhǔn)差%N=106;M=N/10;x=0:1:M;x_=1:M;u1=0.005;u2=0.007;%rand('state',28);X1=rand(1,N);X2=rand

32、(1,N);Y1=sin(X1);Y2=sqrt(-2*log(X1).*cos(2*pi*X2);delta1=u1*Y1;delta2=u2*Y2;delta=delta1+delta2;d_delta=(max(delta)-min(delta)/(M-1); %d_delta為誤差分布的間距delta_n=min(delta):d_delta:max(delta); %delta_n為誤差分布序列% 畫直方圖figure(1)axis(-1,1,0,N)hist(Y1,M);grid on;title('學(xué)號：13S101028 姓名：李明');figure(2)axi

33、s(-1,1,0,N)hist(Y2,M);grid on;title('學(xué)號：13S101028 姓名：李明');figure(3);axis(-1,1,0,N)hist(delta,M);grid on;title('學(xué)號：13S101028 姓名：李明');H=hist(delta,M);%畫包絡(luò)線%figure(4)HH=envelope(x_,H);plot(delta_n,HH,'b:');grid on;title('學(xué)號：13S101028 姓名：李明');hold on;%計算直方圖的面積%S=sum(d_delta*abs(H);% 計算直方圖的面積%s_1表示正向直方圖的每一個單元的面積%s_2表示反向直方圖的每一個單元的面積%s_表示正反向兩兩對稱每一對單元的面積%area表示以中心為對稱軸的累加面積i=1:1:M/2;s_1(i)=d_delta*abs(floo