1、Monte Carlo模擬誤差分析課程設計1 實驗目的1.1了解MATLAB軟件的基本功能和使用1.2 學習不確定度的統(tǒng)計模擬分析方法1.3 研究誤差概率密度函數和Bessel公式獲得擴展不確定度的方法和影響因素2實驗原理在誤差分析的過程中，常用的方法是通過測量方程推導出誤差傳遞方程，再通過不確定度的合成公式獲得間接測量量的標準不確定度和擴展不確定度(GUM)。在有些場合下，測量方程較難獲得，在這種情況下研究誤差的特性就需要借助于模擬統(tǒng)計的方式進行計算。Monte Carlo(MCM)法就是較為常用的數學工具，具體原理相見相關資料。此次課程設計中按照實驗要求產生的隨機數可以模擬測量誤差，通過對

2、這些隨機數的概率密度分布函數的面積、包絡線和概率特征點的求取，可以獲得隨機誤差的標準不確定度(MCM)，并與理論上估計標準不確定度的Bessel公式、極差法作(GUM)比較，完成實驗內容。并以此作為基礎，分析GUM法與MCM法的區(qū)別與聯(lián)系，影響MCM法的參數，自適應MCM法和基于最短包含區(qū)間的MCM法。已知兩項誤差分量服從正態(tài)分布，標準不確定度分別為mV， mV，用統(tǒng)計模擬分析法給出兩項誤差和的分布(誤差分布的統(tǒng)計直方圖，合成的標準差，合成的置信概率 P為99.73%的擴展不確定度)。3實驗內容(1). 利用MATLAB軟件生成0，1區(qū)間的均勻分布的隨機數；(2). 給出誤差分量的隨機值：利用

3、MATLAB，由均勻分布隨機數生成標準正態(tài)分布隨機數，誤差分量隨機數可表示為mV；mV(3). 求和的隨機數：誤差和的隨機數；(4). 重復以上步驟，得誤差和的隨機數系列： ；(5). 作誤差和的統(tǒng)計直方圖：以誤差數值為橫坐標，以頻率為縱坐標作圖。作圖區(qū)間應包含所有數據，按數值將區(qū)間等分為組（盡可能大），每組間隔為，記數各區(qū)間的隨機數的數目，以為底，以為高作第()區(qū)間的矩形，最終構成誤差和的分布直方圖，該圖包絡線線即為實驗的誤差分布曲線。(6). 以頻率為界劃定區(qū)間，該區(qū)間半寬即為測量總誤差的置信概率為95%的擴展不確定度。(7). 合成的標準不確定度：4.實驗流程圖：一實驗1本實驗中隨機數種

4、子為014。并使分別取N為100000點和10000點兩種情況下，得到M值分別為5*N, 2*N, N, N/2, N/5, N/10五種情況下的模擬圖像。1實驗1程序tic;clear;clc;close all;%設定參數值%隨機信號點數N，均值為1，標準差u1,u2%N=105;M=N/10;x=0:1:M;x_=1:M;u1=0.005;u2=0.007; %產生兩個在(0,1)上服從均勻分布的,種子為0,每一次都相同的隨機數X1和X2%rand('state',014);X1=rand(1,N);X2=rand(1,N); %按照Box-Mueller變換方法產生標準

