数理统计课程实验报告Word格式文档下载.docx
- 文档编号:7239313
- 上传时间:2023-05-08
- 格式:DOCX
- 页数:15
- 大小:221.17KB
数理统计课程实验报告Word格式文档下载.docx
《数理统计课程实验报告Word格式文档下载.docx》由会员分享,可在线阅读,更多相关《数理统计课程实验报告Word格式文档下载.docx(15页珍藏版)》请在冰点文库上搜索。
拟合强度系数推导公式为:
^β=(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<
T>
&
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<
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<
if(matrix.nx!
=matrix.ny)
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<
i++)
{
for(j=k;
j<
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;
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]!
for(i=0;
u=i*n+k;
v=i*n+js[k];
l=k*n+k;
matrix.matrix_arry[l]=1.0/matrix.matrix_arry[l];
for(j=0;
if(j!
matrix.matrix_arry[u]=matrix.matrix_arry[u]*matrix.matrix_arry[l];
for(i=0;
if(i!
for(j=0;
if(j!
{
u=i*n+j;
matrix.matrix_arry[u]=matrix.matrix_arry[u]-matrix.matrix_arry[i*n+k]*matrix.matrix_arry[k*n+j];
}
matrix.matrix_arry[u]=-matrix.matrix_arry[u]*matrix.matrix_arry[l];
}
}
for(k=n-1;
k>
=0;
k--)
v=js[k]*n+j;
v=i*n+is[k];
delete[]is;
delete[]js;
}
(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];
inti,j,l,u;
for(i=0;
=nxTmp-1;
i++)
=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;
nx=nxTmp;
ny=nyTmp;
//delete[]delete_tmp;
returntmp_matrix;
(4)计算^β,依次调用函数顺序为:
intSetXMatrix(intrank,intspace);
intSetYMatrix(int*arry,intn);
voidClear();
voidXClear();
具体函数算法见附件“实验”。
voidYClear();
intComputeBeta();
longdouble>
*GetBeta();
注:
要先设置Y,再根据Y设置X,否则会导致失败。
(因为Y中的数据有缺失。
)
(5)计算^Y:
*GetFileMatrix();
(6)计算残差平方和,调用函数:
longGetE();
//获取E具体函数算法见附件“实验”。
(7)计算噪声方差,调用函数:
floatGetYawp();
可执行程序见附件。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数理统计 课程 实验 报告