數(shù)學(xué)實驗-10：數(shù)據(jù)的統(tǒng)計與分析

1、實驗 10 ：數(shù)據(jù)的統(tǒng)計與分析習(xí)題 5：炮彈射擊的目標為一圓形區(qū)域，半徑為100m,彈著點以圓心為中心成二位正態(tài)分布，設(shè)在密度函數(shù)式當(dāng)中，x =80m,y =50m, 相關(guān)系數(shù)r=0.4, 求炮彈命中圓形區(qū)域的概率。1 模型建立設(shè)目標中心為坐標原點。 Rad(radium)=100,則圓形區(qū)域可以表示為：: x2y2rad 2著彈點符合二維正態(tài)分布，記其坐標為（x,y ） ,其概率密度為有：22r xy2p(x, y)1exp12)( x 2y 2 )（1）21 21 r 22(1r11 22其中1 = x ，2 =y ，由于中心在原點，所以上式中不含有期望值（=0）。于是炮彈命中圓形區(qū)域的概

2、率可以利用二重積分求得：Pp( x, y)dxdy（2）以上積分無法用解析訪法求解，可以根據(jù)Monte Carlo 方法通過下式進行運算：(2* rad )2mp( xk , yk )（3）Pp( x, y)dxdynk 1其中 , (2* rad ) 2 表示與圓域外切的正方形區(qū)域的面積,n 為投點次數(shù) , ( xk , yk ) 表示落在區(qū)域中的點的坐標。2 程序設(shè)計（程序部分可直接粘貼運行）：1) 構(gòu)造概率密度函數(shù) ,符合（ 1）式function f=prob(s1,s2,r,x,y)f=1/(2*pi*s1*s2*sqrt(1-r2)*exp(-1/(1-r2)/2*(x2/s12-

3、2*r*x*y/s1/s2 +y2/s22);2) 主函數(shù)clear alls1=80;s2=50;%s1,s2為標準差r=0.4;n=100000;rad=100;x=unifrnd(-rad,rad,1,n);%在（ -100， 100 ）內(nèi)隨機均勻取n組 x,y 值，y=unifrnd(-rad,rad,1,n);sum=0;m=0;tic%計時for k=1:nif x(1,k)2+y(1,k)2<=rad2%實現(xiàn) Monte Carlo方法sum=sum+prob(s1,s2,r,x(1,k),y(1,k);m=m+1;%sum為（ 3 ）式右端和式部分endendtocp=(

4、2*rad)2/n*sum%根據(jù)（ 3 ）式計算概率3 運行結(jié)果及分析：n=10000012345計算結(jié)果0.69620.69770.69650.69800.6967計算時間 （ s）2.1611752.1608332.1965802.1429442.189436n=1000012345計算結(jié)果0.69060.69970.69530.69720.7025計算時間 （ s）0.217610.208250.208710.240710.22048最終結(jié)果為0.7 左右。通過上表還可以看出，隨即試驗的次數(shù)并不能完全的決定最終結(jié)果的準確性。當(dāng)n=1e5時，其結(jié)果比起n=1e4 的結(jié)果相對穩(wěn)定，但是計算時

5、間是后者的10 倍，可以推斷若將本方法應(yīng)用于更大規(guī)模的數(shù)據(jù)處理當(dāng)中，必然產(chǎn)生精度和計算速度的矛盾。以上問題是實際上反映了局部抽樣中必然存在的問題，Monte Carlo 算法的理論基礎(chǔ)是Bernoull 大數(shù)定理， 即：n 次獨立重復(fù)試驗中A 發(fā)生的次數(shù) k,與 A 在每次試驗中發(fā)生的概率p 有如下關(guān)系：kp | 1（）lim P|nn而實際中的試驗次數(shù)必然是有限的，所以最終得到的結(jié)果必然會不能完全符合概率值。但是，在多次重復(fù)試驗中，同分布的隨機變量，其總體期望和方差為：_1nE xExiEX（）n i1_1nDXD xDxi（）n2ni 1可以看出，隨著試驗次數(shù)的增加，總體期望并沒有發(fā)生變化

6、，但是方差變小了。這也就是n=10000 時得到的結(jié)果波動性比n=100000 時要強的原因了，試驗次數(shù)越多，試驗結(jié)果偏離實際概率的程度就越小。可知，在更大的情況下，對應(yīng)著一個精度，在該精度要求下，最終結(jié)果可以認為是和概率值完全符合。這里不再繼續(xù)進行次數(shù)更多的實驗。4 一個錯誤的分析：同課本中例不同的是， 本題的、 是相互關(guān)聯(lián)的，即滿足二維正態(tài)分布，而例種的、坐標是相互獨立的，各自滿足一維正態(tài)分布。因此，在平面上，對兩個坐標的取點就必須考慮其相互影響，、的取值不是關(guān)于坐標軸對稱的（實際是關(guān)于原點對稱的），因此在計算積分時區(qū)域位于四個象限內(nèi)的積分也不是完全相等的。若此時仍采用例算法：用第一象限1

7、 的積分值作為整體積分值， 就會產(chǎn)生錯誤 （得到的結(jié)果為 ）；類似的，若用第二象限2 的積分值作為整體積分值，也會產(chǎn)生錯誤（得到的結(jié)果為）。這都是忽略了、的相互作用造成的。相關(guān)系數(shù)的意義是，當(dāng)=時，、完全不相關(guān)； r=1 時，、成線性關(guān)系； r= 1 時，、成負線性關(guān)系。用以下程序，繪制不同值時的、prob（x,y ）三維圖像：r=0.9;%r=0.1/-0.5;x=-100:0.1:100;y=x;X,Y=meshgrid(x,y);Z=prob(s1,s2,r,X,Y);mesh(X,Y,Z)r=0.9r=0.1r=-0.5可以看到 ,三維圖形也基本符合鐘形分布，但是隨著值的變化，概率密度

8、峰值的出現(xiàn)范圍隨之改變， 紅色橢圓區(qū)域的長軸方向，基本上與走向一致。時，可以推斷，概率密度值講關(guān)于、軸對稱，即不相關(guān)情況（例） ，可以用某一個象限的值計算。5 另一個方法：利用 MATLAB提供的生成符合二維正態(tài)分布的隨機二維向量的函數(shù)，可以直接生成n個符合題目要求分布的點的x、y 坐標，直接計算生成的點位于以圓點為圓心，100m 為半徑的圓內(nèi)的頻率P, 根據(jù)大數(shù)定理（4）式可以得知，當(dāng)n 趨近于無窮時，頻率P 就是概率。s1=80;s2=50;%初始條件若干，同前r=0.4;n=100000;rad=100;m=0;mu=0,0;%期望sigma=s12,s1*s2*r;s1*s2*r,s2

9、2;%協(xié)方差矩陣x=mvnrnd(mu,sigma,n);%生成服從二維分布的隨機二維向量fork=1:n%檢驗在圓域內(nèi)的點的個數(shù)ifx(k,1)2+x(k,2)2<=rad2m=m+1;endendP=m/n結(jié)果： P=0.695350000000.699690000000.697370000000.69934000000 ?？梢钥吹阶罱K的結(jié)果仍然在0.7 左右。采用本方法， 實際上是完全模擬了現(xiàn)實的投彈過程，是一種直接符合大數(shù)定理形勢的MonteCarlo方法。習(xí)題：軋鋼有兩道工序：粗軋和精軋，粗軋鋼坯時由于各種隨機因素的影響，得到的鋼材長度成正態(tài)分布， 其均值可由軋機調(diào)整， 而方差

10、是設(shè)備精度決定的， 不能改變；精軋時將軋得到鋼材軋成規(guī)定的長度（可以認為沒有誤差） 。如果粗軋后的鋼材長度達與規(guī)定長度， 精軋時要把多的部分軋掉， 造成浪費； 如果粗軋后的鋼材長度已經(jīng)小雨規(guī)定長度， 則整根報廢， 浪費更嚴重。 問題是已知的鋼材規(guī)定的長度和粗軋后的鋼材長度的均方差，求可以調(diào)整的粗軋時剛才長度的均值，失蹤的浪費最小。從以下兩種目標函數(shù)種選擇一個，在l=2m,=20cm 條件下求均值：（）每粗軋一根剛才的浪費最?。ǎ]得到一根規(guī)定長度的鋼材浪費最小模型建立本題需要建立反映鋼材浪費程度的目標函數(shù)，并使其最小， 但由于涉及到正態(tài)分布的概率密度函數(shù)， 是一個非線優(yōu)化問題。之后的建模并沒有

11、嚴格按照優(yōu)化問題的步驟進行，而是采用了更簡便的數(shù)值掃描的辦法直接找到最小點。I.浪費程度可以直接用浪費的鋼材的長度表征，建立關(guān)于浪費長度的目標函數(shù)：由于粗軋長度的不同會造成兩種浪費模式，因此對兩種模式分別進行研究：1）粗軋的得到的鋼材長度小于規(guī)定長度，全部浪費；2）粗軋的得到的鋼材長度大于規(guī)定長度，大于規(guī)定長度的部分被浪費；對于模式1，浪費量的“期望”值（實際是不同浪費長度用其概率密度加權(quán)后的和）可用以下方法求得：w1lxp( x) dx(1)x 為粗軋得到的鋼材長度，w1 表示模式1 的浪費總長度的估計值， p(x) 粗軋的概率密度函數(shù)；對于浪費模式2，由于是部分浪費，所以浪費長度由x 本身

12、變?yōu)?x-L ，且積分區(qū)域也將變?yōu)長 右方，有下式：w2(xl ) p( x)dx(2)l本題當(dāng)中， x 服從正態(tài)分布，即 :p(x)1exp(xm)2(3)222其中 m 為本題所求的正態(tài)分布期望m，為方差。寫成優(yōu)化問題的一本形式：min w1w2s.t w1lxp(x)dxw2( xl ) p( x)dx（4）lp(x)1exp(xm)2222l ml3最后的不等式約束條件的原因是：m 為正態(tài)分布的期望，若其小于標準長度L , 很明顯將至少有50%的鋼材由于粗軋后小于 L 而被直接浪費，這顯然不是最優(yōu)的方法，所以有 L<m; 考慮到正態(tài)分布的 3 法則，可以認為位于 m 點左方距離大

