Monte_Carlo 模擬誤差分析課程設計

1、MonteCarlo模擬誤差分析課程設計1.實驗目的1.1學習并掌握MATLAB軟件的基本功能和使用。1.2學習并掌握基于MonteCarloMethod(MCM)分析的不確定度計算方法。1.3研究GuidetotheexpressionofUncertaintyinMeasurement(GUM)法與MCM法的區(qū)別與聯(lián)系和影響因素，自適應MCM方法，基于最短包含區(qū)間的MCM法。2. MonteCarlo模擬誤差分析的實驗原理在誤差分析的過程中，常用的方法是通過測量方程推導出誤差傳遞方程，再通過不確定度的合成公式獲得間接測量量的標準不確定度和擴展不確定度(GUM)。在有些場合下，測量方程較難獲

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

3、法給出兩項誤差和的分布(誤差分布的統(tǒng)計直方圖，合成的標準差，合成的置信概率P為99.73%的擴展不確定度)。3. MonteCarlo模擬誤差分析的實驗內(nèi)容3.1MCM法與GUM法進行模擬誤差分析和不確定度計算(1) .利用MATLAB軟件生成0,1區(qū)間的均勻分布的隨機數(shù)g；(2) .給出誤差分量的隨機值：利用MATLAB，由均勻分布隨機數(shù)g生成標準正態(tài)分布隨機數(shù)耳，誤差分量隨機數(shù)可表示為8=qu=5耳mV；1111同理得8=qu=7耳mV2222(3) .求和的隨機數(shù)：誤差和的隨機數(shù)8=8+8；12(4) .重復以上步驟，得誤差和的隨機數(shù)系列：8=8+8i=1,2,n；i1i2i(5) .作

4、誤差和的統(tǒng)計直方圖：以誤差數(shù)值為橫坐標，以頻率為縱坐標作圖。作圖區(qū)間應包含所有數(shù)據(jù)，按數(shù)值將區(qū)間等分為m組(m盡可能大)，每組間隔為A，記數(shù)各區(qū)間n的隨機數(shù)的數(shù)目n.，以A為底，以j為高作第j(j=1,2m)區(qū)間的矩形，最終的m組jnA矩形構成誤差和的分布直方圖，該圖包絡線線即為實驗的誤差分布曲線。丈n(6) .以頻率4-=99.73%為界劃定區(qū)間，該區(qū)間半寬即為測量總誤差的置信概率為n99.73%的擴展不確定度。(7) .合成的標準不確定度：A.總誤差隨機數(shù)平均值nii=1B.各誤差隨機數(shù)的殘差v=88iiC.按照Bessel公式估計標準不確定度i實驗流程圖：圖6實驗說明：本實驗中隨機數(shù)種子

5、為31。以下為N分別為100000點和500000點兩種情況下M分別等于N/10、N/5、N/2、N、2N、5N六種情況下的模擬圖像。實驗程序：tic;clear;clc;closeall;bdcloseall;%設定參數(shù)值%隨機信號點數(shù)N,均值Mu,標準差Sigma%N=10人5;Mu=1;Sigma=2;M=N/10;x=0:1:M;x_=1:M;u1=0.005;u2=0.007;%產(chǎn)生兩個在(0,1)上服從均勻分布的,種子為31,每一次都相同的隨機數(shù)XI和X2%rand('state',31);X1=rand(1,N);X2=rand(1,N);%按照Box-Muell

