最小二乘法拟合多项式曲线(含C++代码)

原创 2018年04月15日 22:58:03

本文原创,转载请注明来源:https://blog.csdn.net/syc666946/article/details/79954587

最小二乘原理见https://blog.csdn.net/jairuschan/article/details/7517773

关于原文下面评论中的问题:为什么最后一步明明是求解线性方程组AX=Y(X, Y均为向量)中的系数矩阵A,

而A的表达式为A = (X'*X)-1*X'*Y,这是一个超定方程组,A矩阵是奇异矩阵,所以A并不可直接求逆,只能求伪逆。

C++代码如下:https://github.com/shaoyuncen/MyPCL/blob/master/LeastSquare/LeastSquare/class.h

矩阵计算部分用了Eigen库,其中的输入vector中的数据类型MyPoint是我自己定义的X,Y,Z点类型,具体见github。

compute方法输入为同一平面上的点,用于估计系数矩阵A; 具体拟合到几次方可以设置_POWER

#include <iostream>
#include <vector>
#include <math.h>
#include <fstream>
#include <sstream>
#include "Eigen/Eigen"
#include <Eigen/SVD>
int _POWER; //参数个数, _POWER-1等于拟合的最高次数
class LeastSquare {
private:
	Eigen::RowVectorXd coeficients;//系数矩阵
	double Deviation;
	double k, b;

public:
	LeastSquare(const vector<MyPoint> &ptr) {//ptr是一条线的点云数据
		coeficients.setZero();//初始化为0矩阵
		compute(ptr);//计算估计参数矩阵,拟合曲线;
		//compute1(ptr);
	}
	double getk() { return k; }
	double getb() { return b; }
	double power(const double number, int power) 
	{
		double temp = 1;
		while (power--)
		{
			temp *= number;
		}
		return temp;
	}
	void compute1(const vector<MyPoint> &ptr);//compute the k and b;
	void compute(const vector<MyPoint> &ptr);////parameter estimate---compute coeficients
	void filter(const vector<MyPoint> &ptr, vector<MyPoint> &outlier);
	void fucking_print()
	{
		//cout << "y = " << k << "x + " << b << "\t";
		cout << "曲线的多项式系数(形如y = a0 + a1*x + a2*x2 +...am*xm)为:\n" << coeficients << endl;
	}
};

void LeastSquare::compute1(const vector<MyPoint> &ptr) //compute k and b, 这个是之前写的拟合直线的,下面的多项式拟合次数取2和这个方法结果一致
{
	double SumXi = 0.0;
	double SumYi = 0.0;
	double SumXi_2 = 0.0;
	double SumXiYi = 0.0;
	double deviation = 0.0;
	int number = int(ptr.size());//size(返回的是无符号整型,需要强制类型转换后使用)
	for (int i = 0; i < number;i++)
	{
		SumXi += ptr[i].x;
		SumYi += ptr[i].z;
		SumXi_2 += ptr[i].x * ptr[i].x;
		SumXiYi += ptr[i].x * ptr[i].z;
	}

	double temp = (number * SumXi_2 - SumXi * SumXi);//减少计算次数
	k = (number * SumXiYi - SumXi * SumYi) / temp;
	b = (SumXi_2 * SumYi - SumXi * SumXiYi) / temp;
}
void LeastSquare::compute(const vector<MyPoint> &ptr) //多项式拟合的参数估计
{
	if (_POWER < 2) return; //参数至少两个
	const int PointNum = int(ptr.size());
	Eigen::RowVectorXd Xi_1(PointNum);//Xi的一次矩阵,1行n列
	Eigen::RowVectorXd Yi_1(PointNum);//Yi的一次矩阵,1行n列
	Eigen::MatrixXd A(PointNum, _POWER);//Ax = b中的A,结构矩阵,x为系数矩阵coeficients,b为数据矩阵,即Yi向量(1 * n的矩阵),power为次数
	Xi_1.setZero();
	Yi_1.setZero();
	A.setOnes();
	for (unsigned int i = 0; i < ptr.size(); i++)//Ps:行的下标对应文件的行数-2,因为读文件的时候第一行会被忽略,然后下标又是从0开始
	{
		Xi_1[i] = ptr[i].x;
		Yi_1[i] = ptr[i].z;
		for (int cols = 0; cols < _POWER; cols++)
		{
			A(i, cols) = power(ptr[i].x, cols);
		}
		//A(i, _POWER - 1) = ptr[i].x;
	}
	Eigen::MatrixXd temp;
	temp = A.transpose()* A;
	coeficients = (temp.inverse() * A.transpose() ) * Yi_1.transpose();
	
	//coeficients = A.ldlt().solve(Yi_1); // A sym. p.s.d.    #include <Eigen/Cholesky>
	//coeficients = A.llt().solve(Yi_1);  // A sym. p.d.      #include <Eigen/Cholesky>
	//coeficients = A.lu().solve(Yi_1);  // Stable and fast. #include <Eigen/LU>
	//coeficients = A.qr().solve(Yi_1);  // No pivoting.     #include <Eigen/QR>
	//coeficients = A.svd().solve(Yi_1);  // Stable, slowest. #include <Eigen/SVD>
	
}

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/syc666946/article/details/79954587

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

最小二乘拟合曲线原理,以及c++详细代码,最后给出了测试用例,将一组数据拟合成二次曲线。...
  • u010418035
  • u010418035
  • 2015-06-30 16:45:28
  • 8909

C++最小二乘法拟合-(线性拟合和多项式拟合)

在进行曲线拟合时用的最多的是最小二乘法,其中以一元函数(线性)和多元函数(多项式)居多,本文介绍的这个类,用C++封装了专门用于进行多项式拟合和线性拟合的方法,可以根据用户输入的阶次进行多项式拟合,算...
  • czyt1988
  • czyt1988
  • 2014-03-23 20:38:13
  • 19616

多项式曲线拟合最小二乘法

多项式插值法
  • zb1165048017
  • zb1165048017
  • 2015-09-10 21:57:01
  • 6049

最小二乘法 多项式曲线拟合 原理公式理解 Python 实现

曲线拟合方法 最小二乘法 最大似然估计 梯度下降法 多项式拟合 Python代码 数据集征兵抽签1-366号y366个不同的人抽x结果表明生日靠后的人易抽到小号概念最小二乘法多项式曲线拟合,根据给定的...
  • neuldp
  • neuldp
  • 2016-07-27 20:54:52
  • 3345

最小二乘法进行最高3次曲线拟合

最近在做跟踪时,需要预测被跟踪物体的运动轨迹,由于被跟踪物体为车辆,轨迹使用二次曲线基本可以较好的拟合,因此做一下实验。 下面为最小二乘法的核心代码,有需要可以参考: bool CNXMinSqu...
  • spacegrass
  • spacegrass
  • 2016-12-20 13:46:34
  • 1878

最小二乘法多项式曲线拟合原理与实现

概念 最小二乘法多项式曲线拟合,根据给定的m个点,并不要求这条曲线精确地经过这些点,而是曲线y=f(x)的近似曲线y= φ(x)。 原理 [原理部分由个人根据互联网上的资料进行总结,希望对大家能有用...
  • ccwwff
  • ccwwff
  • 2013-12-05 10:41:34
  • 1863

基于最小二乘法的曲线拟合的C++代码的实现

简单思路如下: 1,采用目标函数对多项式系数求偏导,得到最优值条件,组成一个方程组; 2,方程组的解法采用行列式变换(两次变换:普通行列式——三角行列式——对角行列式——求解),行列式的求解算法上优化...
  • maweifei
  • maweifei
  • 2016-07-24 17:28:07
  • 4125

最小二乘法曲线拟合(C++)

  • 2010年10月23日 16:49
  • 6KB
  • 下载

C++实现最小二乘法一元回归和多项式拟合

  • 2014年03月23日 20:50
  • 331KB
  • 下载

最小二乘法多项式拟合的Java实现

背景 由于项目中需要根据磁盘的历史使用情况预测未来一段时间的使用情况,决定采用最小二乘法做多项式拟合,这里简单描述下: 假设给定的数据点和其对应的函数值为 (x1, y1), (x2, y2), ...
  • funnyrand
  • funnyrand
  • 2015-07-03 15:47:37
  • 7607
收藏助手
不良信息举报
您举报文章:最小二乘法拟合多项式曲线(含C++代码)
举报原因:
原因补充:

(最多只允许输入30个字)