最小二乘法的编程实现.docx
- 文档编号:9701412
- 上传时间:2023-05-20
- 格式:DOCX
- 页数:17
- 大小:73.72KB
最小二乘法的编程实现.docx
《最小二乘法的编程实现.docx》由会员分享,可在线阅读,更多相关《最小二乘法的编程实现.docx(17页珍藏版)》请在冰点文库上搜索。
最小二乘法的编程实现
1、最小二乘法:
1)(用
方法计算逆矩阵)
#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的行列式 { 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)^(i+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 { for(j=0;j { 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.0469030.0604030.154991];%由C语言计算出的fi值 fi=[1;-fi']; ts2=idpoly(fi',[]); B=Y(1: N-n); plot(B); compare(B,ts2,1); 2)用构造矩阵 的方法求逆: (列主元高斯消去法) 部分主要程序: 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行消去其他行{ {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]; } } } voidNormal(doublea[n][M]) { inti,j; doubletemp; for(i=0;i { temp=a[i][i]; if(temp! =1.0) { for(j=i;j { a[i][j]=a[i][j]/temp; } } } } voidInverse(doublea[n][M]) { inti,j,k; doubletemp; for(i=n-1;i>=0;i--) { for(j=i-1;j>=0;j--) { temp=a[j][i]; for(k=i;k { a[j][k]=a[j][k]-temp*a[i][k]; } } } } voidGetInverse(doublea[n][M],doubleb[n][n]) { inti,j; for(i=0;i { for(j=0;j { b[i][j]=a[i][j+n]; } } } voidmain() { inti,j;doubletemp; doublesunpot[N];doublesunpot_normal[N];doublematrixr[n]; doublematrix1[n][n]={0};//方阵 doublematrix[n][M]={0};//构造矩阵 doubleinverse[n][n]={0};doublefi[n]; Getdata(sunpot);//从txt文档中读取每年的太阳黑子数据 Normal(sunpot,sunpot_normal);//数据归一化 Matrix(sunpot_normal,matrix1,matrixr);//组建方程 BuildMatrix(matrix,matrix1); for(i=0;i { Judge(matrix,i); Gauss(matrix,i); } for(i=0;i { if(matrix[i][i]==0) { printf("矩阵不可逆! \n"); exit(-1); } } Normal(matrix); Inverse(matrix); GetInverse(matrix,inverse); printf("逆矩阵: \n"); for(i=0;i { for(j=0;j { printf("%10lf",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]); } }
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 最小二乘法 编程 实现
![提示](https://static.bingdoc.com/images/bang_tan.gif)