13、于 3 的點，其概率密度極小，分布函數(shù)值（概率）接近于 0，若 L 位于該區(qū)域，則浪費模式 1 出現(xiàn)的概率基本為 0，失去了討論的價值， 因此有 m<L+3 。w1+w2 的形勢可以利用積分性質(zhì)進行化簡：w1w2lxp(x)dx(xl ) p(x)dxlxp(x)dx lp( x)dx（ 5）lEl *1F (l )ml *1F (l )其中 E 為正態(tài)分布的期望，就是m 的值， F(L) 表示在 x=L 點的分布函數(shù)值。將（ 5）代入（ 4），就可以直接采用掃描m 的辦法找到 w1+w2 的最小值點。II. 第二問是一個條件期望的問題，在第一問的基礎(chǔ)上，利用條件期望公式可以直接得到：E

14、(A| B)E(AB)（ 6）P(B)其中 E ( A | B) 表示在事件B 發(fā)生的條件下A 的期望。對于本題， （ 6）式的意義是：E(每得到一根規(guī)定長度鋼材的浪費)E (每粗軋一根鋼材的浪費)P(每得到一根規(guī)定長度的鋼材)因此，可以直接利用第一問的結(jié)果除以每得到一根規(guī)定長度的鋼材的概率即可。2程序設(shè)計1）第一問clearv=1;forq=2:0.001:2.6%在l<m<l+3*sigma區(qū)間內(nèi)掃描 m值m=q;s=0.2;l=2;Fl=normcdf(l,m,s);%求F(l)p(v,:)=m,m-l*(1-Fl);%(5) 式v=v+1;%用p(v,:)記錄結(jié)果 , 輸出

