一般最小二乘法 c语言,C语言编写最小二乘法

最小二乘法简介：1801年，意大利天文学家朱赛普·皮亚齐发现了第一颗小行星谷神星。经过40天的跟踪观测后，由于谷神星运行至太阳背后，使得皮亚齐失去了谷神星的位置。随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星，但是根据大多数人计算的结果来寻找谷神星都没有结果。时年24岁的高斯也计算了谷神星的轨道。奥地利天文学家海因里希·奥尔伯斯根据高斯计算出来的轨道重新发现了谷神星。

高斯使用的最小二乘法的方法发表于1809年他的著作《天体运动论》中，而法国科学家勒让德于1806年独立发现“最小二乘法”，但因不为时人所知而默默无闻。两人曾为谁最早创立最小二乘法原理发生争执。

1829年，高斯提供了最小二乘法的优化效果强于其他方法的证明，见高斯-马尔可夫定理。计算方法是门联系计算机技术与数学的学科，它提供了用计算机的思维解决数学问题的方法。C语言是应用广泛的编程语言，它虽是高级语言但也可进行底层编程。本代码结合了基本的最小二乘法思想，求解正规方程用的是列主元消去的思想。

#include

#include

int M;

int N;

double x[20],a[10];

double y[20],xx[10][10],yy[10];

void Q_liheF(int m,int n);

void Q_zhi(int m);

void S_chu(int m);

void B_jiao(int m,int k);

void S_qu(int m,int k);

void huidai(int m);

void main()

{

int i;

printf("已知M为拟合的次数，N+1为数据个数,输入M,N的值:\n");

scanf("%d,%d",&M,&N);

for(i=0;i<=N;i++)

scanf("%lf,%lf",&x[i],&y[i]);

printf("\n拟合原始数据:\n");

for(i=0;i<=N;i++)

printf("f(%lf)=%lf\n",x[i],y[i]);

for(i=0;i<=M;i++)

a[i]=0;

Q_liheF(M,N);

Q_zhi(M);

S_chu(M);

}

void S_chu(int m)

{

int i;

printf("\n\tf(x)=");

for(i=0;i<=m;i++)

{

if(a[i]>=0)

{

if(i!=0)

printf("+%.4lfx^%d",a[i],i);

else

printf("%.4lf",a[i]);

}

else

{

if(i!=0)

printf("%.4lfx^%d",a[i],i);

else

printf("%.4lf",a[i]);

}

}

printf("\n\n");

}

void Q_liheF(int m,int n)

{

int i,j,k;

double t1,t2;

float h;

for(i=0;i<=m;i++)

{

t1=0;

t2=0;

for(j=i;j<=m;j++)

{

t1=0;

for(k=0;k<=n;k++)

{

h=i+j;

t1+=pow(x[k],h);

}

xx[i][j]=t1;

xx[j][i]=t1;

}

for(k=0;k<=n;k++)

{

h=i;

t2+=y[k]*pow(x[k],h);

}

yy[i]=t2;

}

printf("\n正规方程组:\n");

for(i=0;i<=m;i++)

{

for(j=0;j<=m;j++)

printf("%.4lf ",xx[i][j]);

printf("**%.4lf\n",yy[i]);

}

}

void Q_zhi(int m)

{

int i,flag;

flag=0;

for(i=0;i<=m-1;i++)

{

B_jiao(m,i);

if(xx[i][i]==0)

{

flag=1;

break;

}

S_qu(m,i);

}

if(flag!=1)

huidai(m);

else

printf("\n******拟合有误无法求解!*****\n");

}

void B_jiao(int m,int k)

{

int i,j,q;

double t;

for(i=k,q=i;i<=m-1;i++)

{

if(fabs(xx[i][k])

q=i+1;

}

if(q!=k)

{

for(j=k;j<=m;j++)

{

t=xx[k][j];

xx[k][j]=xx[q][j];

xx[q][j]=t;

}

t=yy[k];

yy[k]=yy[q];

yy[q]=t;

}

}

void S_qu(int m,int k)

{

int i,j;

double l;

for(i=k+1;i<=m;i++)

{

l=xx[i][k]/xx[k][k];

for(j=k;j<=m;j++)

{

xx[i][j]=xx[i][j]-l*xx[k][j];

}

yy[i]=yy[i]-l*yy[k];

}

}

void huidai(int m)

{

int i,j;

double A=0;

a[m]=yy[m]/xx[m][m];

for(i=m-1;i>=0;i--)

{

for(j=m;j>=i+1;j--)

A+=a[j]*xx[i][j];

a[i]=(yy[i]-A)/xx[i][i];

A=0;

}

}

运行结果：