最小二乘法的多项式拟合代码

     最近的项目中用到了最小二乘法的多项式的拟合，经过网上资料查找首先了解了什么是最小二乘法，后来主要看到了最小二乘法的拟合，如果是一次就是最小二乘法线性拟合，如果是多次就是曲线拟合。这里我不想多说关于最小二乘法的原理。经过对资料的阅读和了解我试着使用c语言实现了多项式拟合，当然这里我没有实现压缩，随着高阶逐渐增加拟合的曲线也会存在病态，所以有时候对比较高的高阶需要压缩后在拟合。这里是没有经过压缩的拟合，贴出来给大家一个参考：

#include "stdafx.h" #include "stdio.h" #include "math.h" #include "stdlib.h" #include "windows.h" double sum(double * dNumarry,int n); double MutilSum(double* dX,double *dY,int n); double RelatePow(double *dx,int n,int ex); double RelateMutiXY(double *dx,double*dy,int n,int ex); void EMatrix(double *dx,double*dy,int n,int ex,double coefficient[]); void CalEquation(int exp,double coefficient[] ); double F(double c[],int l,int m); double Em[10][10]; //所构造的多项式的系数 int main(int argc, char* argv[]) { double arry1[9]={1,3,4,5,6,7,8,9,10}; double arry2[9]={10,5,4,2,1,1,2,3,4}; double coefficient[10]; ZeroMemory(coefficient,sizeof(double)*10); // double dSumarry1= sum(arry1,7); // printf("darry1=%lf/n",dSumarry1); // // double dSumarry2=sum(arry2,7); // printf("darry2=%lf/n",dSumarry2); // // double dMultiarry1=MutilSum(arry1,arry1,7); // printf("dMultiarry1=%lf/n",dMultiarry1); // double dMultiarry12=MutilSum(arry1,arry2,7); // printf("dMultiarry2=%lf/n",dMultiarry12);//,coefficient[4],coefficient[5] //+ %lfx^3 + %lfx^4 EMatrix(arry1,arry2,9,3,coefficient); printf("拟合方程为：y = %lf + %lfx + %lfx^2 /n",coefficient[1],coefficient[2],coefficient[3]); return 0; } double sum(double * dNumarry,int n) { double *dTemp= new double[n]; double dSum=0; dTemp=dNumarry; for (int i=0;i<n;i++) { dSum+=dTemp[i]; } return dSum; } double MutilSum(double* dX,double *dY,int n) { double * dXTemp= new double [n]; double * dYTemp= new double [n]; double dMultiSum=0; dXTemp=dX; dYTemp=dY; for (int i=0;i<n;i++) { dMultiSum += dX[i]*dY[i]; } return dMultiSum; } double RelatePow(double *dx,int n,int ex) { double * dTemp =new double[n]; double ReSum=0; dTemp=dx; for (int j=0;j<n;j++) { ReSum+=pow(dTemp[j],ex); } return ReSum; } double RelateMutiXY(double *dx,double*dy,int n,int ex) { double * dXTemp= new double [n]; double * dYTemp= new double [n]; double dReMultiSum=0; dXTemp=dx; dYTemp=dy; for (int i=0;i<n;i++) { dReMultiSum+=pow(dXTemp[i],ex)*dYTemp[i]; } return dReMultiSum; } //计算方程组的增广矩阵 void EMatrix(double *dx,double*dy,int n,int ex,double coefficient[] ) { double * dXTemp= new double [n]; double * dYTemp= new double [n]; dXTemp=dx; dYTemp=dy; for(int i=1;i<=ex;i++) { for(int j=1;j<=ex;j++) { Em[i][j]=RelatePow(dXTemp,n,i+j-2); } Em[i][ex+1]=RelateMutiXY(dXTemp,dYTemp,n,i-1); } Em[1][1]=n; CalEquation(ex,coefficient); } //求解方程 void CalEquation(int exp,double coefficient[] ) { for(int k=1;k<exp;k++) //消元过程 { for(int i=k+1;i<exp+1;i++) { double p1=0; if(Em[k][k]!=0) p1=Em[i][k]/Em[k][k]; for(int j=k;j<exp+2;j++) Em[i][j]=Em[i][j]-Em[k][j]*p1; } } coefficient[exp]=Em[exp][exp+1]/Em[exp][exp]; for(int l=exp-1;l>=1;l--) //回代求解 coefficient[l]=(Em[l][exp+1]-F(coefficient,l+1,exp))/Em[l][l]; } double F(double c[],int l,int m)//供CalEquation函数调用 { double sum=0; for(int i=l;i<=m;i++) sum+=Em[l-1][i]*c[i]; return sum; }

这段代码你可以用vc建立一个console的程序然后编译运行，我这里主要是拟合二次多项式，简单说明一下EMatrix这个函数的ex要加1就是说如果要求四次就要ex要写5 。其他的就需要自己去阅读代码了，其实知道最小二乘法的原理看代码就会很容易的唷。