最小二乘法 C++实现

#include <iostream>
#include "../common.h"
#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]={1,2,3,4,5};
	double arry2[5]={1.1,1.2,1.3,1.4,1.5};
	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]);
    printf("y6=%lf\n",coefficient[1]+coefficient[2]*6+coefficient[3]*6*6);
	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; 
}

这段代码是转载自https://blog.youkuaiyun.com/lsh_2013/article/details/46697625

#include <iostream>
#include "../common.h"
#include <vector>
#include <cmath>
using namespace std;
 
//最小二乘拟合相关函数定义
float sum(vector<float> Vnum, int n);
float MutilSum(vector<float> Vx, vector<float> Vy, int n);
float RelatePow(vector<float> Vx, int n, int ex);
float RelateMutiXY(vector<float> Vx, vector<float> Vy, int n, int ex);
void EMatrix(vector<float> Vx, vector<float> Vy, int n, int ex, float coefficient[]);
void CalEquation(int exp, float coefficient[]);
float F(float c[],int l,int m);
float Em[11][11];   //修改,防止溢出
 
//主函数，这里将数据拟合成5次曲线
int main(int argc, char* argv[])
{
	float coefficient[10]; //修改,防止溢出
	memset(coefficient,0,sizeof(float)*5);
	vector<float> vx,data;
    data.push_back(1.1);
    data.push_back(1.2);
    data.push_back(1.3);
    data.push_back(1.5);
    data.push_back(1.7);
    data.push_back(1.9);
    data.push_back(2.1);
    data.push_back(2.3);
    for(int i = 0; i<data.size(); ++i)
    {
        vx.push_back(i+1);
    }


	EMatrix(vx,data,data.size(),5,coefficient);
	printf("拟合方程为：y = %lf + %fx + %fx^2 + %fx^3 + %fx^4 + %fx^5\n",coefficient[1],coefficient[2],\
        coefficient[3],coefficient[4],coefficient[5],coefficient[6]);
    int x = 9;
    printf("y(%d)=%lf\n",x,coefficient[1]+coefficient[2]*x+coefficient[3]*pow(x,2)+\
        coefficient[4]*pow(x,3) + coefficient[5]*pow(x,4) + coefficient[6]*pow(x,5));
	return 0;
}
//累加
float sum(vector<float> Vnum, int n)
{
	float dsum=0;
	for (int i=0; i<n; i++)
	{
		dsum+=Vnum[i];
	}
	return dsum;
}
//乘积和
float MutilSum(vector<float> Vx, vector<float> Vy, int n)
{
	float dMultiSum=0;
	for (int i=0; i<n; i++)
	{
		dMultiSum+=Vx[i]*Vy[i];
	}
	return dMultiSum;
}
//ex次方和
float RelatePow(vector<float> Vx, int n, int ex)
{
	float ReSum=0;
	for (int i=0; i<n; i++)
	{
		ReSum+=pow(Vx[i],ex);
	}
	return ReSum;
}
//x的ex次方与y的乘积的累加
float RelateMutiXY(vector<float> Vx, vector<float> Vy, int n, int ex)
{
	float dReMultiSum=0;
	for (int i=0; i<n; i++)
	{
		dReMultiSum+=pow(Vx[i],ex)*Vy[i];
	}
	return dReMultiSum;
}
//计算方程组的增广矩阵
void EMatrix(vector<float> Vx, vector<float> Vy, int n, int ex, float 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, float coefficient[])
{
	for(int k=1;k<exp;k++) //消元过程
	{
		for(int i=k+1;i<exp+1;i++)
		{
			float 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函数调用
float F(float c[],int l,int m)
{
	float sum=0;
	for(int i=l;i<=m;i++)
		sum+=Em[l-1][i]*c[i];
	return sum; 
}

修改之后拟合5阶曲线，原代码当中直接改多阶会有溢出，这里修改了Em的定义以及main()函数中的部分代码，其它地方均为做修改。