最小二乘法拟合多项式原理以及c++实现

转载请注明出处：/lsh_/article/details/46697625

最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。

c++实现代码如下：

#include <iostream>#include <vector>#include <cmath>using namespace std;//最小二乘拟合相关函数定义double sum(vector<double> Vnum, int n);double MutilSum(vector<double> Vx, vector<double> Vy, int n);double RelatePow(vector<double> Vx, int n, int ex);double RelateMutiXY(vector<double> Vx, vector<double> Vy, int n, int ex);void EMatrix(vector<double> Vx, vector<double> Vy, int n, int ex, double coefficient[]);void CalEquation(int exp, double coefficient[]);double F(double c[],int l,int m);double Em[6][4];//主函数，这里将数据拟合成二次曲线int main(int argc, char* argv[]){double arry1[5]={0,0.25,0,5,0.75};double arry2[5]={1,1.283,1.649,2.212,2.178};double coefficient[5];memset(coefficient,0,sizeof(double)*5);vector<double> vx,vy;for (int i=0; i<5; i++){vx.push_back(arry1[i]);vy.push_back(arry2[i]);}EMatrix(vx,vy,5,3,coefficient);printf("拟合方程为：y = %lf + %lfx + %lfx^2 \n",coefficient[1],coefficient[2],coefficient[3]);return 0;}//累加double sum(vector<double> Vnum, int n){double dsum=0;for (int i=0; i<n; i++){dsum+=Vnum[i];}return dsum;}//乘积和double MutilSum(vector<double> Vx, vector<double> Vy, int n){double dMultiSum=0;for (int i=0; i<n; i++){dMultiSum+=Vx[i]*Vy[i];}return dMultiSum;}//ex次方和double RelatePow(vector<double> Vx, int n, int ex){double ReSum=0;for (int i=0; i<n; i++){ReSum+=pow(Vx[i],ex);}return ReSum;}//x的ex次方与y的乘积的累加double RelateMutiXY(vector<double> Vx, vector<double> Vy, int n, int ex){double dReMultiSum=0;for (int i=0; i<n; i++){dReMultiSum+=pow(Vx[i],ex)*Vy[i];}return dReMultiSum;}//计算方程组的增广矩阵void EMatrix(vector<double> Vx, vector<double> Vy, int n, int ex, double coefficient[]){for (int i=1; i<=ex; i++){for (int j=1; j<=ex; j++){Em[i][j]=RelatePow(Vx,n,i+j-2);}Em[i][ex+1]=RelateMutiXY(Vx,Vy,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];}//供CalEquation函数调用double F(double c[],int l,int m){double sum=0;for(int i=l;i<=m;i++)sum+=Em[l-1][i]*c[i];return sum; }