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

實(shí)驗(yàn)10：數(shù)據(jù)的統(tǒng)計(jì)與分析習(xí)題5： 炮彈射擊的目標(biāo)為一圓形區(qū)域，半徑為100m,彈著點(diǎn)以圓心為中心成二位正態(tài)分布，設(shè)在密度函數(shù)式當(dāng)中，=80m,=50m,相關(guān)系數(shù)r=0.4,求炮彈命中圓形區(qū)域的概率。模型建立 設(shè)目標(biāo)中心為坐標(biāo)原點(diǎn)。Rad(radium)=100,則圓形區(qū)域可以表示為：著彈點(diǎn)符合二維正態(tài)分布，記其坐標(biāo)為（x,y）,其概率密度為有： （1）其中=，=，由于中心在原點(diǎn)，所以上式中不含有期望值（=0）。于是炮彈命中圓形區(qū)域的概率可以利用二重積分求得： （2）以上積分無法用解析訪法求解，可以根據(jù)MonteCarlo方法通過下式進(jìn)行運(yùn)算： （3）其中,表示與圓域外切的正方形區(qū)域的面積,n為投點(diǎn)次數(shù),表示落在區(qū)域中的點(diǎn)的坐標(biāo)。程序設(shè)計(jì)（程序部分可直接粘貼運(yùn)行）： 1)構(gòu)造概率密度函數(shù),符合（1）式functionf=prob(s1,s2,r,x,y)f=1/(2*pi*s1*s2*sqrt(1-r^2))*exp(-1/(1-r^2)/2*(x^2/s1^2-2*r*x*y/s1/s2+y^2/s2^2)); 2)主函數(shù)clearalls1=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í)fork=1:n ifx(1,k)^2+y(1,k)^2<=rad^2 %實(shí)現(xiàn)MonteCarlo方法sum=sum+prob(s1,s2,r,x(1,k),y(1,k));m=m+1; %sum為（3）式右端和式部分endendtocp=(2*rad)^2/n*sum %根據(jù)（3）式計(jì)算概率運(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左右。 通過上表還可以看出，隨即試驗(yàn)的次數(shù)并不能完全的決定最終結(jié)果的準(zhǔn)確性。當(dāng)n=1e5時(shí)，其結(jié)果比起n=1e4的結(jié)果相對穩(wěn)定，但是計(jì)算時(shí)間是后者的10倍，可以推斷若將本方法應(yīng)用于更大規(guī)模的數(shù)據(jù)處理當(dāng)中，必然產(chǎn)生精度和計(jì)算速度的矛盾。 以上問題是實(shí)際上反映了局部抽樣中必然存在的問題，MonteCarlo算法的理論基礎(chǔ)是Bernoull大數(shù)定理，即：n次獨(dú)立重復(fù)試驗(yàn)中A發(fā)生的次數(shù)k,與A在每次試驗(yàn)中發(fā)生的概率p有如下關(guān)系： （４）而實(shí)際中的試驗(yàn)次數(shù)必然是有限的，所以最終得到的結(jié)果必然會不能完全符合ｐ概率值。但是，在多次重復(fù)試驗(yàn)中，同分布的隨機(jī)變量，其總體期望和方差為： （５） （６）可以看出，隨著試驗(yàn)次數(shù)的增加，總體期望并沒有發(fā)生變化，但是方差變小了。這也就是n=10000時(shí)得到的結(jié)果波動性比n=100000時(shí)要強(qiáng)的原因了，試驗(yàn)次數(shù)越多，試驗(yàn)結(jié)果偏離實(shí)際概率的程度就越小?？芍冢罡蟮那闆r下，對應(yīng)著一個(gè)精度，在該精度要求下，最終結(jié)果可以認(rèn)為是和概率值完全符合。這里不再繼續(xù)進(jìn)行次數(shù)更多的實(shí)驗(yàn)。一個(gè)錯(cuò)誤的分析： 同課本中例６不同的是，本題的ｘ、ｙ是相互關(guān)聯(lián)的，即滿足二維正態(tài)分布，而例６種的ｘ、ｙ坐標(biāo)是相互獨(dú)立的，各自滿足一維正態(tài)分布。因此，在ｘ－ｙ平面上，對兩個(gè)坐標(biāo)的取點(diǎn)就必須考慮其相互影響，ｘ、ｙ的取值不是關(guān)于坐標(biāo)軸對稱的（實(shí)際是關(guān)于原點(diǎn)對稱的），因此在計(jì)算積分時(shí)區(qū)域位于四個(gè)象限內(nèi)的積分也不是完全相等的。若此時(shí)仍采用例s1=80;s2=50; %初始條件若干，同前r=0.4;n=100000;rad=100;m=0;mu=[0,0]; %期望sigma=[s1^2,s1*s2*r;s1*s2*r,s2^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<=rad^2m=m+1;endendP=m/n結(jié)果：P=0.695350000000.699690000000.697370000000.69934000000。。?？梢钥吹阶罱K的結(jié)果仍然在0.7左右。采用本方法，實(shí)際上是完全模擬了現(xiàn)實(shí)的投彈過程，是一種直接符合大數(shù)定理形勢的MonteCarlo方法。習(xí)題９： 軋鋼有兩道工序：粗軋和精軋，粗軋鋼坯時(shí)由于各種隨機(jī)因素的影響，得到的鋼材長度成正態(tài)分布，其均值可由軋機(jī)調(diào)整，而方差是設(shè)備精度決定的，不能改變；精軋時(shí)將軋得到鋼材軋成規(guī)定的長度（可以認(rèn)為沒有誤差）。如果粗軋后的鋼材長度達(dá)與規(guī)定長度，精軋時(shí)要把多的部分軋掉，造成浪費(fèi)；如果粗軋后的鋼材長度已經(jīng)小雨規(guī)定長度，則整根報(bào)廢，浪費(fèi)更嚴(yán)重。問題是已知的鋼材規(guī)定的長度ｌ和粗軋后的鋼材長度的均方差，求可以調(diào)整的粗軋時(shí)剛才長度的均值ｍ，失蹤的浪費(fèi)最小。從以下兩種目標(biāo)函數(shù)種選擇一個(gè)，在l=2m,=20cm條件下求均值ｍ：（１）每粗軋一根剛才的浪費(fèi)最小（２）沒得到一根規(guī)定長度的鋼材浪費(fèi)最?。保Ｐ徒?本題需要建立反映鋼材浪費(fèi)程度的目標(biāo)函數(shù)，并使其最小，但由于涉及到正態(tài)分布的概率密度函數(shù)，是一個(gè)非線優(yōu)化問題。之后的建模并沒有嚴(yán)格按照優(yōu)化問題的步驟進(jìn)行，而是采用了更簡便的數(shù)值掃描的辦法直接找到最小點(diǎn)。I. 浪費(fèi)程度可以直接用浪費(fèi)的鋼材的長度表征，建立關(guān)于浪費(fèi)長度的目標(biāo)函數(shù)：由于粗軋長度的不同會造成兩種浪費(fèi)模式，因此對兩種模式分別進(jìn)行研究：1）粗軋的得到的鋼材長度小于規(guī)定長度，全部浪費(fèi)；2）粗軋的得到的鋼材長度大于規(guī)定長度，大于規(guī)定長度的部分被浪費(fèi)；對于模式1，浪費(fèi)量的“期望”值（實(shí)際是不同浪費(fèi)長度用其概率密度加權(quán)后的和）可用以下方法求得： (1)x為粗軋得到的鋼材長度，w1表示模式1的浪費(fèi)總長度的估計(jì)值，p(x)粗軋的概率密度函數(shù)；對于浪費(fèi)模式2，由于是部分浪費(fèi)，所以浪費(fèi)長度由x本身變?yōu)閤-L，且積分區(qū)域也將變?yōu)長右方，有下式： (2)本題當(dāng)中，x服從正態(tài)分布，即: (3)其中m為本題所求的正態(tài)分布期望m，為方差。寫成優(yōu)化問題的一本形式： （4）最后的不等式約束條件的原因是： m為正態(tài)分布的期望，若其小于標(biāo)準(zhǔn)長度L,很明顯將至少有50%的鋼材由于粗軋后小于L而被直接浪費(fèi)，這顯然不是最優(yōu)的方法，所以有L<m;考慮到正態(tài)分布的3法則，可以認(rèn)為位于m點(diǎn)左方距離大于3的點(diǎn)，其概率密度極小，分布函數(shù)值（概率）接近于0，若L位于該區(qū)域，則浪費(fèi)模式1出現(xiàn)的概率基本為0，失去了討論的價(jià)值，因此有m<L+3。 w1+w2的形勢可以利用積分性質(zhì)進(jìn)行化簡： （5）其中E為正態(tài)分布的期望，就是m的值，F(xiàn)(L)表示在x=L點(diǎn)的分布函數(shù)值。將（5）代入（4），就可以直接采用掃描m的辦法找到w1+w2的最小值點(diǎn)。 II. 第二問是一個(gè)條件期望的問題，在第一問的基礎(chǔ)上，利用條件期望公式可以直接得到： （6）其中表示在事件B發(fā)生的條件下A的期望。對于本題，（6）式的意義是： 因此，可以直接利用第一問的結(jié)果除以每得到一根規(guī)定長度的鋼材的概率即可。2．程序設(shè)計(jì)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é)果,輸出[m,w1+w2]endpplot(p(:,1),p(:,2)) %繪制浪費(fèi)期望值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點(diǎn)（規(guī)定長度）的分布函數(shù)值 v=v+1;endp3．運(yùn)行結(jié)果 1）第一問mW1+w2mW1+w2212.3350.42892.0010.9972.3360.4292.0020.9942.3370.4292.0030.9912.3380.4292.0040.9882.3390.4291……2.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.4289……2.3310.42892.5970.59982.3320.42892.5980.60082.3330.42892.5990.60172.3340.42892.60.6027表1：數(shù)據(jù)掃描結(jié)果圖1：掃描趨勢曲線 可以看到，最終的結(jié)果m取2.33左右時(shí)，可以使得浪費(fèi)的期望值w1+w2最小(0.4389m). 由曲線可以直觀的看到變化趨勢，結(jié)合本題的實(shí)際意義很好理解。以m=2.33作為起始點(diǎn)，當(dāng)m減小小時(shí)，L(始終位于m左側(cè)，見”模型建立”)靠近m,則產(chǎn)生第一類浪費(fèi)模式（整根報(bào)廢）的鋼材量增加；若m增大，L遠(yuǎn)離m,第一浪費(fèi)模式的鋼材減少，但第二類浪費(fèi)模式（精軋浪費(fèi)）的數(shù)量增加。由于精軋過程的浪費(fèi)一定會少于整根報(bào)廢的情況，所以曲線右端部分相對于左半部分比較平緩。2）第二問mm222.3560.44792.0011.98612.3570.44792.0021.97232.3580.44792.0031.95862.3590.44792.0041.94512.360.448……2.3610.4482.3480.44822.3620