15、 m,w1+w2endpplot(p(:,1),p(:,2)%繪制浪費期望值p 同粗扎期望值m的關(guān)系曲線2）第二問：clearv=1;forq=2:0.01:2.6;m=q;s=0.2;l=2;Fl=normcdf(l,m,s);p(v,:)=m,(m-l*(1-Fl)/(1-Fl);%式（ 6 ），除以每得到一根規(guī)定長度鋼管的概率，即：L 點（規(guī)定長度）的分布函數(shù)值v=v+1;endp3運行結(jié)果1）第一問mW1+w2mW1+w2212.3350.42892.0010.9972.3360.4292.0020.9942.3370.4292.0030.9912.3380.4292.0040.988

16、2.3390.42912.340.42912.3210.42952.3410.42922.3220.42942.3420.42932.3230.42932.3430.42932.3240.42922.3440.42942.3250.42922.3450.42952.3260.42912.3460.42962.3270.4292.3470.42972.3280.4292.3480.42992.3290.4292.3490.432.330.42892.3310.42892.5970.59982.3320.42892.5980.60082.3330.42892.5990.60172.3340.428

17、92.60.6027表 1：數(shù)據(jù)掃描結(jié)果2w+1w1.110.90.80.70.60.5m=2.3330.422.12.22.32.42.52.62.72.82.93m圖 1：掃描趨勢曲線可以看到，最終的結(jié)果m 取 2.33 左右時，可以使得浪費的期望值w1+w2 最小 (0.4389m).由曲線可以直觀的看到變化趨勢，結(jié)合本題的實際意義很好理解。以m=2.33 作為起始點，當(dāng) m 減小小時， L( 始終位于 m 左側(cè)，見 ”模型建立 ”)靠近 m,則產(chǎn)生第一類浪費模式（整根報廢）的鋼材量增加；若 m 增大， L 遠離 m,第一浪費模式的鋼材減少，但第二類浪費模式（精軋浪費） 的數(shù)量增加。由于精軋過程的浪費一定會少于整根報廢的情況，所以曲線右端部分相對于左半部分比較平緩。2）第二問mm222.3560.44792.0011.98612.3570.44792.0021.97232.3580.44792.0031.95862.3590.44792.0041