c++ 多项式拟合polyfit实现

在numpy中有很方便的多项式拟合函数，c++中则没有，需要自己实现。

在网上参考了一些，但是都不能令人满意，其中有借用Opencv的，我就做一个曲线拟合，引入庞大的Opencv有点小题大做了。

也有使用Eigen库的，这个还差强人意，先记录下来。

double polyeval(Eigen::VectorXd coeffs, double x) {
  double result = 0.0;
  for (int i = 0; i < coeffs.size(); i++) {result += coeffs[i]*pow(x, i);}
  return result;
}

Eigen::VectorXd polyfit(Eigen::VectorXd xvals, Eigen::VectorXd yvals, int order) {
  assert(xvals.size() == yvals.size());
  assert(order >= 1 && order <= xvals.size() - 1);
  Eigen::MatrixXd A(xvals.size(), order + 1);
  for (int i = 0; i < xvals.size(); i++) {
    A(i, 0) = 1.0;
  }
  for (int j = 0; j < xvals.size(); j++) {
    for (int i = 0; i < order; i++) {
      A(j, i + 1) = A(j, i) * xvals(j);
    }
  }
  auto Q = A.householderQr();
  auto result = Q.solve(yvals);
  return result;
}

Eigen::VectorXd xvals(6);
Eigen::VectorXd yvals(6);
xvals << 9.261977, -2.06803, -19.6663, -36.868, -51.6263, -66.3482;
yvals << 5.17,     -2.25,    -15.306,  -29.46,  -42.85,   -57.6116;
auto coeffs= polyfit(xvals,yvals,3);
std::cout<< "Y(16)=" << polyeval(coeffs,16) << endl;

这个也可以参考：

#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; 
}

我心中的多项式拟合应该是基于模板类的通用多项式拟合函数。

下面是一段自己写的测试程序：
 

#include <Eigen/Dense>
#include <Eigen/Core>

double polyeval(Eigen::VectorXd coeffs, double x) {
	double result = 0.0;
	for (int i = 0; i < coeffs.size(); i++) { result += coeffs[i] * pow(x, i); }
	return result;
}

Eigen::VectorXd polyfit(Eigen::VectorXd xvals, Eigen::VectorXd yvals, int order) {
	assert(xvals.size() == yvals.size());
	assert(order >= 1 && order <= xvals.size() - 1);
	Eigen::MatrixXd A(xvals.size(), order + 1);
	for (int i = 0; i < xvals.size(); i++) {
		A(i, 0) = 1.0;
	}
	for (int j = 0; j < xvals.size(); j++) {
		for (int i = 0; i < order; i++) {
			A(j, i + 1) = A(j, i) * xvals(j);
		}
	}
	auto Q = A.householderQr();
	auto result = Q.solve(yvals);
	return result;
}

#include <iostream>
#include <vector>
int main(int argc, char** argv)
{
	Eigen::VectorXd xvals(6);
	Eigen::VectorXd yvals(6);
	/*xvals << 9.261977, -2.06803, -19.6663, -36.868, -51.6263, -66.3482;
	yvals << 5.17, -2.25, -15.306, -29.46, -42.85, -57.6116;*/
	Eigen::VectorXf xvalfs(6);
	Eigen::VectorXf yvalfs(6);
	xvalfs << 9.261977, -2.06803, -19.6663, -36.868, -51.6263, -66.3482;
	yvalfs << 18.52, -4.136, -39.333, -73.736, -103.2526, -132.6962;
	xvals = xvalfs.cast<double>();
	yvals = yvalfs.cast<double>();
	std::vector<float> a = { 1, 2, 3, 4 };
	Eigen::VectorXf b = Eigen::Map<Eigen::VectorXf, Eigen::Unaligned>(a.data(), a.size());
	std::cout << "b:\n"<< b << std::endl;
	Eigen::VectorXd d = b.cast<double>();
	std::cout << "d:\n" << d << std::endl;
	auto coeffs = polyfit(xvals, yvals, 1);
	std::cout << coeffs << std::endl;
	std::cout << "Y(16)=" << polyeval(coeffs, 16) << std::endl;
	return 0;
}

等有时间我自己写一个补充上。

 

 

 

 

 

 

参考：C++实现多项式曲线拟合--polyfit

参考：多项式拟合实现（C++）