1、數(shù)理統(tǒng)計部分實驗1：威廉向克斯的背景：1872 年，英國學(xué)者威廉向克斯（18121882）花費了整整二十個年頭把的值算到了小數(shù)點后707 位。向克斯死后，人們紀(jì)念他，就在他的墓碑上刻下了他一生心血的結(jié)晶：的707 位小數(shù)。此后半個多世紀(jì)，人們對威廉向克斯的計算結(jié)果深信不疑，以至于在1937 年巴黎博覽會發(fā)現(xiàn)館的天井里，依然顯赫地刻著向克斯的值。又過了若干年，數(shù)學(xué)家法格遜對向克斯的計算結(jié)果產(chǎn)生懷疑，他認(rèn)為在的數(shù)值式中，各數(shù)碼出現(xiàn)的概率都應(yīng)當(dāng)相等。于是，他統(tǒng)計了威廉向克斯的頭608 位小數(shù)中，各數(shù)碼出現(xiàn)的頻率，發(fā)現(xiàn)其不相等法格遜猜想：在的數(shù)值中各數(shù)碼0,1,9出現(xiàn)的可能性大小應(yīng)當(dāng)相等。實驗問題：（

2、1） 驗證法格遜猜想（2） 進(jìn)一步地，能否驗證任一無理數(shù)的小數(shù)中，任一個數(shù)值出現(xiàn)的可能性都是一樣的。求解過程： (1)計算頻率xx=input(請輸入你的驗證無理數(shù),s); nn=input(請輸出入你要計算的無數(shù)理的小數(shù)位數(shù),n其值不超過25萬);xx=sym(xx);tempq=vpa(xx,nn+1);tempq=char(tempq);loca=find(tempq=.);tempq(1:loca)= ;probn=zeros(3,10);probn(1,:)=0:9;for ii=0:9 probn(2,ii+1)=length(find(tempq=num2str(ii);endp

3、robn(3,:)=probn(2,:)./nn例：xx=pi ,nn=100000結(jié)果為：數(shù)字0123456789頻數(shù)0.10000.10140.09910.10030.09970.10030.10030.10030.09980.0990（2）分布的假設(shè)檢驗方法對其進(jìn)行檢驗。假設(shè)：09這十個數(shù)字的出現(xiàn)是等可能的。設(shè)這十個數(shù)字出現(xiàn)的概率為，則。設(shè)這十個數(shù)字出現(xiàn)的頻數(shù)為。設(shè)為樣本容量。我們可以構(gòu)造統(tǒng)計量：由概率統(tǒng)計知識有：如果原假設(shè)成立，則當(dāng)樣本容量時，的極限分布為自由度為9的分布。本例中，十個數(shù)字出現(xiàn)的頻數(shù)分別為9999，10137，9908，10025，9971，10026，10029，10

4、025，9978，9902。程序如下（接上段）%對頻率進(jìn)行檢驗npi=1./10.*ones(1,10).*nn;vi=probn(2,:);chi2=sum(npi-vi).2 ./npi);alpha=input(請輸入檢驗的顯著性水平alpha ,其值為大于零小于1的數(shù));crichi= chi2inv(1-alpha,9);exp1=strcat(在顯著性水平alpha =,num2str(alpha),的條件下,可以認(rèn)為各數(shù)字出現(xiàn)是等可能的);exp2=strcat(在顯著性水平alpha =,num2str(alpha),的條件下,各數(shù)字出現(xiàn)不是等可能的);if chi2 cric

5、hi disp(exp2);else disp(exp1);end例：xx=pi ,nn=10000，顯著性水平=0.05的情況下。結(jié)果為：接受原假設(shè)。能否進(jìn)一步猜想，無理數(shù)的小數(shù)中，任何一個k位數(shù)的出現(xiàn)都是等可能的。能否用一些無理數(shù)來驗證。 實驗2：用回歸找規(guī)律 某個煉鋼廠為了測定某種合金鋼的抗拉強度y(kg/mm2)與鋼中的含碳時x(%)之間的關(guān)系。收集了92組煉鋼中碳含量及強度的數(shù)據(jù)（見表1）。為滿足某個用戶的要求，此種合金鋼的抗拉強度y要大于32(kg/mm2)。若要以90%的把握滿足此要求，應(yīng)把含碳量控制在什么范圍內(nèi)？表1：抗拉強度Y（kg/mm2）與含碳量X的對應(yīng)關(guān)系序號x(%)y

6、(kg/mm2)序號x(%)y(kg/mm2)序號x(%)y(kg/mm2)10.0340.5320.143.5630.1347.520.0441.5330.140.5640.1349.530.0438340.144650.144940.0542.5350.142.5660.144150.0540360.141.5670.144360.0541370.137680.1447.570.0540380.143690.154680.0643390.141.5700.154990.0643.5400.145710.1539.5100.0739.5410.141720.1555110.0743420.1

7、142.5730.1548120.0742.5430.1142740.1648.5130.0842440.1142750.1651140.0842450.1146760.1648150.0842460.1145.5770.1753160.0841.5470.1249780.1850170.0842480.1242.5790.252.5180.0841.5490.1244800.255.5190.0842500.1242810.257200.0942.5510.1243820.2156210.0939.5520.1246.5830.2152.5220.0943.5530.1246.5840.21

8、56230.0939540.1343850.2360240.0942.5550.1346860.2456250.0942560.1343870.2453260.0943570.1344.5880.2453270.0943580.1349.5890.2554.5280.0944.5590.1343900.2661.5290.0943600.1345.5910.2959.5300.0945610.1344.5920.3264310.0945.5620.1346問題分析：這是一個用數(shù)據(jù)來近似規(guī)律的問題。初看起來，大家覺得用回歸的辦法對數(shù)據(jù)進(jìn)行擬合就行了。其實分析時要考慮相關(guān)問題。首先，我分析數(shù)據(jù)反映

9、的是確定性規(guī)律還是不確定性規(guī)律（這里主要指相關(guān)關(guān)系）。如果是確定性規(guī)律，常用插值和擬合的方式來近似。如果是相關(guān)關(guān)系，常用回歸來處理。其次，回歸方法是有條件的。要求試驗誤差項相互獨立同分布、誤差服從正態(tài)分布、回歸關(guān)系和系數(shù)要顯著等等。由于我們對背景了解少，所有的只是大量的實驗數(shù)據(jù)，我們可以認(rèn)為這是一個相關(guān)關(guān)系。來找Y對X的回歸方程?；炯僭O(shè)：（1）抗拉強度y（kg/mm2）主要受含碳量x有影響，其他因素的綜合影響視為隨機(jī)因素。（2）每組數(shù)據(jù)的采集是相互獨立的。（1） 抗拉強度y（kg/mm2）與含碳量x的關(guān)系為相關(guān)關(guān)系。（2） 誤差服從正態(tài)分布.模型的建立（1）作出數(shù)據(jù)的散布圖，找到近似的回歸函

10、數(shù)輸入為：xx=0.03, 0.04, 0.04, 0.05, 0.05, 0.05, 0.05, 0.06, 0.06, 0.07, 0.07, 0.07, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.09, 0.09, 0.09, 0.09, 0.09, 0.09, 0.09, 0.09, 0.09, 0.09, 0.09, 0.09, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.11, 0.11, 0.11, 0.11, 0.11, 0.12, 0.12, 0.12

11、, 0.12, 0.12, 0.12, 0.12, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.14, 0.14, 0.14, 0.14, 0.15, 0.15, 0.15, 0.15, 0.15, 0.16, 0.16, 0.16, 0.17, 0.18, 0.20, 0.20, 0.20, 0.21, 0.21,0.21, 0.23, 0.24, 0.24, 0.24, 0.25, 0.26, 0.29, 0.32 ;yy=40.5, 41.5, 38.0, 42.50, 40.0, 41.00,

12、40.0, 43.00, 43.5, 39.5, 43.0, 42.5, 42.0, 42.0, 42.0, 41.5, 42.0, 41.5, 42.0, 42.5, 39.5, 43.5, 39.0, 42.5, 42.0, 43.0, 43.0, 44.5, 43.0, 45.0, 45.5, 43.5, 40.5, 44.0, 42.5, 41.5, 37.0, 43.0, 41.5, 45.0, 41.0, 42.5, 42.0, 42.0, 46.0, 45.5, 49.0, 42.5, 44.0, 42.0, 43.0, 46.5, 46.5, 43.0, 46.0, 43.0,

13、 44.5, 49.5, 43.0, 45.5, 44.5, 46.0, 47.5, 49.5, 49.0, 41.0, 43.0, 47.5, 46.0, 49.0, 39.5, 55.0, 48.0, 48.5, 51.0, 48.0, 53.0, 50.0, 52.5, 55.5, 57.0, 56.0, 52.5, 56.0, 60.0, 56.0, 53.0, 53.0, 54.5, 61.5, 59.5, 64.0 ;plot(xx,yy,.);輸出為： 圖2：抗拉強度y kg/mm2）對含碳量x的散布圖（2）找近似規(guī)律 從上面的散布圖看，我們可以認(rèn)為回歸關(guān)系是近似于線性的。從而用

14、可用線性回歸來近似。即。 輸入為 aa=ones(length(xx),1); XX=aa,xx; coe=XXyy; 輸出為:coe=coe = 34.7699 87.9259 所以回歸的經(jīng)驗公式為。（3）回歸的擬合優(yōu)度檢驗 對于數(shù)據(jù)，在擬合關(guān)系下，我們可以計算殘差，和的方差的估計值，從而可以計算出標(biāo)準(zhǔn)殘差 如果標(biāo)準(zhǔn)殘差主要分布在區(qū)間-2，2上，而且是完全隨機(jī)的，則我們認(rèn)為擬合是成功的。這一點可以從標(biāo)準(zhǔn)殘差圖來觀察。輸入為：ei=yy-XX*coe;dive=ei*ei/(length(ei)-2);rei=ei/sqrt(dive);plot(yy,rei,o);輸出為：圖3：回歸標(biāo)準(zhǔn)殘差

15、re隨抗拉強度Y的散布圖對散布圖中，我們可以認(rèn)為擬合是合理的，因為絕大多數(shù)數(shù)據(jù)點落在區(qū)間-2，2 上，而且數(shù)據(jù)沒有呈現(xiàn)出其它的規(guī)律性。（4）的方差的齊次性檢驗 對于數(shù)據(jù)，求擬合關(guān)系，我們假設(shè)了殘差方差都。但實際數(shù)據(jù)并不一定是這樣的，從而需要我們對進(jìn)診斷。如果殘差是完全隨機(jī)的，而且它的上下界沒有擴(kuò)大的穩(wěn)勢，則我們認(rèn)為方差都。這一點可以從殘差圖來觀察。輸入為：plot(yy,ei,o) 輸出為：圖4：回歸殘差ei隨抗拉強度Y的散布圖 從上圖以可得到：殘差的方差沒有變大的趨勢，我們可以認(rèn)為隨機(jī)誤差的方差都是。(5)對回歸系數(shù)b的顯著性檢驗假設(shè)：在原假設(shè)成立的前提下，回歸離差平方和的分布服從。殘差平方

16、和有分布服從, 且與相互獨立。 所以可以取檢驗統(tǒng)計量，在顯著性水平的情況下作假設(shè)檢驗。 輸入為：ssr=norm(XX*coe-mean(yy).2;sse=ei*ei;ranF=ssr/(sse/length(ei)-2);alpha=0.05; %顯著性水平alphacrit=finv(1-alpha,1, length(ei)-2); 得到檢驗統(tǒng)計量= 500.87，而顯著性水平的情況下，上側(cè)分位數(shù)為3.9469。 所以拒絕，從而我們認(rèn)為回歸系數(shù)是顯著的。 從而我們的回歸模型是成立的。即 （6）計算控制量x的值 合金鋼的抗拉強度Y的值實際上是一個隨機(jī)變量，它的分布服從從而有其中。即求x使

17、對于。分布的上側(cè)面分位數(shù)為= 1.645.即求x的取值，使輸入為：syms xlxx=(xx-mean(xx)*(xx-mean(xx);dlta =sqrt(1 + 1length(ei) + (x - mean(xx)2/lxx)*sqrt(dive);alpha=0.1; %保證抗拉強度0.9的可能性在32以上tacrit=tinv(1-alpha, length(ei)-2);dlta=dlta.*tacrit;yyxx=1,x*coe;yu = yyxx + dlta;yl = yyxx - dlta;ezplot(yu,-20,20)hold onezplot(yl,-20,20)

18、hold offx0 = fzero(char(yl-32),mean(xx)輸出為：圖5：和的關(guān)系圖，中間的曲線為，上面的為yu,下面的為yl.x0 = 0.3380所以我們可以得到結(jié)論：為滿足用戶的要求，此種合金鋼的抗拉強度y應(yīng)大于32，若要以90%的把握滿足此要求，應(yīng)把含碳量x大于0.3380。進(jìn)一步分析：我們上面的討論是在認(rèn)為每個數(shù)據(jù)是相互獨立的條件下做的擬合，我們的數(shù)據(jù)不一定是相互獨立的。我們能不能檢驗。如果我們的有些檢驗是不成立的，我們有沒有辦法對數(shù)據(jù)或擬合方法進(jìn)行修正，使它也能用回歸的方法來求解模型。 實驗3估計和檢驗 我們常說一個班上的學(xué)生的成績服從正態(tài)分布。但卻鮮有人去舉個例

19、子來驗證。 下面是我們考查我校一個班（160名學(xué)生）的某一門課程的成績?nèi)缦隆?0 95 60 60 60 78 43 63 92 91 95 75 81 80 84 30 71 44 44 88 70 73 83 15 66 71 60 15 83 61 89 87 72 43 63 49 87 80 91 78 60 81 84 97 73 92 81 76 73 69 88 82 70 95 88 62 68 63 81 68 85 70 84 41 80 83 94 72 60 77 83 65 61 79 84 79 61 47 51 90 67 60 44 95 78 63 86 6

20、9 80 65 90 34 37 90 41 65 86 72 82 62 61 41 60 77 91 75 83 67 85 87 79 27 84 84 80 70 74 44 71 61 76 79 72 81 26 51 60 85 91 70 76 71 64 61 90 78 35 87 80 79 82 80 91 32 60 69 71 50 80 44 81 61 60 70 50 76 66 63 44 99實驗問題：1)計算均值、標(biāo)準(zhǔn)差、極差、偏度、峰度，畫出直方圖；2)檢驗分布的正態(tài)性；3)若檢驗符合正態(tài)分布，估計正態(tài)分布的參數(shù)并檢驗參數(shù)。求解過程：（1）理論基礎(chǔ)設(shè)總

21、體為X，是一組容量為n的樣本。1、 表示取值集中的位置的統(tǒng)計量平均值和中位數(shù) 平均值（或均值，數(shù)學(xué)期望）： 中位數(shù)：將數(shù)據(jù)由小到大排序后位于中間位置的那個數(shù)值.2、 表示變異程度的統(tǒng)計量標(biāo)準(zhǔn)差、方差和極差 標(biāo)準(zhǔn)差： 它是各個數(shù)據(jù)與均值偏離程度的度量. 方差：標(biāo)準(zhǔn)差的平方. 極差：樣本中最大值與最小值之差. 3. 表示分布形狀的統(tǒng)計量偏度和峰度偏度： 偏度反映分布的對稱性，g1 0稱為右偏態(tài)，此時數(shù)據(jù)位于均值右邊的比位于左邊的多；g1 0稱為左偏態(tài)，情況相反；而g1接近0，則可認(rèn)為分布是對稱的。峰度： 峰度是分布形狀的另一種度量，正態(tài)分布的峰度為3，若g2比3大很多，表示分布有沉重的尾巴，說明樣

22、本中含有較多遠(yuǎn)離均值的數(shù)據(jù)，因而峰度可用作衡量偏離正態(tài)分布的尺度之一.4. 作頻率直方圖。在直角坐標(biāo)系的橫軸上，標(biāo)出各分點，分別以為底邊，作高為（其中是樣本值落在的頻率，）的矩形，,即得頻率直方圖。（2）計算xx=80 95 60 60 60 78 43 63 92 91 95 75 81 80 84 30 71 44 44 88 70 73 83 15 66 71 60 15 83 61 89 87 72 43 63 49 87 80 91 78 60 81 84 97 73 92 81 76 73 69 88 82 70 95 88 62 68 63 81 68 85 70 84 41 8

23、0 83 94 72 60 77 83 65 61 79 84 79 61 47 51 90 67 60 44 95 78 63 86 69 80 65 90 34 37 90 41 65 86 72 82 62 61 41 60 77 91 75 83 67 85 87 79 27 84 84 80 70 74 44 71 61 76 79 72 81 26 51 60 85 91 70 76 71 64 61 90 78 35 87 80 79 82 80 91 32 60 69 71 50 80 44 81 61 60 70 50 76 66 63 44 99;xxbar=mean(xx); %計算樣本均值s_gu=std(xx); %計算樣本標(biāo)準(zhǔn)差jicha=max(xx)-min(xx); %計算極差g1= skewness(xx); %計算偏度g2= kurtosis(xx); %計算偏度xxarea=min(xx)+jicha.*(1:15)./15; %確定直方度的劃分區(qū)間為15個pinglv, xi=hist(xx,xxarea); %計算各區(qū)間的頻數(shù)和區(qū)間中點xipinglv= pinglv ./ 160; %計算各區(qū)間的頻率bar(xi,pinglv); %繪制頻率直方圖圖1：學(xué)生成績頻率直方圖 結(jié)果