在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++)