數(shù)學(xué)建模算法數(shù)據(jù)的統(tǒng)計描述和分析及咖啡店創(chuàng)業(yè)計劃書

第十章數(shù)據(jù)的統(tǒng)計描述和分析數(shù)理統(tǒng)計研究的對象是受隨機因素影響的數(shù)據(jù)，以下數(shù)理統(tǒng)計就簡稱統(tǒng)計，統(tǒng)計是以概率論為基礎(chǔ)的一門應(yīng)用學(xué)科。數(shù)據(jù)樣本少則幾個，多則成千上萬，人們希望能用少數(shù)幾個包含其最多相關(guān)信息的數(shù)值來體現(xiàn)數(shù)據(jù)樣本總體的規(guī)律。描述性統(tǒng)計就是搜集、整理、加工和分析統(tǒng)計數(shù)據(jù)，使之系統(tǒng)化、條理化，以顯示出數(shù)據(jù)資料的趨勢、特征和數(shù)量關(guān)系。它是統(tǒng)計推斷的基礎(chǔ)，實用性較強，在統(tǒng)計工作中經(jīng)常使用。面對一批數(shù)據(jù)如何進行描述與分析，需要掌握參數(shù)估計和假設(shè)檢驗這兩個數(shù)理統(tǒng)計的最基本方法。我們將用Matlab的統(tǒng)計工具箱(StatisticsToolbox)來實現(xiàn)數(shù)據(jù)的統(tǒng)計描述和分析。§1統(tǒng)計的基本概念1.1總體和樣本總體是人們研究對象的全體，又稱母體，如工廠一天生產(chǎn)的全部產(chǎn)品（按合格品及廢品分類），學(xué)校全體學(xué)生的身高。總體中的每一個基本單位稱為個體，個體的特征用一個變量（如）來表示，如一件產(chǎn)品是合格品記，是廢品記；一個身高170（cm）的學(xué)生記。從總體中隨機產(chǎn)生的若干個個體的集合稱為樣本，或子樣，如件產(chǎn)品，100名學(xué)生的身高，或者一根軸直徑的10次測量。實際上這就是從總體中隨機取得的一批數(shù)據(jù)，不妨記作，稱為樣本容量。簡單地說，統(tǒng)計的任務(wù)是由樣本推斷總體。1.2頻數(shù)表和直方圖一組數(shù)據(jù)（樣本）往往是雜亂無章的，作出它的頻數(shù)表和直方圖，可以看作是對這組數(shù)據(jù)的一個初步整理和直觀描述。將數(shù)據(jù)的取值范圍劃分為若干個區(qū)間，然后統(tǒng)計這組數(shù)據(jù)在每個區(qū)間中出現(xiàn)的次數(shù)，稱為頻數(shù)，由此得到一個頻數(shù)表。以數(shù)據(jù)的取值為橫坐標(biāo)，頻數(shù)為縱坐標(biāo)，畫出一個階梯形的圖，稱為直方圖，或頻數(shù)分布圖。若樣本容量不大，能夠手工作出頻數(shù)表和直方圖，當(dāng)樣本容量較大時則可以借助Matlab這樣的軟件了。讓我們以下面的例子為例，介紹頻數(shù)表和直方圖的作法。例1學(xué)生的身高和體重學(xué)校隨機抽取100名學(xué)生，測量他們的身高和體重，所得數(shù)據(jù)如表身高體重身高體重身高體重身高體重身高體重17275171621666216055155571735816655170631675317360178601737316347165661706016350172571826317159177641695516867168651756717664168501614916963171611786417766170581736717259170621725917758176681756818470169641655216459173741726916952173571736116670163571705616065165581776616963176601776717256165561664917165169621705817264169581677217576164591666316954167541796217663182691867716676169721735916965171711674716865165641685717657170571585116562172531696616958172501625217575174661676316650174641686217059（=1\*romani）數(shù)據(jù)輸入數(shù)據(jù)輸入通常有兩種方法，一種是在交互環(huán)境中直接輸入，如果在統(tǒng)計中數(shù)據(jù)量比較大，這樣作不太方便；另一種辦法是先把數(shù)據(jù)寫入一個純文本數(shù)據(jù)文件data.txt中，格式如例1的表格，有20行、10列，數(shù)據(jù)列之間用空格鍵或Tab鍵分割，該數(shù)據(jù)文件data.txt存放在matlab\work子目錄下，在Matlab中用load命令讀入數(shù)據(jù)，具體作法是：loaddata.txt這樣在內(nèi)存中建立了一個變量data，它是一個包含有個數(shù)據(jù)的矩陣。為了得到我們需要的100個身高和體重各為一列的矩陣，應(yīng)做如下的改變：high=data(:,1:2:9);high=high(:)weight=data(:,2:2:10);weight=weight(:)（=2\*romanii）作頻數(shù)表及直方圖用hist命令實現(xiàn)，其用法是：[N,X]=hist(Y,M)數(shù)組（行、列均可）Y的頻數(shù)表。它將區(qū)間[min(Y),max(Y)]等分為M份（缺省時M設(shè)定為10），N返回M個小區(qū)間的頻數(shù)，X返回M個小區(qū)間的中點。hist(Y,M)數(shù)組Y的直方圖。對于例1的數(shù)據(jù)，編寫程序如下：loaddata.txt;high=data(:,1:2:9);high=high(:);weight=data(:,2:2:10);weight=weight(:);[n1,x1]=hist(high)%下面語句與hist命令等價%n1=[length(find(high<158.1)),...%length(find(high>=158.1&high<161.2)),...%length(find(high>=161.2&high<164.5)),...%length(find(high>=164.5&high<167.6)),...%length(find(high>=167.6&high<170.7)),...%length(find(high>=170.7&high<173.8)),...%length(find(high>=173.8&high<176.9)),...%length(find(high>=176.9&high<180)),...%length(find(high>=180&high<183.1)),...%length(find(high>=183.1))][n2,x2]=hist(weight)subplot(1,2,1)hist(high)subplot(1,2,2)hist(weight)計算結(jié)果略，直方圖如下圖所示：從直方圖上可以看出，身高的分布大致呈中間高、兩端低的鐘形；而體重則看不出什么規(guī)律。要想從數(shù)值上給出更確切的描述，需要進一步研究反映數(shù)據(jù)特征的所謂“統(tǒng)計量”。直方圖所展示的身高的分布形狀可看作正態(tài)分布，當(dāng)然也可以用這組數(shù)據(jù)對分布作假設(shè)檢驗。例2統(tǒng)計下列五行字符串中字符a、g、c、t出現(xiàn)的頻數(shù)1.aggcacggaaaaacgggaataacggaggaggacttggcacggcattacacggagg2.cggaggacaaacgggatggcggtattggaggtggcggactgttcgggga3.gggacggatacggattctggccacggacggaaaggaggacacggcggacataca4.atggataacggaaacaaaccagacaaacttcggtagaaatacagaagctta5.cggctggcggacaacggactggcggattccaaaaacggaggaggcggacggaggc解把上述五行復(fù)制到一個純文本數(shù)據(jù)文件shuju.txt中，放在matlab\work子目錄下，編寫如下程序：clcfid1=fopen('shuju.txt','r');i=1;while(~feof(fid1))data=fgetl(fid1);a=length(find(data==97));b=length(find(data==99));c=length(find(data==103));d=length(find(data==116));e=length(find(data>=97&data<=122));f(i,:)=[abcdea+b+c+d];i=i+1;endfhe=[sum(f(:,1))sum(f(:,2))sum(f(:,3))sum(f(:,4))...sum(f(:,5))sum(f(:,6))]fid2=fopen('pinshu.txt','w');fprintf(fid2,'%8d%8d%8d%8d%8d%8d\n',f');fclose(fid1);fclose(fid2);我們把統(tǒng)計結(jié)果最后寫到一個純文本文件pinshu.txt中，在程序中多引進了幾個變量，是為了檢驗字符串是否只包含a、g、c、t四個字符。1.3統(tǒng)計量假設(shè)有一個容量為的樣本（即一組數(shù)據(jù)），記作，需要對它進行一定的加工，才能提出有用的信息，用作對總體（分布）參數(shù)的估計和檢驗。統(tǒng)計量就是加工出來的、反映樣本數(shù)量特征的函數(shù)，它不含任何未知量。下面我們介紹幾種常用的統(tǒng)計量。（=1\*romani）表示位置的統(tǒng)計量—算術(shù)平均值和中位數(shù)算術(shù)平均值（簡稱均值）描述數(shù)據(jù)取值的平均位置，記作，（1）中位數(shù)是將數(shù)據(jù)由小到大排序后位于中間位置的那個數(shù)值。Matlab中mean(x)返回x的均值，median(x)返回中位數(shù)。（=2\*romanii）表示變異程度的統(tǒng)計量—標(biāo)準(zhǔn)差、方差和極差標(biāo)準(zhǔn)差定義為（2）它是各個數(shù)據(jù)與均值偏離程度的度量，這種偏離不妨稱為變異。方差是標(biāo)準(zhǔn)差的平方。極差是的最大值與最小值之差。Matlab中std(x)返回x的標(biāo)準(zhǔn)差，var(x)返回方差，range(x)返回極差。你可能注意到標(biāo)準(zhǔn)差的定義（2）中，對個的平方求和，卻被除，這是出于無偏估計的要求。若需要改為被除，Matlab可用std(x,1)和var(x,1)來實現(xiàn)。（=3\*romaniii）中心矩、表示分布形狀的統(tǒng)計量—偏度和峰度隨機變量的階中心矩為。隨機變量的偏度和峰度指的是的標(biāo)準(zhǔn)化變量的三階中心矩和四階中心矩：偏度反映分布的對稱性，稱為右偏態(tài)，此時數(shù)據(jù)位于均值右邊的比位于左邊的多；稱為左偏態(tài)，情況相反；而接近0則可認為分布是對稱的。峰度是分布形狀的另一種度量，正態(tài)分布的峰度為3，若比3大得多，表示分布有沉重的尾巴，說明樣本中含有較多遠離均值的數(shù)據(jù)，因而峰度可以用作衡量偏離正態(tài)分布的尺度之一。Matlab中moment(x,order)返回x的order階中心矩，order為中心矩的階數(shù)。skewness(x)返回x的偏度，kurtosis(x)返回峰度。在以上用Matlab計算各個統(tǒng)計量的命令中，若x為矩陣，則作用于x的列，返回一個行向量。對例1給出的學(xué)生身高和體重，用Matlab計算這些統(tǒng)計量，程序如下：clcloaddata.txt;high=data(:,1:2:9);high=high(:);weight=data(:,2:2:10);weight=weight(:);shuju=[highweight];jun_zhi=mean([highweight])zhong_wei_shu=median(shuju)biao_zhun_cha=std(shuju)ji_cha=range(shuju)pian_du=skewness(shuju)feng_du=kurtosis(shuju)統(tǒng)計量中最重要、最常用的是均值和標(biāo)準(zhǔn)差，由于樣本是隨機變量，它們作為樣本的函數(shù)自然也是隨機變量，當(dāng)用它們?nèi)ネ茢嗫傮w時，有多大的可靠性就與統(tǒng)計量的概率分布有關(guān)，因此我們需要知道幾個重要分布的簡單性質(zhì)。1.4統(tǒng)計中幾個重要的概率分布1.4.1分布函數(shù)、密度函數(shù)和分位數(shù)隨機變量的特性完全由它的（概率）分布函數(shù)或（概率）密度函數(shù)來描述。設(shè)有隨機變量，其分布函數(shù)定義為的概率，即。若是連續(xù)型隨機變量，則其密度函數(shù)與的關(guān)系為.分位數(shù)是下面常用的一個概念，其定義為：對于，使某分布函數(shù)的，成為這個分布的分位數(shù)，記作。我們前面畫過的直方圖是頻數(shù)分布圖，頻數(shù)除以樣本容量，稱為頻率，充分大時頻率是概率的近似，因此直方圖可以看作密度函數(shù)圖形的（離散化）近似。1.4.2統(tǒng)計中幾個重要的概率分布（=1\*romani）正態(tài)分布正態(tài)分布隨機變量的密度函數(shù)曲線呈中間高兩邊低、對稱的鐘形，期望（均值），方差，記作，稱均方差或標(biāo)準(zhǔn)差，當(dāng)時稱為標(biāo)準(zhǔn)正態(tài)分布，記作。正態(tài)分布完全由均值和方差決定，它的偏度為0，峰度為3。正態(tài)分布可以說是最常見的（連續(xù)型）概率分布，成批生產(chǎn)時零件的尺寸，射擊中彈著點的位置，儀器反復(fù)量測的結(jié)果，自然界中一種生物的數(shù)量特征等，多數(shù)情況下都服從正態(tài)分布，這不僅是觀察和經(jīng)驗的總結(jié)，而且有著深刻的理論依據(jù)，即在大量相互獨立的、作用差不多大的隨機因素影響下形成的隨機變量，其極限分布為正態(tài)分布。鑒于正態(tài)分布的隨機變量在實際生活中如此地常見，記住下面3個數(shù)字是有用的：68%的數(shù)值落在距均值左右1個標(biāo)準(zhǔn)差的范圍內(nèi)，即；95%的數(shù)值落在距均值左右2個標(biāo)準(zhǔn)差的范圍內(nèi)，即；99.7%的數(shù)值落在距均值左右3個標(biāo)準(zhǔn)差的范圍內(nèi)，即.（=2\*romanii）分布(Chisquare)若為相互獨立的、服從標(biāo)準(zhǔn)正態(tài)分布的隨機變量，則它們的平方和服從分布，記作，稱自由度，它的期望，方差。（=3\*romaniii）分布若，，且相互獨立，則服從分布，記作，稱自由度。分布又稱學(xué)生氏(Student)分布。 分布的密度函數(shù)曲線和曲線形狀相似。理論上時，，實際上當(dāng)時它與就相差無幾了。（=4\*romaniv）分布若，，且相互獨立，則服從分布，記作，稱自由度。1.4.3Matlab統(tǒng)計工具箱(Toolbox\Stats)中的概率分布Matlab統(tǒng)計工具箱中有20種概率分布，這里只對上面所述4種分布列出命令的字符：norm正態(tài)分布；chi2分布；t分布f分布工具箱對每一種分布都提供5類函數(shù)，其命令的字符是：pdf概率密度；cdf分布函數(shù)；inv分布函數(shù)的反函數(shù)；stat均值與方差；rnd隨機數(shù)生成當(dāng)需要一種分布的某一類函數(shù)時，將以上所列的分布命令字符與函數(shù)命令字符接起來，并輸入自變量（可以是標(biāo)量、數(shù)組或矩陣）和參數(shù)就行了，如：p=normpdf(x,mu,sigma)均值mu、標(biāo)準(zhǔn)差sigma的正態(tài)分布在的密度函數(shù)(mu=0，sigma=1時可缺省)。p=tcdf(x,n)分布（自由度n）在x的分布函數(shù)。x=chi2inv(p,n)分布（自由度n）使分布函數(shù)F(x)=p的x(即p分位數(shù))。[m,v]=fstat(n1,n2)分布（自由度n1,n2）的均值m和方差v。幾個分布的密度函數(shù)圖形就可以用這些命令作出，如：x=-6:0.01:6;y=normpdf(x);z=normpdf(x,0,2);plot(x,y,x,z),gtext('N(0,1)'),gtext('N(0,2^2)')分布函數(shù)的反函數(shù)的意義從下例看出：x=chi2inv(0.9,10)x=15.9872如果反過來計算，則P=chi2cdf(15.9872,10)P=0.90001.5正態(tài)總體統(tǒng)計量的分布用樣本來推斷總體，需要知道樣本統(tǒng)計量的分布，而樣本又是一組與總體同分布的隨機變量，所以樣本統(tǒng)計量的分布依賴于總體的分布。當(dāng)總體服從一般的分布時，求某個樣本統(tǒng)計量的分布是很困難的，只有在總體服從正態(tài)分布時，一些重要的樣本統(tǒng)計量（均值、標(biāo)準(zhǔn)差）的分布才有便于使用的結(jié)果。另一方面，現(xiàn)實生活中需要進行統(tǒng)計推斷的總體，多數(shù)可以認為服從（或近似服從）正態(tài)分布，所以統(tǒng)計中人們在正態(tài)總體的假定下研究統(tǒng)計量的分布，是必要的與合理的。設(shè)總體，為一容量的樣本，其均值和標(biāo)準(zhǔn)差由式（1）、（2）確定，則用和構(gòu)造的下面幾個分布在統(tǒng)計中是非常有用的?；颍?）（4）（5）設(shè)有兩個總體和，及由容量分別為，的兩個樣本確定的均值和標(biāo)準(zhǔn)差，則（6）（7）其中，（8）對于（7）式，假定，但它們未知，于是用代替。在下面的統(tǒng)計推斷中我們要反復(fù)用到這些分布。§2參數(shù)估計利用樣本對總體進行統(tǒng)計推斷的一類問題是參數(shù)估計，即假定已知總體的分布，通常是，估計參數(shù)的分布，如。參數(shù)估計分點估計和區(qū)間估計兩種。2.1點估計點估計是用樣本統(tǒng)計量確定總體參數(shù)的一個數(shù)值。評價估計優(yōu)劣的標(biāo)準(zhǔn)有無偏性、最小方差性、有效性等，估計的方法有矩法、極大似然法等。最常用的是對總體均值和方差（或標(biāo)準(zhǔn)差）作點估計。讓我們暫時拋開評價標(biāo)準(zhǔn)，當(dāng)從一個樣本按照式（1）、（2）算出樣本均值和方差后，對和（或）一個自然、合理的點估計顯然是（在字母上加^表示它的估計值），，（9）2.2區(qū)間估計點估計雖然給出了待估參數(shù)的一個數(shù)值，卻沒有告訴我們這個估計值的精度和可信程度。一般地，總體的待估參數(shù)記作（如），由樣本算出的的估計量記作，人們常希望給出一個區(qū)間，使以一定的概率落在此區(qū)間內(nèi)。若有，（10）則稱為的置信區(qū)間，分別稱為置信下限和置信上限，稱為置信概率或置信水平，稱為顯著性水平。給出的置信水平為的置信區(qū)間，稱為的區(qū)間估計。置信區(qū)間越小，估計的精度越高；置信水平越大，估計的可信程度越高。但是這兩個指標(biāo)顯然是矛盾的，通常是在一定的置信水平下使置信區(qū)間盡量小。通俗地說，區(qū)間估計給出了點估計的誤差范圍。2.3參數(shù)估計的Matlab實現(xiàn)Matlab統(tǒng)計工具箱中，有專門計算總體均值、標(biāo)準(zhǔn)差的點估計和區(qū)間估計的函數(shù)。對于正態(tài)總體，命令是[mu,sigma,muci,sigmaci]=normfit(x,alpha)其中x為樣本（數(shù)組或矩陣），alpha為顯著性水平（alpha缺省時設(shè)定為0.05），返回總體均值和標(biāo)準(zhǔn)差的點估計mu和sigma，及總體均值和標(biāo)準(zhǔn)差的區(qū)間估計muci和sigmaci。當(dāng)x為矩陣時返回行向量。Matlab統(tǒng)計工具箱中還提供了一些具有特定分布總體的區(qū)間估計的命令，如expfit，poissfit，gamfit，你可以從這些字頭猜出它們用于哪個分布，具體用法參見幫助系統(tǒng)?！?假設(shè)檢驗統(tǒng)計推斷的另一類重要問題是假設(shè)檢驗問題。在總體的分布函數(shù)完全未知或只知其形式但不知其參數(shù)的情況，為了推斷總體的某些性質(zhì)，提出某些關(guān)于總體的假設(shè)。例如，提出總體服從泊松分布的假設(shè)，又如對于正態(tài)總體提出數(shù)學(xué)期望等于的假設(shè)等。假設(shè)檢驗就是根據(jù)樣本對所提出的假設(shè)做出判斷：是接受還是拒絕。這就是所謂的假設(shè)檢驗問題。3.1單個總體均值的檢驗原假設(shè)（或零假設(shè)）為：。備選假設(shè)有三種可能：；；。3.1.1已知，關(guān)于的檢驗（檢驗）在Matlab中檢驗法由函數(shù)ztest來實現(xiàn)，命令為[h,p,ci]=ztest(x,mu,sigma,alpha,tail)其中輸入?yún)?shù)x是樣本，mu是中的，sigma是總體標(biāo)準(zhǔn)差，alpha是顯著性水平（alpha缺省時設(shè)定為0.05），tail是對備選假設(shè)的選擇：為時用tail=0（可缺省）；為時用tail=1；為時用tail=-1。輸出參數(shù)h=0表示接受，h=1表示拒絕，p表示在假設(shè)下樣本均值出現(xiàn)的概率，p越小越值得懷疑，ci是的置信區(qū)間。例3某車間用一臺包裝機包裝糖果。包得的袋裝糖重是一個隨機變量，它服從正態(tài)分布。當(dāng)機器正常時，其均值為0.5公斤，標(biāo)準(zhǔn)差為0.015公斤。某日開工后為檢驗包裝機是否正常，隨機地抽取它所包裝的糖9袋，稱得凈重為（公斤）：0.4970.5060.5180.5240.4980.5110.5200.5150.512問機器是否正常？解總體已知，，未知。于是提出假設(shè)和。Matlab實現(xiàn)如下：x=[0.4970.5060.5180.5240.498...0.5110.5200.5150.512];[h,p,ci]=ztest(x,0.5,0.015)求得h=1，p=0.0248，說明在0.05的水平下，可拒絕原假設(shè)，即認為這天包裝機工作不正常。3.1.2未知，關(guān)于的檢驗（檢驗）在Matlab中檢驗法由函數(shù)ttest來實現(xiàn)，命令為[h,p,ci]=ttest(x,mu,alpha,tail)例4某種電子元件的壽命(以小時計)服從正態(tài)分布,均未知.現(xiàn)得16只元件的壽命如下:159280101212224379179264222362168250149260485170問是否有理由認為元件的平均壽命大于225(小時)?解按題意需檢驗，取。Matlab實現(xiàn)如下：x=[159280101212224379179264...222362168250149260485170];[h,p,ci]=ttest(x,225,0.05,1)求得h=0，p=0.2570，說明在顯著水平為0.05的情況下，不能拒絕原假設(shè)，認為元件的平均壽命不大于225小時。3.2兩個正態(tài)總體均值差的檢驗（檢驗）還可以用檢驗法檢驗具有相同方差的2個正態(tài)總體均值差的假設(shè)。在Matlab中由函數(shù)ttest2實現(xiàn)，命令為：[h,p,ci]=ttest2(x,y,alpha,tail)與上面的ttest相比，不同處只在于輸入的是兩個樣本x,y（長度不一定相同），而不是一個樣本和它的總體均值；tail的用法與ttest相似，可參看幫助系統(tǒng)。例5在平爐上進行一項試驗以確定改變操作方法的建議是否會增加鋼的得率,試驗是在同一平爐上進行的。每煉一爐鋼時除操作方法外,其它條件都可能做到相同。先用標(biāo)準(zhǔn)方法煉一爐,然后用建議的新方法煉一爐,以后交換進行,各煉了10爐,其得率分別為1°標(biāo)準(zhǔn)方法78.172.476.274.377.478.476.075.676.777.32°新方法79.181.077.379.180.079.179.177.380.282.1設(shè)這兩個樣本相互獨立且分別來自正態(tài)總體和,均未知，問建議的新方法能否提高得率?(取。)解（=1\*romani）需要檢驗假設(shè)，.(=2\*romanii)Matlab實現(xiàn)x=[78.172.476.274.377.478.476.075.676.777.3];y=[79.181.077.379.180.079.179.177.380.282.1];[h,p,ci]=ttest2(x,y,0.05,-1)求得h=1,p=2.2126×10-4。表明在的顯著水平下，可以拒絕原假設(shè)，即認為建議的新操作方法較原方法優(yōu)。3.3分布擬合檢驗在實際問題中，有時不能預(yù)知總體服從什么類型的分布，這時就需要根據(jù)樣本來檢驗關(guān)于分布的假設(shè)。下面介紹檢驗法和專用于檢驗分布是否為正態(tài)的“偏峰、峰度檢驗法”。3.3.1檢驗法：總體的分布函數(shù)為,:總體的分布函數(shù)不是在用下述檢驗法檢驗假設(shè)時，若在假設(shè)下的形式已知，但其參數(shù)值未知，這時需要先用極大似然估計法估計參數(shù)，然后作檢驗。檢驗法的基本思想如下：將隨機試驗可能結(jié)果的全體分為個互不相容的事件。于是在假設(shè)下,我們可以計算（或），。在次試驗中,事件出現(xiàn)的頻率與()往往有差異,但一般來說,若為真,且試驗的次數(shù)又甚多時,則這種差異不應(yīng)該很大?；谶@種想法，皮爾遜使用（11）作為檢驗假設(shè)的統(tǒng)計量。并證明了以下定理。定理若充分大,則當(dāng)為真時(不論中的分布屬什么分布),統(tǒng)計量（11）總是近似地服從自由度為的分布,其中是被估計的參數(shù)的個數(shù)。于是,若在假設(shè)下算得（11）有在顯著性水平下拒絕,否則就接受。注意：在使用檢驗法時，要求樣本容量不小于50,以及每個都不小于5,而且最好是在5以上。否則應(yīng)適當(dāng)?shù)睾喜?，以滿足這個要求。例6下面列出了84個伊特拉斯坎（Etruscan）人男子的頭顱的最大寬度（mm），試檢驗這些數(shù)據(jù)是否來自正態(tài)總體（取。141148132138154142150146155158150140147148144150149145149158143141144144126140144142141140145135147146141136140146142137148154137139143140131143141149148135148152143144141143147146150132142142143153149146149138142149142137134144146147140142140137152145解編寫Matlab程序如下：clcx=[141148132138154142150146155158...150140147148144150149145149158...143141144144126140144142141140...145135147146141136140146142137...148154137139143140131143141149...148135148152143144141143147146...150132142142143153149146149138...142149142137134144146147140142...140137152145];min(x),max(x)%求數(shù)據(jù)中的最小數(shù)和最大數(shù)hist(x,8)%畫直方圖fi=[length(find(x<135)),...length(find(x>=135&x<138)),...length(find(x>=138&x<142)),...length(find(x>=142&x<146)),...length(find(x>=146&x<150)),...length(find(x>=150&x<154)),...length(find(x>=154))]%各區(qū)間上出現(xiàn)的頻數(shù)mu=mean(x),sigma=std(x)%均值和標(biāo)準(zhǔn)差fendian=[135,138,142,146,150,154]%區(qū)間的分點p0=normcdf(fendian,mu,sigma)%分點處分布函數(shù)的值p1=diff(p0)%中間各區(qū)間的概率p=[p0(1),p1,1-p0(6)]%所有區(qū)間的概率chi=(fi-84*p).^2./(84*p)chisum=sum(chi)%皮爾遜統(tǒng)計量的值x_a=chi2inv(0.9,4)%chi2分布的0.9分位數(shù)求得皮爾遜統(tǒng)計量chisum=1.9723，，故在水平0.1下接受，即認為數(shù)據(jù)來自正態(tài)分布總體。3.3.2偏度、峰度檢驗（留作習(xí)題1）3.4其它非參數(shù)檢驗Matlab還提供了一些非參數(shù)方法。3.4.1Wilcoxon秩和檢驗在Matlab中，秩和檢驗由函數(shù)ranksum實現(xiàn)。命令為：[p,h]=ranksum(x,y,alpha)其中x，y可為不等長向量，alpha為給定的顯著水平，它必須為0和1之間的數(shù)量。p返回產(chǎn)生兩獨立樣本的總體是否相同的顯著性概率，h返回假設(shè)檢驗的結(jié)果。如果x和y的總體差別不顯著，則h為零；如果x和y的總體差別顯著，則h為1。如果p接近于零，則可對原假設(shè)質(zhì)疑。例7某商店為了確定向公司或公司購買某種產(chǎn)品，將公司以往各次進貨的次品率進行比較，數(shù)據(jù)如下所示，設(shè)兩樣本獨立。問兩公司的商品的質(zhì)量有無顯著差異。設(shè)兩公司的商品的次品的密度最多只差一個平移，取。：7.03.59.68.16.25.110.44.02.010.5：5.73.24.211.09.76.93.64.85.68.410.15.512.3解分別以記公司的商品次品率總體的均值。所需檢驗的假設(shè)是，.Matlab實現(xiàn)如下：a=[7.03.59.68.16.25.110.44.02.010.5];b=[5.73.24.211.09.76.93.64.85.68.410.15.512.3];[p,h]=ranksum(a,b)求得p=0.8041，h=0，表明兩樣本總體均值相等的概率為0.8041，并不很接近于零，且h=0說明可以接受原假設(shè)，即認為兩個公司的商品的質(zhì)量無明顯差異。3.5中位數(shù)檢驗在假設(shè)檢驗中還有一種檢驗方法為中位數(shù)檢驗,在一般的教學(xué)中不一定介紹,但在實際中也是被廣泛應(yīng)用到的。在Matlab中提供了這種檢驗的函數(shù)。函數(shù)的使用方法簡單，下面只給出函數(shù)介紹。3.5.1signrank函數(shù)signrankWilcoxon符號秩檢驗[p,h]=signrank(x,y,alpha)其中p給出兩個配對樣本x和y的中位數(shù)相等的假設(shè)的顯著性概率。向量x，y的長度必須相同，alpha為給出的顯著性水平，取值為0和1之間的數(shù)。h返回假設(shè)檢驗的結(jié)果。如果這兩個樣本的中位數(shù)之差幾乎為0，則h=0；若有顯著差異，則h=1。3.5.2signtest函數(shù)signtest符號檢驗[p,h]=signtest(x,y,alpha)其中p給出兩個配對樣本x和y的中位數(shù)相等的假設(shè)的顯著性概率。x和y若為向量，二者的長度必須相同；y亦可為標(biāo)量，在此情況下，計算x的中位數(shù)與常數(shù)y之間的差異。alpha和h同上。習(xí)題十1.試用偏度、峰度檢驗法檢驗例6中的數(shù)據(jù)是否來自正態(tài)總體（?。?。2.下面列出的是某工廠隨機選取的20只部件的裝配時間（分）：9.8,10.4,10.6,9.6,9.7,9.9,10.9,11.1,9.6,10.2,10.3,9.6,9.9,11.2,10.6,9.8,10.5,10.1,10.5,9.7。設(shè)裝配時間的總體服從正態(tài)分布，是否可以認為裝配時間的均值顯著地大于10（取）？3.下表分別給出兩個文學(xué)家馬克﹒吐溫(MarkTwain)的八篇小品文及斯諾特格拉斯（Snodgrass）的10篇小品文中由3個字母組成的詞的比例。馬克﹒吐溫0.2250.2620.2170.2400.2300.2290.2350.217斯諾特格拉斯0.2090.2050.1960.2100.2020.2070.2240.2230.2200.201設(shè)兩組數(shù)據(jù)分別來自正態(tài)總體，且兩總體方差相等。兩樣本相互獨立，問兩個作家所寫的小品文中包含由3個字母組成的詞的比例是否有顯著的差異（取）？第十一章方差分析我們已經(jīng)作過兩個總體均值的假設(shè)檢驗，如兩臺機床生產(chǎn)的零件尺寸是否相等，病人和正常人的某個生理指標(biāo)是否一樣。如果把這類問題推廣一下，要檢驗兩個以上總體的均值彼此是否相等，仍然用以前介紹的方法是很難做到的。而你在實際生產(chǎn)和生活中可以舉出許多這樣的問題：從用幾種不同工藝制成的燈泡中，各抽取了若干個測量其壽命，要推斷這幾種工藝制成的燈泡壽命是否有顯著差異；用幾種化肥和幾個小麥品種在若干塊試驗田里種植小麥，要推斷不同的化肥和品種對產(chǎn)量有無顯著影響?？梢钥吹?，為了使生產(chǎn)過程穩(wěn)定，達到優(yōu)質(zhì)、高產(chǎn)，需要對影響產(chǎn)品質(zhì)量的因素進行分析，找出有顯著影響的那些因素，除了從機理方面進行研究外，常常要作許多試驗，對結(jié)果作分析、比較，尋求規(guī)律。用數(shù)理統(tǒng)計分析試驗結(jié)果、鑒別各因素對結(jié)果影響程度的方法稱為方差分析（AnalysisOfVariance），記作ANOVA。人們關(guān)心的試驗結(jié)果稱為指標(biāo)，試驗中需要考察、可以控制的條件稱為因素或因子，因素所處的狀態(tài)稱為水平。上面提到的燈泡壽命問題是單因素試驗，小麥產(chǎn)量問題是雙因素試驗。處理這些試驗結(jié)果的統(tǒng)計方法就稱為單因素方差分析和雙因素方差分析?！?單因素方差分析只考慮一個因素對所關(guān)心的指標(biāo)的影響，取幾個水平，在每個水平上作若干個試驗，試驗過程中除外其它影響指標(biāo)的因素都保持不變（只有隨機因素存在），我們的任務(wù)是從試驗結(jié)果推斷，因素對指標(biāo)有無顯著影響，即當(dāng)取不同水平時指標(biāo)有無顯著差別。取某個水平下的指標(biāo)視為隨機變量，判斷取不同水平時指標(biāo)有無顯著差別，相當(dāng)于檢驗若干總體的均值是否相等。1.1數(shù)學(xué)模型設(shè)取個水平，在水平下總體服從正態(tài)分布，，這里未知，可以互不相同，但假定有相同的方差。又設(shè)在每個水平下都作了次獨立試驗，即從中抽取容量為的樣本，記作，服從，且相互獨立。將這些數(shù)據(jù)列成下表（單因素試驗數(shù)據(jù)表）的形式：…1…2……將第列稱為第組數(shù)據(jù)。判斷的個水平對指標(biāo)有無顯著影響，相當(dāng)于要作以下的假設(shè)檢驗；不全相等由于的取值既受不同水平的影響，又受固定下隨機因素的影響，所以將它分解為，，（1）其中，且相互獨立。記，，（2）是總均值，是水平對指標(biāo)的效應(yīng)。由（1）、（2）模型可表為（3）原假設(shè)為（以后略去備選假設(shè)）（4）1.2統(tǒng)計分析記，（5）是第組數(shù)據(jù)的組平均值，是總平均值?？疾烊w數(shù)據(jù)對的偏差平方和（6）經(jīng)分解可得記（7）（8）則（9）是各組均值對總方差的偏差平方和，稱為組間平方和；是各組內(nèi)的數(shù)據(jù)對均值偏差平方和的總和。反映不同水平間的差異，則表示在同一水平下隨機誤差的大小。對和作進一步分析可得（10）（11）當(dāng)成立時（12）可知若成立，只反映隨機波動，而若不成立，那它就還反映了的不同水平的效應(yīng)。單從數(shù)值上看，當(dāng)成立時，由（10）、（12）對于一次試驗應(yīng)有而當(dāng)不成立時這個比值將遠大于1。當(dāng)成立時，該比值服從自由度，的分布，即（13）為檢驗，給定顯著性水平，記分布的分位數(shù)為，檢驗規(guī)則為時接受，否則拒絕。以上對的分析相當(dāng)于對組間、組內(nèi)等方差的分析，所以這種假設(shè)檢驗方法稱方差分析。1.3方差分析表將試驗數(shù)據(jù)按上述分析、計算的結(jié)果排成下表的形式，稱為單因素方差分析表。方差來源平方和自由度平方均值值概率因素誤差總和最后一列給出的概率相當(dāng)于。方差分析一般用的顯著性水平是：取，拒絕，稱因素的影響（或各水平的差異）非常顯著；取，不拒絕，但取，拒絕，稱因素的影響顯著；取，不拒絕，稱因素?zé)o顯著影響。1.4Matlab實現(xiàn)Matlab統(tǒng)計工具箱中單因素方差分析的命令是anoval，用法為：p=anoval(x)返回值p是一個概率，當(dāng)p時接受，x為的數(shù)據(jù)矩陣（如上面的單因素試驗數(shù)據(jù)表形式），x的每一列是一個水平的數(shù)據(jù)。另外，還給出一個方差表和一個Box圖。例1為考察5名工人的勞動生產(chǎn)率是否相同，記錄了每人4天的產(chǎn)量，并算出其平均值，如下表。你能從這些數(shù)據(jù)推斷出他們的生產(chǎn)率有無顯著差別嗎？工人天1256254250248236224233027728025232802902303052204298295302289252平均產(chǎn)量269.00292.25264.75280.50240.00解編寫程序如下：x=[256254250248236242330277280252280290230305220298295302289252];p=anova1(x)求得，故接受，即5名工人的生產(chǎn)率沒有顯著差異。方差表對應(yīng)于上面的單因素方差分析表的列，是分布的分位數(shù)，可以驗證fcdf(2.262,4,15)=0.8891=1-pBox圖反映了各組數(shù)據(jù)的特征。注：接受，是將5名工人的生產(chǎn)率作為一個整體進行假設(shè)檢驗的結(jié)果，并不表明取其中2個工人的生產(chǎn)率作兩總體的均值檢驗時，也一定接受均值相等的假設(shè)。實際上，讀者可以用ttest2對本題作的檢驗，看看會得到什么結(jié)果。1.5非均衡數(shù)據(jù)的方差分析上面所討論的情況是個樣本的容量即各組數(shù)據(jù)個數(shù)相等，稱為均衡數(shù)據(jù)。若各組數(shù)據(jù)個數(shù)不等，稱非均衡數(shù)據(jù)。非均衡數(shù)據(jù)的方差分析，其數(shù)學(xué)模型和統(tǒng)計分析的思路和方法與上面一樣。anova1也能處理非均衡數(shù)據(jù)，與處理均衡數(shù)據(jù)的區(qū)別僅在于數(shù)據(jù)輸入的不同：p=anova1(x,group)x為數(shù)組，從第1組到第組數(shù)據(jù)依次排列；group為與x同長度的數(shù)組，標(biāo)志x中數(shù)據(jù)的組別（在與x第組數(shù)據(jù)相對應(yīng)的位置處輸入整數(shù)）。例2用4種工藝生產(chǎn)燈泡，從各種工藝制成的燈泡中各抽出了若干個測量其壽命，結(jié)果如下表，試推斷這幾種工藝制成的燈泡壽命是否有顯著差異。工藝序號116201580146015002167016001540155031700164016201610417501720168051800解編寫程序如下：x=[1620158014601500167016001540155017001640162016101750172016801800];x=[x(1:4),x(16),x(5:8),x(9:11),x(12:15)];g=[ones(1,5),2*ones(1,4),3*ones(1,3),4*ones(1,4)];p=anova1(x,g)求得0.01<p=0.0331<0.05，所以幾種工藝制成的燈泡壽命有顯著差異。1.6多重比較在燈泡壽命問題中，為了確定哪幾種工藝制成的燈泡壽命有顯著差異，我們先算出各組數(shù)據(jù)的均值：工藝均值1708163515401585雖然的均值最大，但要判斷它與其它幾種有顯著差異，尚需作兩總體均值的假設(shè)檢驗。用ttest2檢驗的結(jié)果如下：原假設(shè) 0110.14590.02020.0408即與有顯著差異（，但與無顯著差異，要想進一步比較優(yōu)劣，應(yīng)增加試驗數(shù)據(jù)。以上作的幾個兩總體均值的假設(shè)檢驗，是多重比較的一部分。一般多重比較要對所有個總體作兩兩對比，分析相互間的差異。根據(jù)問題的具體情況可以減少對比次數(shù)?！?雙因素方差分析如果要考慮兩個因素對指標(biāo)的影響，各劃分幾個水平，對每一個水平組合作若干次試驗，對所得數(shù)據(jù)進行方差分析，檢驗兩因素是否分別對指標(biāo)有顯著影響，或者還要進一步檢驗兩因素是否對指標(biāo)有顯著的交互影響。2.1數(shù)學(xué)模型設(shè)取個水平，取個水平，在水平組合下總體服從正態(tài)分布，，。又設(shè)在水平組合下作了個試驗，所得結(jié)果記作，服從，，，，且相互獨立。將這些數(shù)據(jù)列成下表的形式：……… …將分解為，，，（14）其中，且相互獨立。記，，，，（15）是總均值，是水平對指標(biāo)的效應(yīng)，是水平對指標(biāo)的效應(yīng)，是水平與對指標(biāo)的交互效應(yīng)。模型表為（16）原假設(shè)為2.2無交互影響的雙因素方差分析如果根據(jù)經(jīng)驗或某種分析能夠事先判定兩因素之間沒有交互影響，每組試驗就不必重復(fù)，即可令，過程大為簡化。2.3Matlab實現(xiàn)統(tǒng)計工具箱中用anova2作雙因素方差分析。命令為p=anova2(x,reps)其中x不同列的數(shù)據(jù)表示單一因素的變化情況，不同行中的數(shù)據(jù)表示另一因素的變化情況。如果每種行—列對（“單元”）有不止一個的觀測值，則用參數(shù)reps來表明每個“單元”多個觀測值的不同標(biāo)號，即reps給出重復(fù)試驗的次數(shù)。下面的矩陣中，列因素有3種水平，行因素有兩種水平，但每組水平有兩組樣本，相應(yīng)地用下標(biāo)來標(biāo)識：例3一火箭使用了4種燃料，3種推進器作射程試驗，每種燃料與每種推進器的組合各發(fā)射火箭2次，得到結(jié)果如下：58.2,52.656.2,41.265.3,60.849.1,42.854.1,50.551.6,48.460.1,58.370.9,73.239.2,40.775.8,71.558.2,51.048.7,41.4試在水平0.05下，檢驗不同燃料（因素）、不同推進器（因素）下的射程是否有顯著差異？交互作用是否顯著？解編寫程序如下：clc,clearx0=[58.2,52.6 56.2,41.2 65.3,60.849.1,42.8 54.1,50.5 51.6,48.460.1,58.3 70.9,73.2 39.2,40.775.8,71.5 58.2,51.0 48.7,41.4];x1=x0(:,1:2:5);x2=x0(:,2:2:6);fori=1:4x(2*i-1,:)=x1(i,:);x(2*i,:)=x2(i,:);endp=anova2(x,2)求得p=0.00350.02600.001，表明各試驗均值相等的概率都為小概率，故可拒絕均值相等假設(shè)。即認為不同燃料（因素）、不同推進器（因素）下的射程有顯著差異，交互作用也是顯著的。習(xí)題十一1.將抗生素注入人體會產(chǎn)生抗生素與血漿蛋白質(zhì)結(jié)合的現(xiàn)象，以致減少了藥效。下表列出5種常用的抗生素注入到牛的體內(nèi)時，抗生素與血漿蛋白質(zhì)結(jié)合的百分比。試在水平下檢驗這些百分比的均值有無顯著的差異。設(shè)各總體服從正態(tài)分布，且方差相同。青霉素四環(huán)素鏈霉素紅霉素氯霉素29.624.328.532.027.332.630.834.85.86.211.08.321.617.418.319.029.232.825.024.22.為分析4種化肥和3個小麥品種對小麥產(chǎn)量的影響，把一塊試驗田等分成36小塊，對種子和化肥的每一種組合種植3小塊田，產(chǎn)量如下表所示（單位公斤），問品種、化肥及二者的交互作用對小麥產(chǎn)量有無顯著影響?；势贩N173,172,173174,176,178177,179,176172,173,174175,173,176178,177,179174,175,173170,171,172177,175,176174,174,175174,173,174169,169,170第十二章回歸分析前面我們講過曲線擬合問題。曲線擬合問題的特點是，根據(jù)得到的若干有關(guān)變量的一組數(shù)據(jù)，尋找因變量與（一個或幾個）自變量之間的一個函數(shù)，使這個函數(shù)對那組數(shù)據(jù)擬合得最好。通常，函數(shù)的形式可以由經(jīng)驗、先驗知識或?qū)?shù)據(jù)的直觀觀察決定，要作的工作是由數(shù)據(jù)用最小二乘法計算函數(shù)中的待定系數(shù)。從計算的角度看，問題似乎已經(jīng)完全解決了，還有進一步研究的必要嗎?從數(shù)理統(tǒng)計的觀點看，這里涉及的都是隨機變量，我們根據(jù)一個樣本計算出的那些系數(shù)，只是它們的一個（點）估計，應(yīng)該對它們作區(qū)間估計或假設(shè)檢驗，如果置信區(qū)間太大，甚至包含了零點，那么系數(shù)的估計值是沒有多大意義的。另外也可以用方差分析方法對模型的誤差進行分析，對擬合的優(yōu)劣給出評價。簡單地說，回歸分析就是對擬合問題作的統(tǒng)計分析。具體地說，回歸分析在一組數(shù)據(jù)的基礎(chǔ)上研究這樣幾個問題：（=1\*romani）建立因變量與自變量之間的回歸模型（經(jīng)驗公式）；（=2\*romanii）對回歸模型的可信度進行檢驗；（=3\*romaniii）判斷每個自變量對的影響是否顯著；（=4\*romaniv）診斷回歸模型是否適合這組數(shù)據(jù)；（=5\*romanv）利用回歸模型對進行預(yù)報或控制?！?多元線性回歸回歸分析中最簡單的形式是，均為標(biāo)量，為回歸系數(shù)，稱一元線性回歸。它的一個自然推廣是為多元變量，形如(1)，或者更一般地（2）其中，是已知函數(shù)。這里對回歸系數(shù)是線性的，稱為多元線性回歸。不難看出，對自變量作變量代換，就可將（2）化為（1）的形式，所以下面以（1）為多元線性回歸的標(biāo)準(zhǔn)型。1.1模型在回歸分析中自變量是影響因變量的主要因素，是人們能控制或能觀察的，而還受到隨機因素的干擾，可以合理地假設(shè)這種干擾服從零均值的正態(tài)分布，于是模型記作（3）其中未知。現(xiàn)得到個獨立觀測數(shù)據(jù)，，由（3）得（4）記，（5），（4）表為（6）1.2參數(shù)估計用最小二乘法估計模型（3）中的參數(shù)。由（4）式這組數(shù)據(jù)的誤差平方和為（7）求使最小，得到的最小二乘估計，記作，可以推出（8）將代回原模型得到的估計值（9）而這組數(shù)據(jù)的擬合值為，擬合誤差稱為殘差，可作為隨機誤差的估計，而（10）為殘差平方和（或剩余平方和），即。1.3統(tǒng)計分析不加證明地給出以下結(jié)果：（=1\*romani）是的線性無偏最小方差估計。指的是是的線性函數(shù)；的期望等于；在的線性無偏估計中，的方差最小。（=2\*romanii）服從正態(tài)分布（11）（=3\*romaniii）對殘差平方和，，且（12）由此得到的無偏估計（13）是剩余方差（殘差的方差），稱為剩余標(biāo)準(zhǔn)差。（=4\*romaniv）對總平方和進行分解，有，（14）其中是由（10）定義的殘差平方和，反映隨機誤差對的影響，稱為回歸平方和，反映自變量對的影響。1.4回歸模型的假設(shè)檢驗因變量與自變量之間是否存在如模型（1）所示的線性關(guān)系是需要檢驗的，顯然，如果所有的都很小，與的線性關(guān)系就不明顯，所以可令原假設(shè)為當(dāng)成立時由分解式（14）定義的滿足(15)在顯著性水平下有分位數(shù)，若，接受；否則，拒絕。注意拒絕只說明與的線性關(guān)系不明顯，可能存在非線性關(guān)系，如平方關(guān)系。還有一些衡量與相關(guān)程度的指標(biāo)，如用回歸平方和在總平方和中的比值定義（16）稱為相關(guān)系數(shù)，越大，與相關(guān)關(guān)系越密切，通常，大于0.8（或0.9）才認為相關(guān)關(guān)系成立。1.5回歸系數(shù)的假設(shè)檢驗和區(qū)間估計當(dāng)上面的被拒絕時，不全為零，但是不排除其中若干個等于零。所以應(yīng)進一步作如下個檢驗：由（11）式，，是對角線上的元素，用代替，由（11）~（13）式，當(dāng)成立時（17）對給定的，若，接受；否則，拒絕。(17)式也可用于對作區(qū)間估計（），在置信水平下，的置信區(qū)間為（18）其中。1.6利用回歸模型進行預(yù)測當(dāng)回歸模型和系數(shù)通過檢驗后，可由給定的預(yù)測，是隨機的，顯然其預(yù)測值（點估計）為（19）給定可以算出的預(yù)測區(qū)間（區(qū)間估計），結(jié)果較復(fù)雜，但當(dāng)較大且接近平均值時，的預(yù)測區(qū)間可簡化為（20）其中是標(biāo)準(zhǔn)正態(tài)分布的分位數(shù)。對的區(qū)間估計方法可用于給出已知數(shù)據(jù)殘差的置信區(qū)間，服從均值為零的正態(tài)分布，所以若某個的置信區(qū)間不包含零點，則認為這個數(shù)據(jù)是異常的，可予以剔除。1.7Matlab實現(xiàn)Matlab統(tǒng)計工具箱用命令regress實現(xiàn)多元線性回歸，用的方法是最小二乘法，用法是：b=regress(Y,X)其中Y,X為按（5）式排列的數(shù)據(jù)，b為回歸系數(shù)估計值。[b,bint,r,rint,stats]=regress(Y,X,alpha)這里Y,X同上，alpha為顯著性水平（缺省時設(shè)定為0.05），b,bint為回歸系數(shù)估計值和它們的置信區(qū)間，r,rint為殘差（向量）及其置信區(qū)間，stats是用于檢驗回歸模型的統(tǒng)計量，有三個數(shù)值，第一個是（見（16）式），第二個是（見（15）式），第3個是與對應(yīng)的概率，拒絕，回歸模型成立。殘差及其置信區(qū)間可以用rcoplot(r,rint)畫圖。例1合金的強度與其中的碳含量有比較密切的關(guān)系，今從生產(chǎn)中收集了一批數(shù)據(jù)如下表：0.100.110.120.130.140.150.160.170.1842.041.545.045.545.047.549.055.050.0試先擬合一個函數(shù)，再用回歸分析對它進行檢驗。解先畫出散點圖：x=0.1:0.01:0.18;y=[42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0];plot(x,y,'+')可知與大致上為線性關(guān)系。設(shè)回歸模型為（21）用regress和rcoplot編程如下：clc,clearx1=[0.1:0.01:0.18]';y=[42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0]';x=[ones(9,1),x1];[b,bint,r,rint,stats]=regress(y,x);b,bint,stats,rcoplot(r,rint)得到b=27.4722137.5000bint=18.685136.259475.7755199.2245stats=0.798527.74690.0012即，，的置信區(qū)間是[18.6851,36.2594]，的置信區(qū)間是[75.7755,199.2245]；，，?？芍Ｐ停?1）成立。觀察命令rcoplot(r,rint)所畫的殘差分布，除第8個數(shù)據(jù)外其余殘差的置信區(qū)間均包含零點，第8個點應(yīng)視為異常點，將其剔除后重新計算，可得b=30.7820109.3985bint=26.280535.283476.9014141.8955stats=0.918867.85340.0002應(yīng)該用修改后的這個結(jié)果。例2某廠生產(chǎn)的一種電器的銷售量與競爭對手的價格和本廠的價格有關(guān)。下表是該商品在10個城市的銷售記錄。元120140190130155175125145180150元10011090150210150250270300250個10210012077469326696585試根據(jù)這些數(shù)據(jù)建立與和的關(guān)系式，對得到的模型和系數(shù)進行檢驗。若某市本廠產(chǎn)品售價160（元），競爭對手售價170（元），預(yù)測商品在該市的銷售量。解分別畫出關(guān)于和關(guān)于的散點圖，可以看出與有較明顯的線性關(guān)系，而與之間的關(guān)系則難以確定，我們將作幾種嘗試，用統(tǒng)計分析決定優(yōu)劣。設(shè)回歸模型為（22）編寫如下程序：x1=[120140190130155175125145180150]';x2=[10011090150210150250270300250]';y=[10210012077469326696585]';x=[ones(10,1),x1,x2];[b,bint,r,rint,stats]=regress(y,x);b,bint,stats得到b=66.51760.4139-0.2698bint=-32.5060165.5411-0.20211.0296-0.4611-0.0785stats=0.65276.57860.0247可以看出結(jié)果不是太好：，取時回歸模型（22）可用，但取則模型不能用；較??；的置信區(qū)間包含了零點。下面將試圖用的二次函數(shù)改進它。1.8多項式回歸如果從數(shù)據(jù)的散點圖上發(fā)現(xiàn)與呈較明顯的二次（或高次）函數(shù)關(guān)系，或者用線性模型（1）的效果不太好，就可以選用多項式回歸。1.8.1一元多項式回歸可用命令polyfit實現(xiàn)。例3將17至29歲的運動員每兩歲一組分為7組，每組兩人測量其旋轉(zhuǎn)定向能力，以考察年齡對這種運動能力的影響?，F(xiàn)得到一組數(shù)據(jù)如下表：年齡17192123252729第一人20.4825.1326.1530.026.120.319.35第二人24.3528.1126.331.426.9225.721.3試建立二者之間的關(guān)系。解數(shù)據(jù)的散點圖明顯地呈現(xiàn)兩端低中間高的形狀，所以應(yīng)擬合一條二次曲線。選用二次模型（23）編寫如下程序：x0=17:2:29;x0=[x0,x0];y0=[20.4825.1326.1530.026.120.319.35...24.3528.1126.331.426.9225.721.3];[p,s]=polyfit(x0,y0,2);p得到p=-0.20218.9782-72.2150即，，。上面的s是一個數(shù)據(jù)結(jié)構(gòu)，用于計算函數(shù)值，如[y,delta]=polyconf(p,x0,s);y得到的擬合值，及預(yù)測值的置信區(qū)間半徑delta。用polytool(x0,y0,2)，可以得到一個如上圖的交互式畫面，在畫面中綠色曲線為擬合曲線，它兩側(cè)的紅線是的置信區(qū)間。你可以用鼠標(biāo)移動圖中的十字線來改變圖下方的值，也可以在窗口內(nèi)輸入，左邊就給出的預(yù)測值及其置信區(qū)間。通過左下方的Export下拉式菜單，可以輸出回歸系數(shù)等。這個命令的用法與下面將介紹的rstool相似。1.8.2多元二項式回歸統(tǒng)計工具箱提供了一個作多元二項式回歸的命令rstool，它也產(chǎn)生一個交互式畫面，并輸出有關(guān)信息，用法是rstool(x,y,model,alpha)其中輸入數(shù)據(jù)x,y分別為矩陣和維向量，alpha為顯著性水平（缺省時設(shè)定為0.05），model由下列4個模型中選擇1個（用字符串輸入，缺省時設(shè)定為線性模型）：linear(線性)：purequadratic(純二次)：interaction（交叉）：quadratic(完全二次)：我們再作一遍例2商品銷售量與價格問題，選擇純二次模型，即（24）編程如下：x1=[120140190130155175125145180150]';x2=[10011090150210150250270300250]';y=[10210012077469326696585]';x=[x1x2];rstool(x,y,'purequadratic')得到一個如圖所示的交互式畫面，左邊是（=151）固定時的曲線及其置信區(qū)間，右邊是（=188）固定時的曲線及其置信區(qū)間。用鼠標(biāo)移動圖中的十字線，或在圖下方窗口內(nèi)輸入，可改變。圖左邊給出的預(yù)測值及其置信區(qū)間，就用這種畫面可以回答例2提出的“若某市本廠產(chǎn)品售價160（元），競爭對手售價170（元），預(yù)測該市的銷售量”問題。圖的左下方有兩個下拉式菜單，一個菜單Export用以向Matlab工作區(qū)傳送數(shù)據(jù)，包括beta(回歸系數(shù))，rmse（剩余標(biāo)準(zhǔn)差），residuals(殘差)。模型（24）的回歸系數(shù)和剩余標(biāo)準(zhǔn)差為beta=-312.58717.2701-1.7337-0.02280.0037rmse=16.6436另一個菜單model用以在上述4個模型中選擇，你可以比較一下它們的剩余標(biāo)準(zhǔn)差，會發(fā)現(xiàn)以模型（24）的rmse=16.6436最小?！?非線性回歸和逐步回歸本節(jié)介紹怎樣用Matlab統(tǒng)計工具箱實現(xiàn)非線性回歸和逐步回歸。2.1非線性回歸非線性回歸是指因變量對回歸系數(shù)（而不是自變量）是非線性的。Matlab統(tǒng)計工具箱中的nlinfit，nlparci，nlpredci，nlintool，不僅給出擬合的回歸系數(shù)，而且可以給出它的置信區(qū)間，及預(yù)測值和置信區(qū)間等。下面通過例題說明這些命令的用法。例4在研究化學(xué)動力學(xué)反應(yīng)過程中，建立了一個反應(yīng)速度和反應(yīng)物含量的數(shù)學(xué)模型，形式為其中是未知的參數(shù)，是三種反應(yīng)物（氫，戊烷，異構(gòu)戊烷）的含量，是反應(yīng)速度。今測得一組數(shù)據(jù)如下表，試由此確定參數(shù)，并給出其置信區(qū)間。的參考值為（0.1，0.05，0.02，1，2）。序號反應(yīng)速度氫戊烷異構(gòu)戊烷18.554703001023.79285801034.8247030012040.024708012052.754708010614.391001901072.54100806584.3547019065913.0010030054108.50100300120110.05100801201211.3228530010133首先，以回歸系數(shù)和自變量為輸入變量，將要擬合的模型寫成函數(shù)文件huaxue.m：functionyhat=huaxue(beta,x);yhat=(beta(4)*x(:,2)-x(:,3)/beta(5))./(1+beta(1)*x(:,1)+...beta(2)*x(:,2)+beta(3)*x(:,3));然后，用nlinfit計算回歸系數(shù)，用nlparci計算回歸系數(shù)的置信區(qū)間，用nlpredci計算預(yù)測值及其置信區(qū)間，編程如下：clc,clearx0=[18.554703001023.79285801034.82470