最小二乘法的编程实现.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;itemp=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;imatrix2[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（（ik）
（
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;kif（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行消去其他行{
升级会员 