最小二乘法的编程实现.docx
- 文档编号:3105310
- 上传时间:2022-11-17
- 格式:DOCX
- 页数:14
- 大小:31.52KB
最小二乘法的编程实现.docx
《最小二乘法的编程实现.docx》由会员分享,可在线阅读,更多相关《最小二乘法的编程实现.docx(14页珍藏版)》请在冰豆网上搜索。
最小二乘法的编程实现
1、最小二乘法:
1)(用Aat方法计算逆矩阵)
IA
#include
#include
#include
#include
#include
#defineN200
#definen9
voidGetdata(doublesun[N])//从txt文档中读取数据(小数)
{
chardata;
charsunpot[10]={0000000000};//为了预防结果出现'烫'字inti=0,j=0;
doubled;
FILE*fp=fopen("新建文本文档.txt","r〞);
if(!
fp)
{
printf("can'topenfile\n");
}
while(!
feof(fp))
{
data=fgetc(fp);
if(data!
='\n')
{
sunpot[i]=data;
i++;
}
elseif(data=='\n')
{
sunpot[i]='\0';//给定结束符
d=atof(sunpot);//将字符串转换成浮点数sun[j]=d;
j++;
i=0;//将i复位
}
}
}
voidNormal(doublesun[N],doublesun1[N])//将数据进行标准化{
doublemean,temp=0,variance=0;
inti;
for(i=0;i temp=temp+sun[i]; mean=temp/N; for(i=0;i { sun1[i]=sun[i]-mean; } for(i=0;i { variance=variance+sun1[i]*sun1[i]; } variance=variance/N; variance=sqrt(variance); for(i=0;i { sun1[i]=sun[i]/variance;//减去均值除以均方差进行回一化 } } voidMatrix(doublesun[N],doublematrixl[n][n],doublematrixr[n])//组建方程的左右两个矩阵 { doublematrix1[N-n][n];//方程左边系数矩阵 doublematrix_trans[n][N-n]; doublematrix2[N-n];//方程右边系数矩阵 inti,j,k; doubletemp; for(i=0;i { for(j=0;j { matrix1[i][j]=sun[n+i-j-1];//此处与matlab不同,matlab是从1开始的 } } for(i=0;i matrix2[i]=sun[i+n]; for(i=0;i { for(j=0;j { matrix_trans[i][j]=matrix1[j][i]; for(i=0;i ( for(k=0;k ( temp=0; for(j=0;j ( temp+=matrix_trans[i][j]*matrix1[j][k]; ) matrixl[i][k]=temp; ) ) for(i=0;i ( temp=0; for(j=0;j ( temp+=matrix_trans[i][j]*matrix2[j]; ) matrixr[i]=temp; ) ) voidExchangeLine(doublea[n][n],inti,intj,intp)//交换矩阵第i行和第j行 ( intk; doubletemp; for(k=0;k<=p-1;k++) ( temp=a[j][k]; a[j][k]=a[i][k]; a[i][k]=temp; ) ) voidGauss(doublea[n][n],intk,intp)//高斯消元,用第k行消去其他行第一个非零元素 ( inti,j; doublem; for(i=k+1;i ( if(a[i][k]==0) i++; m=a[i][k]/a[k][k]; for(j=k;j ( a[i][j]=a[i][j]-m*a[k][j]; ) ) ) doubleDeterminant(doublea[n][n],intp)//求a的行歹U式 ( intk,i,j; doubleresult=1.0; for(k=0;k ( i=k;j=k; while(a[i][j]==0) i++; if((i k) ( ExchangeLine(a,k,i,p); result=(-1)*result;//换行之后行列式变号一次 ) elseif(i>=p)//即找不到可以交换的行 ( result=0; break; ) Gauss(a,k,p);//对该列进行消元 ) for(i=0;i ( result=result*a[i][i]; ) returnresult; ) voidAdjoint_matrix(doublea[n][n],doubleb[n][n],intp,inti,intj)//求a[i][j]元素的伴随矩阵 ( doubletemp[n][n]=({0,0,0),(0,0,0),(0,0,0});//实际上只要n-1xn-1维,但为了了调用方便,此处写为了nxn 维 intk,m,l,r; b[j][i]=1; for(k=0,l=0;k if(l==i)l++; for(m=0,r=0;m ( if(r==j)r++; temp[k][m]=a[l][r]; ) ) b[j][i]=b[j][i]*Determinant(temp,n-1); if((i+j)%2==1)//(-1)1i+j) b[j][i]=-b[j][i]; ) voidInverse(doublea[n][n],doubleb,doublec[n][n])//求逆最后一步运算 ( inti,j; for(i=0;i ( for(j=0;j ( c[i][j]=a[i][j]/b; ) ) ) voidmain() ( intp=n,i,j; doublesunpot[N],sunpot_normal[N]; doublefi[n]; doublematrix[n][n];//要换成返回的方阵 doublematrixr[n]; doubleadjoint[n][n],inverse[n][n]; doubledeterminant,temp; Getdata(sunpot);//从txt文档中读取每年的太阳黑子数据 Normal(sunpot,sunpot_normal);//数据回一化 Matrix(sunpot_normal,matrix,matrixr);//组建方程 for(i=0;i ( Adjoint_matrix(matrix,adjoint,p,i,j); ) ) determinant=Determinant(matrix,p); if(determinant==0) printf("该矩阵不存在逆矩阵! \n"); printf("行列式: \n"); printf("%10lf\n",determinant); Inverse(adjoint,determinant,inverse); printf("逆矩阵: \n"); for(i=0;i ( for(j=0;j ( printf("%12lf〞,inverse[i][j]); ) printf("\n"); ) printf("fi: \n"); for(i=0;i ( temp=0; for(j=0;j ( temp+=inverse[i][j]*matrixr[j]; ) fi[i]=temp; printf("%10lf〞,fi[i]); ) printf("\n"); ) 计算结果: fi: 1.327258-0.6226830.0185710.129074-0.1412780.103665-0.0469030.0604030.154991 下面用matlab进行拟合度计算: clear M=load('太阳黑子数.txt'); sunpot=(M(: 2)); N=200; n=9; [Y,mu,sigma]=zscore(sunpot); B(1: N-n)=Y(n+1: N); B=B'; fi=[1.327258-0.6226830.0185710.129074-0.1412780.103665-0.046903 0.0604030.154991];%由C语言计算出的fi值 fi=[1;-fi']; ts2=idpoly(fi',[]); B=Y(1: N-n); plot(B); compare(B,ts2,1); 2)用构造矩阵[A;I][I;A1]的方法求逆: (列主元高斯消去法) 局部主要程序: voidBuildMatrix(doublea[n][M],doubleb[n][n]) { inti,j; for(i=0;i { for(j=0;j { a[i][j]=b[i][j]; } a[i][i+n]=1.0; voidExchangeLine(doublea[n][M],inti,intj) ( intk; doubletemp; for(k=0;k ( temp=a[i][k]; a[i][k]=a[j][k]; a[j][k]=temp; } } voidJudge(doublea[n][M],inti) ( intj,m=0;//m=0来表示是否需要换行 doubletemp=a[i][i]; for(j=i+1;j ( if(fabs(a[j][i])>fabs(temp)) ( temp=a[j][i]; m=j; } } if(m) ExchangeLine(a,i,m); } voidGauss(doublea[n][M],intk)〃高斯消元法,用第k行消去其他行{
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 最小二乘法 编程 实现
![提示](https://static.bdocx.com/images/bang_tan.gif)