数理统计课程实验报告.docx
《数理统计课程实验报告.docx》由会员分享，可在线阅读，更多相关《数理统计课程实验报告.docx（15页珍藏版）》请在冰豆网上搜索。
数理统计课程实验报告
高等数理统计课程实验
【摘要】本实验报告描述了用最小二成估计算法解决实际问题中参数估计的过程。
包括引言、实验原理、实验过程、实验结果及分析，同时给出了在实验过程中所遇到的问题描述，以及问题是否解决及待改进的地方。
本次实验所采用的编程工具为Visualstudio2008，编程语言采用C++。
1引言：
实验目的：
应用参数估计方法解决实际问题。
实验意义：
通过本次实验，更加熟烂的掌握最小二成估计算法。
使用实验中给出的数据选用适当的函数（如适当阶次的多项式、高斯势函数），用LS估计方法，拟合给定数据，给出拟合强度系数以及噪声方差（设为独立高斯噪声）。
2原理：
y=a+bx+，其中y、x可测，是均值为0的随机变量，a、b为未知参数。
通过n次实验，得到测量数据yi和xi,i=1,2,…,n，
确定未知参数a、b。
使
的估计称为最小二乘（LS）估计，即残差平方和最小的估计。
基本模型：
写为向量形式为：
写为矩阵形式为：
其中：
拟合强度系数推导公式为：
^β=（X’X）X’Y
所以拟合后的函数值为：
^Y=X^β
残差平方和计算公式如下：
噪声方差计算公式：
Yawp=J（a）/（p-n）
其中p为矩阵Y的行数，n为所使用的阶数。
3实验结果及分析：
构造X为：
X=
Y=
共167个数据
输入不同的N进行实验，观察不同的N值所对应的残差平方和及噪声方差：
下图中黑线为原始数据所对应的函数图，红色为N阶拟合函数图。
以下列举几个比较具有代表性的N值所对应的拟合函数及对应的残差平方和与噪声方差。
（1）取N=3
实验结果如下：
（2）取N=9
（3）取N=13
（4）取N=17
（5）取N=40
观察以上N阶拟合函数，发现当N=17时拟合效果最好，即在N=17时残差平方和最小。
4小结
试验中遇到的一些问题：
（1）在写求矩阵的逆矩阵的算法时，要先判断该矩阵的行列式是否为0，由于逻辑错误，导致程序进入死循环。
解决方法：
不在程序中判断矩阵的行列式是否为0，改为在实验过程中保证所涉及的矩阵行列式都不为0，再进行运算。
（2）描述一个矩阵时要用一个数组及x,y来描述，但是这样曾加了结构复杂度，导致整体结构复杂。
解决方法：
用一个结构体封装这个矩阵，结构体里包含存放数据的数组及表示行数列数的x,y。
（3）开始把数据类型定义为DOUBLE，但是在计算N很大时有可能发出溢出。
解决方法：
使用第三方高精度浮点数运算库函数。
但是由于能力有限，该错误还是存在，例如上面当N=17时，残差平方和很小，但拟合函数在后期却显示出与原函数偏差很大，估计就是由这一未解决的问题引起的。
5参考文献
[1]孙荣恒.应用数理统计[M].北京：
科学出版社，2003.
[2]夏普（英）.Visualstudio2008从入门到精通[M].北京：
清华大学出版社，2009.
[3]同济大学应用数学系.线性代数[M].北京：
高等教育出版社，2003.
6附录：
主要程序代码：
（1）求矩阵转置的算法：
taticintMatrixAlgo:
:
Transpose（Matrix&matrix）
{
intnxTmp;//转置后的x
intnyTmp;//转置后的y
nxTmp=matrix.ny;
nyTmp=matrix.nx;
T*tmp_matrix_arry=newT[nxTmp*nyTmp];
if（tmp_matrix_arry==NULL）
return0;
for（intx=0;x{
for（inty=0;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&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];
//开始计算逆矩阵
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*MatrixAlgo:
:
Multiplication（Matrix&matrix_lift,Matrix&matrix_right）
{
Matrix*tmp_matrix=newMatrix;
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）计算^β，依次调用函数顺序为：
intSetXMatrix（intrank,intspace）;
intSetYMatrix（int*arry,intn）;
voidClear（）;
voidXClear（）;具体函数算法见附件“实验”。
voidYClear（）;
intComputeBeta（）;
Matrix*GetBeta（）;
注：
要先设置Y，再根据Y设置X，否则会导致失败。
（因为Y中的数据有缺失。
）
（5）计算^Y：
Matrix*GetFileMatrix（）;
（6）计算残差平方和，调用函数：
longGetE（）;//获取E具体函数算法见附件“实验”。
（7）计算噪声方差，调用函数：
floatGetYawp（）;
可执行程序见附件。
升级会员 