數(shù)理統(tǒng)計(jì)課程實(shí)驗(yàn)報(bào)告

高等數(shù)理統(tǒng)計(jì)課程實(shí)驗(yàn)【摘要】本實(shí)驗(yàn)報(bào)告描述了用最小二成估計(jì)算法解決實(shí)際問題中參數(shù)估計(jì)的過程。包括引言、實(shí)驗(yàn)原理、實(shí)驗(yàn)過程、實(shí)驗(yàn)結(jié)果及分析，同時(shí)給出了在實(shí)驗(yàn)過程中所遇到的問題描述，以及問題是否解決及待改進(jìn)的地方。本次實(shí)驗(yàn)所采用的編程工具為Visualstudio2008，編程語言采用C++。1引言：實(shí)驗(yàn)?zāi)康模簯?yīng)用參數(shù)估計(jì)方法解決實(shí)際問題。實(shí)驗(yàn)意義：通過本次實(shí)驗(yàn)，更加熟爛的掌握最小二成估計(jì)算法。使用實(shí)驗(yàn)中給出的數(shù)據(jù)選用適當(dāng)?shù)暮瘮?shù)（如適當(dāng)階次的多項(xiàng)式、高斯勢(shì)函數(shù)），用LS估計(jì)方法，擬合給定數(shù)據(jù)，給出擬合強(qiáng)度系數(shù)以及噪聲方差（設(shè)為獨(dú)立高斯噪聲）。2原理：y=a+bx+e，其中y、x可測(cè)，e是均值為0的隨機(jī)變量，a、b為未知參數(shù)。通過n次實(shí)驗(yàn)，得到測(cè)量數(shù)據(jù)yi和xi,i=1,2,…,n，確定未知參數(shù)a、b。使的估計(jì)稱為最小二乘（LS）估計(jì)，即殘差平方和最小的估計(jì)?；灸Ｐ停簩憺橄蛄啃问綖椋簩憺榫仃囆问綖椋浩渲校簲M合強(qiáng)度系數(shù)推導(dǎo)公式為：^β=(X’X)X’Y所以擬合后的函數(shù)值為：^Y=X^β殘差平方和計(jì)算公式如下：噪聲方差計(jì)算公式：Yawp=J(a)/(p-n)其中p為矩陣Y的行數(shù)，n為所使用的階數(shù)。3實(shí)驗(yàn)結(jié)果及分析：11010110102103……10n120202203……20n130302303…….30n…………….116716721673........167n....構(gòu)造X為：X=-152-152-148-148-166-158。。。.-115-108-120..Y=共167個(gè)數(shù)據(jù)輸入不同的N進(jìn)行實(shí)驗(yàn)，觀察不同的N值所對(duì)應(yīng)的殘差平方和及噪聲方差：下圖中黑線為原始數(shù)據(jù)所對(duì)應(yīng)的函數(shù)圖，紅色為N階擬合函數(shù)圖。以下列舉幾個(gè)比較具有代表性的N值所對(duì)應(yīng)的擬合函數(shù)及對(duì)應(yīng)的殘差平方和與噪聲方差。(1)取N=3實(shí)驗(yàn)結(jié)果如下：(2)取N=9(3)取N=13(4)取N=17(5)取N=40觀察以上N階擬合函數(shù)，發(fā)現(xiàn)當(dāng)N=17時(shí)擬合效果最好，即在N=17時(shí)殘差平方和最小。4小結(jié)試驗(yàn)中遇到的一些問題：(1)在寫求矩陣的逆矩陣的算法時(shí)，要先判斷該矩陣的行列式是否為0，由于邏輯錯(cuò)誤，導(dǎo)致程序進(jìn)入死循環(huán)。解決方法：不在程序中判斷矩陣的行列式是否為0，改為在實(shí)驗(yàn)過程中保證所涉及的矩陣行列式都不為0，再進(jìn)行運(yùn)算。(2)描述一個(gè)矩陣時(shí)要用一個(gè)數(shù)組及x,y來描述，但是這樣曾加了結(jié)構(gòu)復(fù)雜度，導(dǎo)致整體結(jié)構(gòu)復(fù)雜。解決方法：用一個(gè)結(jié)構(gòu)體封裝這個(gè)矩陣，結(jié)構(gòu)體里包含存放數(shù)據(jù)的數(shù)組及表示行數(shù)列數(shù)的x,y。(3)開始把數(shù)據(jù)類型定義為DOUBLE，但是在計(jì)算N很大時(shí)有可能發(fā)出溢出。解決方法：使用第三方高精度浮點(diǎn)數(shù)運(yùn)算庫(kù)函數(shù)。但是由于能力有限，該錯(cuò)誤還是存在，例如上面當(dāng)N=17時(shí)，殘差平方和很小，但擬合函數(shù)在后期卻顯示出與原函數(shù)偏差很大，估計(jì)就是由這一未解決的問題引起的。5參考文獻(xiàn)[1]孫榮恒.應(yīng)用數(shù)理統(tǒng)計(jì)[M].北京：科學(xué)出版社，2003.[2]夏普(英).Visualstudio2008從入門到精通[M].北京：清華大學(xué)出版社，2009.[3]同濟(jì)大學(xué)應(yīng)用數(shù)學(xué)系.線性代數(shù)[M].北京：高等教育出版社，2003.6附錄：主要程序代碼：(1)求矩陣轉(zhuǎn)置的算法：taticintMatrixAlgo::Transpose(Matrix<T>&matrix){ intnxTmp;//轉(zhuǎn)置后的x intnyTmp;//轉(zhuǎn)置后的y nxTmp=matrix.ny; nyTmp=matrix.nx; T*tmp_matrix_arry=newT[nxTmp*nyTmp]; if(tmp_matrix_arry==NULL) return0; for(intx=0;x<matrix.nx;++x) { for(inty=0;y<matrix.ny;++y) { tmp_matrix_arry[y*nyTmp+x]=matrix.matrix_arry[x*matrix.ny+y]; } } T*delete_tmp=matrix.matrix_arry; matrix.matrix_arry=tmp_matrix_arry; matrix.nx=nxTmp; matrix.ny=nyTmp; delete[]delete_tmp; return1;}(2)求矩陣逆矩陣的算法：staticintMatrixAlgo::Inverse(Matrix<T>&matrix){ if(matrix.nx!=matrix.ny) return0; intn=matrix.nx; int*is,*js,i,j,k,l,u,v; Td,p; is=newint[n]; js=newint[n];//開始計(jì)算逆矩陣 for(k=0;k<=n-1;k++) { d=0.0; for(i=k;i<=n-1;i++) { for(j=k;j<=n-1;j++) { l=i*n+j; p=fabs(matrix.matrix_arry[l]); if(p>d) { d=p; is[k]=i; js[k]=j; } } } if(is[k]!=k) { for(j=0;j<=n-1;j++) { u=k*n+j; v=is[k]*n+j; p=matrix.matrix_arry[u]; matrix.matrix_arry[u]=matrix.matrix_arry[v]; matrix.matrix_arry[v]=p; } } if(js[k]!=k) { for(i=0;i<=n-1;i++) { u=i*n+k; v=i*n+js[k]; p=matrix.matrix_arry[u]; matrix.matrix_arry[u]=matrix.matrix_arry[v]; matrix.matrix_arry[v]=p; } } l=k*n+k; matrix.matrix_arry[l]=1.0/matrix.matrix_arry[l]; for(j=0;j<=n-1;j++) { if(j!=k) { u=k*n+j; matrix.matrix_arry[u]=matrix.matrix_arry[u]*matrix.matrix_arry[l]; } } for(i=0;i<=n-1;i++) { if(i!=k) { for(j=0;j<=n-1;j++) { if(j!=k) { u=i*n+j; matrix.matrix_arry[u]=matrix.matrix_arry[u]-matrix.matrix_arry[i*n+k]*matrix.matrix_arry[k*n+j]; } } } } for(i=0;i<=n-1;i++) { if(i!=k) { u=i*n+k; matrix.matrix_arry[u]=-matrix.matrix_arry[u]*matrix.matrix_arry[l]; } } } for(k=n-1;k>=0;k--) { if(js[k]!=k) { for(j=0;j<=n-1;j++) { u=k*n+j; v=js[k]*n+j; p=matrix.matrix_arry[u]; matrix.matrix_arry[u]=matrix.matrix_arry[v]; matrix.matrix_arry[v]=p; } } if(is[k]!=k) for(i=0;i<=n-1;i++) { u=i*n+k; v=i*n+is[k]; p=matrix.matrix_arry[u]; matrix.matrix_arry[u]=matrix.matrix_arry[v]; matrix.matrix_arry[v]=p; } } delete[]is;delete[]js; return1;}(3)矩陣乘法算法：staticMatrix<T>*MatrixAlgo::Multiplication(Matrix<T>&matrix_lift,Matrix<T>&matrix_right){ Matrix<T>*tmp_matrix=newMatrix<T>; if(matrix_lift.ny!=matrix_right.nx)//非法 { returnNULL; } intnxTmp=matrix_lift.nx;//乘法后的x intnyTmp=matrix_right.ny;//乘法后的y T*tmp_matrix_arry=newT[nxTmp*nyTmp]; if(tmp_matrix_arry==NULL) { returnNULL; } inti,j,l,u; for(i=0;i<=nxTmp-1;i++) { for(j=0;j<=nyTmp-1;j++) { u=i*nyTmp+j; tmp_matrix_arry[u]=0; for(l=0;l<=matrix_lift.ny-1;l++) { tmp_matrix_arry[u]=(tmp_matrix_arry[u]+matrix_lift.matrix_arry[i*matrix_lift.ny+l]*matrix_right.matrix_arry[l*nyTmp+j]); } } } //T*delete_tmp=matrix.matrix_arry; tmp_matrix->matrix_arry=tmp_matrix_arry; tmp_matrix->nx=nxTmp; tmp_matrix->ny=nyTmp; //delete[]delete_tmp; returntmp_matrix;}(4)計(jì)算^β，依次調(diào)用函數(shù)順序?yàn)椋?intSetXMatrix(intrank,intspace); intSetYMatrix(int*arry,intn); voidClear(); voidXClear();具體函數(shù)算法見附