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

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

2、率可以利用二重積分求得：Pp( x, y)dxdy（2）以上積分無(wú)法用解析訪法求解，可以根據(jù)Monte Carlo 方法通過(guò)下式進(jìn)行運(yùn)算：(2* rad )2mp( xk , yk )（3）Pp( x, y)dxdynk 1其中 , (2* rad ) 2 表示與圓域外切的正方形區(qū)域的面積,n 為投點(diǎn)次數(shù) , ( xk , yk ) 表示落在區(qū)域中的點(diǎn)的坐標(biāo)。2 程序設(shè)計(jì)（程序部分可直接粘貼運(yùn)行）：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為標(biāo)準(zhǔn)差r=0.4;n=100000;rad=100;x=unifrnd(-rad,rad,1,n);%在（ -100， 100 ）內(nèi)隨機(jī)均勻取n組 x,y 值，y=unifrnd(-rad,rad,1,n);sum=0;m=0;tic%計(jì)時(shí)for k=1:nif x(1,k)2+y(1,k)2<=rad2%實(shí)現(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 ）式計(jì)算概率3 運(yùn)行結(jié)果及分析：n=10000012345計(jì)算結(jié)果0.69620.69770.69650.69800.6967計(jì)算時(shí)間 （ s）2.1611752.1608332.1965802.1429442.189436n=1000012345計(jì)算結(jié)果0.69060.69970.69530.69720.7025計(jì)算時(shí)間 （ s）0.217610.208250.208710.240710.22048最終結(jié)果為0.7 左右。通過(guò)上表還可以看出，隨即試驗(yàn)的次數(shù)并不能完全的決定最終結(jié)果的準(zhǔn)確性。當(dāng)n=1e5時(shí)，其結(jié)果比起n=1e4 的結(jié)果相對(duì)穩(wěn)定，但是計(jì)算時(shí)

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

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

7、 的積分值作為整體積分值， 就會(huì)產(chǎn)生錯(cuò)誤 （得到的結(jié)果為 ）；類(lèi)似的，若用第二象限2 的積分值作為整體積分值，也會(huì)產(chǎn)生錯(cuò)誤（得到的結(jié)果為）。這都是忽略了、的相互作用造成的。相關(guān)系數(shù)的意義是，當(dāng)=時(shí)，、完全不相關(guān)； r=1 時(shí)，、成線性關(guān)系； r= 1 時(shí)，、成負(fù)線性關(guān)系。用以下程序，繪制不同值時(shí)的、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ū)域的長(zhǎng)軸方向，基本上與走向一致。時(shí)，可以推斷，概率密度值講關(guān)于、軸對(duì)稱(chēng)，即不相關(guān)情況（例） ，可以用某一個(gè)象限的值計(jì)算。5 另一個(gè)方法：利用 MATLAB提供的生成符合二維正態(tài)分布的隨機(jī)二維向量的函數(shù)，可以直接生成n個(gè)符合題目要求分布的點(diǎn)的x、y 坐標(biāo)，直接計(jì)算生成的點(diǎn)位于以圓點(diǎn)為圓心，100m 為半徑的圓內(nèi)的頻率P, 根據(jù)大數(shù)定理（4）式可以得知，當(dāng)n 趨近于無(wú)窮時(shí)，頻率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);%生成服從二維分布的隨機(jī)二維向量fork=1:n%檢驗(yàn)在圓域內(nèi)的點(diǎn)的個(gè)數(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 左右。采用本方法， 實(shí)際上是完全模擬了現(xiàn)實(shí)的投彈過(guò)程，是一種直接符合大數(shù)定理形勢(shì)的MonteCarlo方法。習(xí)題：軋鋼有兩道工序：粗軋和精軋，粗軋鋼坯時(shí)由于各種隨機(jī)因素的影響，得到的鋼材長(zhǎng)度成正態(tài)分布， 其均值可由軋機(jī)調(diào)整， 而方差

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

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

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

13、于 3 的點(diǎn)，其概率密度極小，分布函數(shù)值（概率）接近于 0，若 L 位于該區(qū)域，則浪費(fèi)模式 1 出現(xiàn)的概率基本為 0，失去了討論的價(jià)值， 因此有 m<L+3 。w1+w2 的形勢(shì)可以利用積分性質(zhì)進(jìn)行化簡(jiǎn)：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 點(diǎn)的分布函數(shù)值。將（ 5）代入（ 4），就可以直接采用掃描m 的辦法找到 w1+w2 的最小值點(diǎn)。II. 第二問(wèn)是一個(gè)條件期望的問(wèn)題，在第一問(wèn)的基礎(chǔ)上，利用條件期望公式可以直接得到：E

14、(A| B)E(AB)（ 6）P(B)其中 E ( A | B) 表示在事件B 發(fā)生的條件下A 的期望。對(duì)于本題， （ 6）式的意義是：E(每得到一根規(guī)定長(zhǎng)度鋼材的浪費(fèi))E (每粗軋一根鋼材的浪費(fèi))P(每得到一根規(guī)定長(zhǎng)度的鋼材)因此，可以直接利用第一問(wèn)的結(jié)果除以每得到一根規(guī)定長(zhǎng)度的鋼材的概率即可。2程序設(shè)計(jì)1）第一問(wèn)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)%繪制浪費(fèi)期望值p 同粗扎期望值m的關(guān)系曲線2）第二問(wèn)：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ī)定長(zhǎng)度鋼管的概率，即：L 點(diǎn)（規(guī)定長(zhǎng)度）的分布函數(shù)值v=v+1;endp3運(yùn)行結(jié)果1）第一問(wèn)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：掃描趨勢(shì)曲線可以看到，最終的結(jié)果m 取 2.33 左右時(shí)，可以使得浪費(fèi)的期望值w1+w2 最小 (0.4389m).由曲線可以直觀的看到變化趨勢(shì)，結(jié)合本題的實(shí)際意義很好理解。以m=2.33 作為起始點(diǎn)，當(dāng) m 減小小時(shí)， L( 始終位于 m 左側(cè)，見(jiàn) ”模型建立 ”)靠近 m,則產(chǎn)生第一類(lèi)浪費(fèi)模式（整根報(bào)廢）的鋼材量增加；若 m 增大， L 遠(yuǎn)離 m,第一浪費(fèi)模式的鋼材減少，但第二類(lèi)浪費(fèi)模式（精軋浪費(fèi)） 的數(shù)量增加。由于精軋過(guò)程的浪費(fèi)一定會(huì)少于整根報(bào)廢的情況，所以曲線右端部分相對(duì)于左半部分比較平緩。2）第二問(wèn)mm222.3560.44792.0011.98612.3570.44792.0021.97232.3580.44792.0031.95862.3590.44792.0041