5、正態(tài)分布Y1和Y2%Y1=sqrt(-2*log(X1).*cos(2*pi*X2);Y2=sqrt(-2*log(X1).*sin(2*pi*X2);% 為做直方圖先定義好X軸的坐標數據%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(5*Y1)p

6、lot(Y1);grid on;figure(2),axis(0,N,-max(5*Y2),max(5*Y2)plot(Y2);grid on;% 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;%作圖%高斯隨機信號% subplot(2,2,2)% axis(0,N,-6,6)% plot(Z1);grid on;% hold on% plot(x,

7、Mu,'k');% 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;%正態(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,M);hist(delta,M);grid on;%畫包

8、絡線%figure(6)HH=envelope(x_,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: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)

9、+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+1),delta_n(M/2+iii),H(M/2-iii+1),H(M/2+iii),'ro');grid on; delta_n_u=(delta_n(floor(M/2+iii)-delta_n(floor(M/2-iii+1)/6;%理

10、論計算標準不確定度%delta_mean=mean(delta);delta_cancha=delta-delta_mean;s=sqrt(sum(delta_cancha.2)/(N-1);%toc;2. 實驗1程序運行結果圖（1）當M=N/10時Figure 1Figure 2Figure3Figure4Figure 5Figure 6（2）當更改N與M不同的倍數關系時，可得到不同的計算結果，如以下個圖所示：圖1.1 N=105, M=N*5，s=0.0086，detla_n_u=0.0087圖1.2 N=105, M=N*2，s=0.0086，detla_n_u=0.0087圖1.3 N

11、=105, M=N，s=0.0086，detla_n_u=0.0087圖1.4 N=105, M=N/2，s=0.0086，detla_n_u=0.0087圖1.5 N=105, M=N/5，s=0.0086，detla_n_u=0.0086圖1.6 N=105, M=N/10，s=0.0086，detla_n_u=0.0085圖1.7 N=104, M=N*5，s=0.0085，detla_n_u=0.0087圖1.8 N=104, M=N*2，s=0.0085，detla_n_u=0.0087圖1.9 N=104, M=N，s=0.0085，detla_n_u=0.0084圖1.10 N=

12、104, M=N/2，s=0.0085，detla_n_u=0.0082圖1.11 N=104, M=N/5，s=0.0085，detla_n_u=0.0078圖1.12 N=104, M=N/10，s=0.0085，detla_n_u=0.0074表2 N=105時，N與M成不同倍數k時，直方圖計算結果與理論計算結果的差異k=N/M1/51/212510s0.00860.00860.00860.00860.00860.0086delta_n_u0.00870.00870.00870.00870.00860.0085|delta_n_us|0.00010.00010.00010.000100.

13、0001表2 N=104時，N與M成不同倍數k時，直方圖計算結果與理論計算結果的差異k=N/M1/51/212510s0.00850.00850.00850.00850.00850.0085delta_n_u0.00870.00870.00840.00820.00780.0074|delta_n_us|0.00020.00020.00010.00030.00070.00113 實驗需要討論的問題(1). N(采樣點數)，M(劃分的區(qū)間數)與直方圖的關系(平滑，Y軸的高度)。有圖1.11.12可知：當N固定的情況下，隨著M值得增大直方圖的平滑性變差，Y軸高度下降。其中，M<N時，Besse

14、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)定，且不能改善誤差值。二實驗2自適應MCM法在執(zhí)行自適應蒙特卡洛方法的基本過程中，蒙特卡洛

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

16、); Result1(j,jj)=shiyan(a(j),b(jj); endendfigure(1),surfc(a,b,Result1);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;toc;2. 實驗2程序運行結果圖Figure 1 logspace(1,6)Figure 2 logspace(1,6)圖2.1 logspace(1,5)圖2.1 logspace(1,4)3 實

17、驗需要討論的問題如何根據三維網格曲線和三維曲線獲得最佳的N，M組合。本實驗中，N、M的取值，是由logspace(a,b,n)函數確定的。該函數的意義為生成從10的a次方到10的b次方之間按對數等分的n個元素的行向量。當logspace(a,b)時，則Matlab中默認n=50。 由以上實驗可知，當M、N的取值越大，即logspace(a,b)中a與b的取值越大時，其誤差越小，但是由于取值越大，要求的計算量也隨之增加，所以N、M的取值使figure 2中藍色圖線隨著橫坐標的增加越快趨于平滑越好，且越接近原點處其斜率越大，波動越小越好。但是隨著N、M值的增加，其計算量也會隨著增加，這就要求計算機

18、達到一定的配置，所以要在計算機現有配置的基礎上，選擇最佳N、M的值。三實驗3基于最短包含區(qū)間的MCM法如果PDF不對稱，可采用最短包含區(qū)間。確定，使得，可得最短包含區(qū)間。(1). 先確定q(P=0.9973,M=1000000,q=PM=10000)(2). 重新計算包絡線下的面積(不是對稱的時候)(3). 根據算法：，計算r*(4). r*對應的區(qū)間長度極為最短包含區(qū)間1實驗3程序%x為均勻分布，正態(tài)分布，反正弦分布%y=sin(x)為何種分布tic;clear;clc;close all; %設定參數值%隨機信號點數N，均值1，標準差%N=106;M=N/10;x=0:1:M;x_=1:M

19、;u1=0.005;u2=0.007;%rand('state',014);X1=rand(1,N);X2=rand(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)h

20、ist(Y1,M);grid on;figure(2)axis(-1,1,0,N)hist(Y2,M);grid on;figure(3);axis(-1,1,0,N)hist(delta,M);grid on;H=hist(delta,M); %畫包絡線%figure(4)HH=envelope(x_,H);plot(delta_n,HH,'b:');grid on;hold on;%計算直方圖的面積%S=sum(d_delta*abs(H); % 計算直方圖的面積%s_1表示正向直方圖的每一個單元的面積%s_2表示反向直方圖的每一個單元的面積%s_表示正反向兩兩對稱每一對單

21、元的面積%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+1),delta_n(M/2+iii),H(M/2-iii+1),H(M/2+iii),'ro');grid on; delta_n_u1=(delta_n(floor(M/2+iii)