6、er變換方法產(chǎn)生標準正態(tài)分布Y1和Y2%Y1=sqrt(-2*log(X1).*cos(2*pi*X2);Y2=sqrt(-2*log(X1).*sin(2*pi*X2);%為做直方圖先定義好X軸的坐標數(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為誤差分布序列%作圖%高斯隨機信號%figure(1),axis(0,N,-max(5*Y1),max

7、(5*Y1)plot(Y1);gridon;figure(2),axis(0,N,-max(5*Y2),max(5*Y2)plot(Y2);gridon;%holdon%plot(x,0,'k');gridon;%plot(x,1,'r-');gridon;%plot(x,-1,'r-');gridon;%holdon%變換為任意均值和方差的正態(tài)分布%Z1=Sigma*Y1+Mu;%作圖%高斯隨機信號%subplot(2,2,2)%axis(0,N,-6,6)%plot(Z1);gridon;%holdon%plot(x,Mu,'k

8、9;);%plot(x,Mu+Sigma,'r-');gridon;%plot(x,Mu-Sigma,'r-');gridon;%holdon%正態(tài)分布誤差1幅度直方圖%figure(3)axis(-1,1,0,N)hist(delta1,M);gridon;%正態(tài)分布誤差2幅度直方圖%figure(4)axis(-1,1,0,N)hist(delta2,M);gridon;%合成誤差幅度直方圖%figure(5)axis(-1,1,0,N)H=hist(delta,M);hist(delta,M);gridon;%畫包絡線%figure(6)HH=envelo

9、pe(x_,H);plot(delta_n,HH,'b:');gridon;holdon;%計算直方圖的面積%S=sum(d_delta*abs(H);%計算直方圖的面積%s_1表示正向直方圖的每一個單元的面積%s_2表示反向直方圖的每一個單元的面積%s_表示正反向兩兩對稱每一對單元的面積%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);fo

10、rii=1:M/2-1area(ii+1)=area(ii)+s_(ii);end%計算99.73%的直方圖面積foriii=1:M/2;area(iii);if(area(iii)-0.9973*S)>=0;breakendendplot(delta_n(M/2-iii+1),delta_n(M/2+iii),H(M/2-iii),H(M/2+iii),'ro');gridon;delta_n_u=(delta_n(floor(M/2+iii)-delta_n(floor(M/2-iii+1)/6;%理論計算標準不確定度%delta_mean=mean(delta);d

11、elta_cancha=delta-delta_mean;s=sqrt(sum(delta_cancha.A2)/(N-l);%toc;程序運行結果(1)當N=1OOOOO,M=N/1O時:FigurelFigure2Figure3Figure5Figure6計算結果：s=O.OO86,delta_n_u=0.0089。(2)當更改N與M不同的倍數(shù)關系時，可得到不同的計算結果，如下各圖所示:35302520151017-0.04-0.03-0.02-0.0100.010.020.030.040.05圖3.1N=100000,M=N/5;s=0.0086,delta_n_u=0.0090圖3.2

12、N=100000,M=N/2;s=0.0086,deltanu=0.0091圖3.3N=100000,M=N;s=0.0086,delta_n_u=0.0091圖3.4N=100000,M=2N;s=0.0086,delta_n_u=0.0091圖3.5N=100000,M=5N;s=0.0086;deltanu=0.0091圖3.6N=500000,M=N/5;s=0.0086,delta_n_u=0.0086圖3.7N=500000,M=N/2;s=0.0086,delta_n_u=0.0086圖3.8N=500000,M=N;s=0.0086,delta_n_u=0.0086圖3.9N=

13、500000,M=2N;s=0.0086,deltanu=0.008604oLjKOMLT!fWWWVW5n;BI:BI:BI;s!;j888888883SSWSRSSSSyi3Ioo.G02o01o02-03Io-0.140圖3.10N=500000,M=5N;s=0.0086,delta_n_u=0.0086表1N=100000時，N與M成不同倍數(shù)k時，直方圖計算結果與理論計算結果的差異k=N/M105211/21/5s0.00860.00860.00860.00860.00860.0086deltanu0.00890.00900.00910.00910.00910.0091ldeltan

14、u-si0.00030.00040.00050.00050.00050.0005表2N=500000時，N與M成不同倍數(shù)k時，直方圖計算結果與理論計算結果的差異k=N/M105211/21/5s0.00860.00860.00860.00860.00860.0086deltanu0.00860.00860.00860.00860.00860.0086|deltanu-s|000000實驗需要討論的問題：(1)N(采樣點數(shù))，M(劃分的區(qū)間數(shù))與直方圖的關系(平滑，Y軸的高度)。根據(jù)以上各圖分析知：當N固定的情況下，隨著M值得增大直方圖的平滑性變差，Y軸高度下降。其中，M<N時，Besse

15、l公式計算的標準不確定度與99.73%直方圖面積的擴展不確定度兩者之間的誤差隨著M的增大而逐漸減小。當N改變時，即當N增大時可使直方圖更為精細，且此時不改變直方圖包絡的基本形狀。(2)Bessel公式計算的標準不確定度與99.73%直方圖面積的擴展不確定度兩者之間會存在誤差，這個誤差與哪些因素有關(N,M,N>=M)此誤差的大小和M、N的相對大小值有關。當N>=M時，由于對離散的誤差值統(tǒng)計運算存在舍入誤差導致誤差，此誤差隨著M的增大可消除此項舍入誤差。當M>N時，增大M值，誤差值穩(wěn)定，但不能改善誤差值。3.2自適應MCM法在執(zhí)行自適應蒙特卡洛方法的基本過程中，蒙特卡洛試驗次數(shù)

16、不斷增加，直至所需要的各種結果達到統(tǒng)計意義上的穩(wěn)定。如果某結果的兩倍標準偏差小于標準不確定度的數(shù)值容差時，則認定該數(shù)值結果穩(wěn)定。(1) .基于前一個實驗，構建衡量MonteCarlo法和GUM法計算得到標準不確定度差值的函數(shù)。(2) .將隨機數(shù)個數(shù)N,分割區(qū)間數(shù)M分別作為該函數(shù)的自變量，定義自變量的取值范圍，從而獲得相應的函數(shù)值。(3) .分別進行三維網(wǎng)格作圖和三維曲線作圖，通過觀察曲線獲得最佳的N,M組合。實驗示例程序：tic;warningoff;a,b=meshgrid(logspace(1,6);forj=1:max(size(a);forjj=1:max(size(b);Result

17、1(j,jj)=shiyan(a(j),b(jj);endendfigure(1),surfc(a,b,Result1);c=logspace(1,6);d=logspace(1,6);forjjj=1:max(size(c);Result2(jjj)=shiyan(c(jjj),d(jjj);endfigure(2),plot3(c,d,Result2);gridon;toc;其中shiyan()子程序為：functiony=shiyan(N,M)%clear;clc;closeall;%bdcloseall;%設定參數(shù)值%隨機信號點數(shù)N,均值為1,標準差ul,u2%N=10人5;Mu=1;

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

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

20、plot(x,1,'r-');gridon;%plot(x,-1,'r-');gridon;%holdon%變換為任意均值和方差的正態(tài)分布%Z1=Sigma*Y1+Mu;%作圖%高斯隨機信號%subplot(2,2,2)%axis(0,N,-6,6)%plot(Z1);gridon;%holdon%plot(x,Mu,'k');%plot(x,Mu+Sigma,'r-');gridon;%plot(x,Mu-Sigma,'r-');gridon;%holdon%正態(tài)分布誤差1幅度直方圖%figure(3)%axis

21、(-1,1,0,N)%hist(delta1,M);gridon;%正態(tài)分布誤差2幅度直方圖%figure(4)%axis(-1,1,0,N)%hist(delta2,M);gridon;%合成誤差幅度直方圖%figure(5)%axis(-1,1,0,N)H=hist(delta,floor(M);%hist(delta,M);gridon;%畫包絡線%figure(6)%HH=envelope(x_,H);%plot(delta_n,HH,'b:');gridon;%holdon;%計算直方圖的面積%S=sum(d_delta.*abs(H);%計算直方圖的面積%s_1表示

22、正向直方圖的每一個單元的面積%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);area(1)=s_(1);forii=1:floor(M./2)-1area(ii+1)=area(ii)+s_(ii);end%計算99.73%的直方圖面積foriii=1:floor(M./2);ar

23、ea(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');gridon;delta_n_u=(delta_n(floor(M./2+iii)-delta_n(floor(M./2-iii+1)./6;%理論計算標準不確定度%delta_mean=mean(delta);delta_cancha=delta-delta_mean;s=sqrt(sum(delta_cancha.人2)./(floor(N)

24、-l);y=abs(delta_n_u-s);程序運行結果：l23實驗需要討論的問題：如何根據(jù)三維網(wǎng)格曲線和三維曲線獲得最佳的N，M組合。根據(jù)shiyan()子函數(shù)知：程序返回值為y=abs(delta_n_u-s)；顯然，當y=0時即可獲得N,M的最佳組合，即三維網(wǎng)格曲線和三維曲線的Z坐標為0時的N,M。使得3.3基于最短包含區(qū)間的MCM法如果PDF不對稱，可采用最短100p%包含區(qū)間。確定r*，yy5yy，r=1,Mq，可得最短100p%包含區(qū)間y,y(r*+q)(r*)(r+Q)(r)L(”*)(+q)-(1) .先確定q(P=0.9973,M=1000000,q=PM=10000)(2

25、) .重新計算包絡線下的面積(不是對稱的時候)(3) .根據(jù)算法：yy<yy，r=1,Mq，計算r*(r*+q)(r*)(r+q)(r)(4) .r*對應的區(qū)間長度極為最短包含區(qū)間實驗示例程序：%x為均勻分布，正態(tài)分布，反正弦分布%y=sin(x)為何種分布tic;clear;clc;closeall;%設定參數(shù)值%隨機信號點數(shù)N,均值1,標準差%N=10A6;M=N/10;x=0:1:M;x_=1:M;u1=0.005;u2=0.007;%rand('state',31);X1=rand(1,N);X2=rand(1,N);Y1=sin(X1);Y2=sqrt(-2*l

26、og(X1).*cos(2*pi*X2);delta1=u1*Y1;%d_delta為誤差分布的間距%delta_n為誤差分布序列delta2=u2*Y2;delta=delta1+delta2;d_delta=(max(delta)-min(delta)/(M-1);delta_n=min(delta):d_delta:max(delta);%畫直方圖figure(1)axis(-1,1,0,N)hist(Y1,M);gridon;figure(2)axis(-1,1,0,N)hist(Y2,M);gridon;figure(3);axis(-1,1,0,N)hist(delta,M);gr

27、idon;H=hist(delta,M);%畫包絡線%figure(4)HH=envelope(x_,H);plot(delta_n,HH,'b:');gridon;holdon;%計算直方圖的面積%S=sum(d_delta*abs(H);%計算直方圖的面積%s_l表示正向直方圖的每一個單元的面積%s_2表示反向直方圖的每一個單元的面積%s_表示正反向兩兩對稱每一對單元的面積%area表示以中心為對稱軸的累加面積i=l:l:M/2;s_l(i)=d_delta*abs(floor(H(floor(M/2+i);s_2(i)=d_delta*abs(floor(H(floor(M/2-i+l);s_(i)=s_l(i)+s_2(i);area(l)=s_(l);forii=l:M/2-larea(ii+l)=area(ii)+s_(ii);end%計算99.73%的直方圖面積foriii=l:M/