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

最近在做自车轨迹预测的工作，遇到曲线拟合、多项式拟合、最小二乘法这些概念有点不清晰，
做一些概念区别的总结：

曲线拟合用于查找一系列数据点的“最佳拟合”线或曲线。大多数情况下，曲线拟合将产生一个函数，可用于在曲线的任何位置找到点。在某些情况下，也可能不关心找到函数，而是只想使用曲线拟合来平滑数据并改善绘图的外观。

简而言之，曲线拟合就是在受到一定约束条件的情况下，构建可以表示一系列数据点的曲线或数学函数的过程。更加数学化一点的话，对于一个给定的数据集，曲线拟合是以方程形式引入因变量和自变量之间的数学关系的过程。

多项式拟合是假设函数模型是多项式，是曲线拟合中的夹着数据符合的函数模型，还有其他曲线模型拟合，如线性回归，指数型函数，其他非线性曲线拟合；

曲线拟合结果可以有很多，也就是说解有很多，拟合结果的好坏如何评价，如果找到最优解，最优解的评判标准是什么。那么“最佳”的准则是什么？可以是所有观测点到直线的距离和最小，也可以是所有观测点到直线的误差（真实值-理论值）绝对值和最小，也可以是其它。

请添加图片描述

早在19世纪，勒让德就认为让“误差的平方和最小”估计出来的模型是最接近真实情形的。这就是最小二乘法的思想，所谓“二乘”就是平方的意思。从这里我们可以看到，所以最小二乘法其实就是用来做函数拟合的一种思想或者准则。至于怎么求出具体的参数那就是另外一个问题了，理论上可以用导数法、几何法，工程上可以用梯度下降法。

因此最小二乘法（又称最小平方法）是一种数学优化思想。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。

参考下方 C 代码改写实现接口：
实现的 2阶、 3阶最小二乘多项式拟合函数 C 语言接口：

2阶多项式
void QuadraticSpline(Point_xy *points, uint32 n, float32 *coeffs)
{coeffs[3] = 0;// 计算矩阵A和向量bfloat32 sum_x = 0.0, sum_x2 = 0.0, sum_y = 0.0, sum_xy = 0.0, sum_x3 = 0.0, sum_x4 = 0.0, sum_x2y = 0.0;for (uint32 i = 0; i < n; i++){float32 xi = points[i].x, yi = points[i].y;sum_x += xi;sum_x2 += xi * xi;sum_y += yi;sum_xy += xi * yi;sum_x3 += xi * xi * xi;sum_x4 += xi * xi * xi * xi;sum_x2y += xi * xi * yi;}// 构造矩阵A和向量bfloat32 A[3][3] = {{n, sum_x, sum_x2},{sum_x, sum_x2, sum_x3},{sum_x2, sum_x3, sum_x4}};float32 b[3] = {sum_y, sum_xy, sum_x2y};// 求解线性方程组Ax=bfor (uint8 k = 0; k < 3 - 1; k++){for (uint8 i = k + 1; i < 3; i++){float32 p = A[i][k] / A[k][k];for (uint8 j = k + 1; j < 3; j++){A[i][j] -= p * A[k][j];}b[i] -= p * b[k];}}coeffs[2] = b[2] / A[2][2];coeffs[1] = (b[1] - A[1][2] * coeffs[2]) / A[1][1];coeffs[0] = (b[0] - A[0][1] * coeffs[1] - A[0][2] * coeffs[2]) / A[0][0];
}

3阶多项式
void CubicSpline(Point_xy* points, uint32 n, float32* coeffs)
{coeffs[3] = 0;// 计算矩阵A和向量bfloat32 sum_x = 0.0, sum_x2 = 0.0, sum_x3 = 0.0, sum_x4 = 0.0, sum_x5 = 0.0, sum_x6 = 0.0, sum_y = 0.0, sum_xy = 0.0, sum_x2y = 0.0, sum_x3y = 0.0;for (uint32 i = 0; i < n; i++){float32 xi = points[i].x, yi = points[i].y;sum_x  += xi;sum_x2 += xi * xi;sum_x3 += xi * xi * xi;sum_x4 += xi * xi * xi * xi;sum_x5 += xi * xi * xi * xi * xi;sum_x6 += xi * xi * xi * xi * xi * xi;sum_y  += yi;sum_xy += xi * yi;sum_x2y += xi * xi * yi;sum_x3y += xi * xi * xi * yi;}// 构造矩阵A和向量bfloat32 A[4][4] = { {n, sum_x, sum_x2, sum_x3},{sum_x, sum_x2, sum_x3, sum_x4},{sum_x2, sum_x3, sum_x4,sum_x5},{sum_x3, sum_x4, sum_x5,sum_x6}};float32 b[4] = { sum_y, sum_xy, sum_x2y,sum_x3y};// 求解线性方程组Ax=bfor (uint8 k = 0; k < 4 - 1; k++){for (uint8 i = k + 1; i < 4; i++){float32 p = A[i][k] / A[k][k];for (uint8 j = k + 1; j < 4; j++){A[i][j] -= p * A[k][j];}b[i] -= p * b[k];}}coeffs[3] = b[3] / A[3][3];coeffs[2] = (b[2] - A[2][3] * coeffs[3]) / A[2][2];coeffs[1] = (b[1] - A[1][2] * coeffs[2] - A[1][3] * coeffs[3]) / A[1][1];coeffs[0] = (b[0] - A[0][1] * coeffs[1] - A[0][2] * coeffs[2] - A[0][3] * coeffs[3]) / A[0][0];
}

C ++ 实现：
代码出处

#include <stdio.h>
#include "stdlib.h"
#include "math.h"
#include <vector>
using namespace std;struct point
{double x;double y;
};typedef vector<double> doubleVector;
vector<point> getFile(char *File);  //获取文件数据
doubleVector getCoeff(vector<point> sample, int n);   //矩阵方程void main()
{int i, n;char *File = "XY.txt";vector<point> sample;doubleVector  coefficient;sample = getFile(File);printf("拟合多项式阶数n=");scanf_s("%d", &n);coefficient = getCoeff(sample, n);printf("\n拟合矩阵的系数为：\n");for (i = 0; i < coefficient.size(); i++)printf("a%d = %lf\n", i, coefficient[i]);}
//矩阵方程
doubleVector getCoeff(vector<point> sample, int n)
{vector<doubleVector> matFunX;  //公式3左矩阵vector<doubleVector> matFunY;  //公式3右矩阵doubleVector temp;double sum;int i, j, k;//公式3左矩阵for (i = 0; i <= n; i++){temp.clear();for (j = 0; j <= n; j++){sum = 0;for (k = 0; k < sample.size(); k++)sum += pow(sample[k].x, j + i);temp.push_back(sum);}matFunX.push_back(temp);}//打印matFunX矩阵printf("matFunX矩阵：\n");for (i = 0;i < matFunX.size();i++) {for (j = 0;j < matFunX.size();j++)printf("%f\t", matFunX[i][j]);printf("\n");}printf("matFunX.size=%d\n", matFunX.size());//printf("matFunX[3][3]=%f\n", matFunX[3][3]);//公式3右矩阵for (i = 0; i <= n; i++){temp.clear();sum = 0;for (k = 0; k < sample.size(); k++)sum += sample[k].y*pow(sample[k].x, i);temp.push_back(sum);matFunY.push_back(temp);}printf("matFunY.size=%d\n", matFunY.size());//打印matFunY的矩阵printf("matFunY的矩阵：\n");for (i = 0;i < matFunY.size();i++) {printf("%f\n", matFunY[i][0]);}//矩阵行列式变换，将matFunX矩阵变为下三角矩阵,将matFunY矩阵做怎样的变换呢？//AX=Y在将X变换为下三角矩阵X'时，是左右两边同乘ratiodouble num1, num2, ratio;for (i = 0; i < matFunX.size() - 1; i++){num1 = matFunX[i][i];for (j = i + 1; j < matFunX.size(); j++){num2 = matFunX[j][i];ratio = num2 / num1;for (k = 0; k < matFunX.size(); k++)matFunX[j][k] = matFunX[j][k] - matFunX[i][k] * ratio;matFunY[j][0] = matFunY[j][0] - matFunY[i][0] * ratio;}}//打印matFunX行列变化之后的矩阵printf("matFunX行列变换之后的矩阵：\n");for (i = 0;i < matFunX.size();i++) {for (j = 0;j < matFunX.size();j++)printf("%f\t", matFunX[i][j]);printf("\n");}//打印matFunY行列变换之后的矩阵printf("matFunY行列变换之后的矩阵：\n");for (i = 0;i < matFunY.size();i++) {printf("%f\n", matFunY[i][0]);}//计算拟合曲线的系数doubleVector coeff(matFunX.size(), 0);for (i = matFunX.size() - 1; i >= 0; i--){if (i == matFunX.size() - 1)coeff[i] = matFunY[i][0] / matFunX[i][i];else{for (j = i + 1; j < matFunX.size(); j++)matFunY[i][0] = matFunY[i][0] - coeff[j] * matFunX[i][j];coeff[i] = matFunY[i][0] / matFunX[i][i];}}return coeff;
}//获取文件数据
vector<point> getFile(char *File)
{int i = 1;vector<point> dst;FILE *fp = fopen(File, "r");if (fp == NULL){printf("Open file error!!!\n");exit(0);}point temp;double num;while (fscanf(fp, "%lf", &num) != EOF){if (i % 2 == 0){temp.y = num;dst.push_back(temp);}elsetemp.x = num;i++;}fclose(fp);return dst;
}

改成C 语言：
代码出处

#include <stdio.h>
#include "stdlib.h"
#include "math.h"
//#include <vector>
//using namespace std;#define MAX_EXP 4 /* x最高次幂 */
#define MATRIX_DIM (MAX_EXP + 1)   /* 矩阵阶数 */
#define SMPNUM 5    /* 采样个数 */struct point
{double x;double y;
};/* 采样点想， y */
float sampleX[SMPNUM] = {0.0};
float sampleY[SMPNUM] = {0.0};void getFile(char *File);  //获取文件数据
void getCoeff(float sampleX[SMPNUM], float sampleY[SMPNUM], float coeff[MATRIX_DIM]);   //获取系数int main()
{int i;char *File = "XY.txt";//vector<point> sample;//doubleVector  coefficient;//sample = getFile(File);getFile(File);printf("拟合多项式阶数n=");//scanf("%d", &n);// coefficient = getCoeff(sample, n);//n = 3;float coeff[MATRIX_DIM] = {0};getCoeff(sampleX, sampleY, coeff);printf("\n拟合矩阵的系数为：\n");for (i = 0; i < MATRIX_DIM; i++){printf("a%d = %lf\n", i, coeff[i]);}return 0;
}
//矩阵方程
void getCoeff(float sampleX[SMPNUM], float sampleY[SMPNUM], float coeff[MATRIX_DIM])
{double sum;int i, j, k;float matX[MATRIX_DIM][MATRIX_DIM];  //公式3左矩阵float matY[MATRIX_DIM][MATRIX_DIM];  //公式3右矩阵float temp2[MATRIX_DIM];//公式3左矩阵for (i = 0; i < MATRIX_DIM; i++){for (j = 0; j < MATRIX_DIM; j++){sum = 0.0;for (k = 0; k < SMPNUM; k++){sum += pow(sampleX[k], j + i);}matX[i][j] = sum;}}//打印matFunX矩阵printf("matFunX矩阵：\n");for (i = 0; i < MATRIX_DIM; i++){for (j = 0; j < MATRIX_DIM; j++){printf("%f\t", matX[i][j]);}printf("\n");}//公式3右矩阵for (i = 0; i < MATRIX_DIM; i++){//temp.clear();sum = 0;for (k = 0; k < SMPNUM; k++){sum += sampleY[k] * pow(sampleX[k], i);}matY[i][0] = sum;}//printf("matFunY.size=%d\n", matFunY.size());//打印matFunY的矩阵printf("matFunY的矩阵：\n");for (i = 0; i < MATRIX_DIM; i++){printf("%f\n", matY[i][0]);}//矩阵行列式变换，将matFunX矩阵变为下三角矩阵,将matFunY矩阵做怎样的变换呢？//AX=Y在将X变换为下三角矩阵X'时，是左右两边同乘ratiodouble num1, num2, ratio;for (i = 0; i < MATRIX_DIM; i++){num1 = matX[i][i];for (j = i + 1; j < MATRIX_DIM; j++){num2 = matX[j][i];ratio = num2 / num1;for (k = 0; k < MATRIX_DIM; k++){matX[j][k] = matX[j][k] - matX[i][k] * ratio;}matY[j][0] = matY[j][0] - matY[i][0] * ratio;}}//打印matFunX行列变化之后的矩阵printf("matFunX行列变换之后的矩阵：\n");for (i = 0; i < MATRIX_DIM; i++){for (j = 0; j < MATRIX_DIM; j++){printf("%f\t", matX[i][j]);}printf("\n");}//打印matFunY行列变换之后的矩阵printf("matFunY行列变换之后的矩阵：\n");for (i = 0; i < MATRIX_DIM; i++){printf("%f\n", matY[i][0]);}//计算拟合曲线的系数//doubleVector coeff(n + 1, 0);for (i = MATRIX_DIM - 1; i >= 0; i--){if (i == MATRIX_DIM - 1){coeff[i] = matY[i][0] / matX[i][i];}else{for (j = i + 1; j < MATRIX_DIM; j++){matY[i][0] = matY[i][0] - coeff[j] * matX[i][j];}coeff[i] = matY[i][0] / matX[i][i];}}return;
}//获取文件数据
void getFile(char *File)
{int i = 0, j = 0, k = 0;//vector<point> dst;FILE *fp = fopen(File, "r");if (fp == NULL){printf("Open file error!!!\n");exit(0);}point temp;double num;while (fscanf(fp, "%lf", &num) != EOF){if (i % 2 == 0){temp.x = num;sampleX[j++] = num;}else{temp.y = num;sampleY[k++] = num;}i++;}fclose(fp);return;//return dst;
}

C语言消元法求解代码实现：
代码出处

#define RANK_  	3  		多项式次数
/*
*********************************************************************************************************
*   函 数 名: polyfit
*   功能说明: 多项式曲线拟合（与matlab效果一致）
*   形    参: d_X	输入的数据的x值d_Y	输入的数据的y值d_N	输入的数据的个数rank  多项式的次数coeff 输出的系数
*   返 回 值: 无
*********************************************************************************************************
*/
//原理：At * A * C = At * Y	，其中 At 为 A 转置矩阵，C 为系数矩阵
void polyfit(double *d_X, double *d_Y, int d_N, int rank, double *coeff)
{if(RANK_ != rank)	//判断次数是否合法return;int i,j,k;	double aT_A[RANK_ + 1][RANK_ + 1] = {0};double aT_Y[RANK_ + 1] = {0};for(i = 0; i < rank + 1; i++)	//行{for(j = 0; j < rank + 1; j++)	//列{for(k = 0; k < d_N; k++)	{aT_A[i][j] +=  pow(d_X[k], i+j);		//At * A 线性矩阵}}}for(i = 0; i < rank + 1; i++){for(k = 0; k < d_N; k++){aT_Y[i] += pow(d_X[k], i) * d_Y[k];		//At * Y 线性矩阵}}//以下为高斯列主元素消去法解线性方程组for(k = 0; k < rank + 1 - 1; k++){int row = k;double mainElement = fabs(aT_A[k][k]);double temp = 0.0;//找主元素for(i = k + 1; i < rank + 1 - 1; i++){if( fabs(aT_A[i][i]) > mainElement ){mainElement = fabs(aT_A[i][i]);row = i;}}//交换两行if(row != k){for(i = 0; i < rank + 1; i++){temp = aT_A[k][i];aT_A[k][i] = aT_A[row][i];aT_A[row][i] = temp;}	temp = aT_Y[k];aT_Y[k] = aT_Y[row];aT_Y[row] = temp;}//消元过程for(i = k + 1; i < rank + 1; i++){for(j = k + 1; j < rank + 1; j++){aT_A[i][j] -= aT_A[k][j] * aT_A[i][k] / aT_A[k][k];}aT_Y[i] -= aT_Y[k] * aT_A[i][k] / aT_A[k][k];}}	//回代过程for(i = rank + 1 - 1; i >= 0; coeff[i] /= aT_A[i][i], i--){for(j = i + 1, coeff[i] = aT_Y[i]; j < rank + 1; j++){coeff[i] -= aT_A[i][j] * coeff[j];}}return;	
}

python 代码实现多项式曲线拟合：
代码出处

'''多项式：yi = w0 + w1*xi^1 + w2*xi^2 + ... + wn*xi^mN为数据点个数，M为阶数先用数据点（xa、ya）求出未知参数，然后用参数带入后的公式求解给定值(xxa)
'''
import matplotlib.pyplot as plt
import numpy
import randomfig = plt.figure()
ax = fig.add_subplot(111)# 在这里给出拟合多项式的阶数
order = 9# 1. 生成曲线上的各个点
x = numpy.arange(-1,1,0.02)
y = [((a*a-1)*(a*a-1)*(a*a-1)+0.5)*numpy.sin(a*2) for a in x]   # y=[(x^2-1)^3]*sin(2x)，阶数为6？？？# ax.plot(x,y,color='r',linestyle='-',marker='')
# ,label="(a*a-1)*(a*a-1)*(a*a-1)+0.5"
plt.plot(x,y,c='red')# 2. 生成的曲线上的各个点偏移一下，并放入到x_data,y_data中去
i=0
x_data=[]
y_data=[]
for xx in x:yy=y[i]d=float(random.randint(60,140))/100#ax.plot([xx*d],[yy*d],color='m',linestyle='',marker='.')i+=1x_data.append(xx*d)y_data.append(yy*d)ax.plot(x_data,y_data,color='m',linestyle='',marker='.')# 3. 计算Ax=b中，矩阵A、b
# 存储从0次到m次的所有冥方和
bigMat=[]
for j in range(0, 2*order+1):sum = 0for i in range(0,len(xa)):sum += (xa[i]**j)bigMat.append(sum)# 计算线性方程组系数矩阵：A
matA=[]
for rowNum in range(0,order+1):row=bigMat[rowNum:rowNum+order+1]matA.append(row)matA=numpy.array(matA)matB=[]
for i in range(0,order+1):ty=0.0for k in range(0,len(xa)):ty+=ya[k]*(xa[k]**i)matB.append(ty)matB=numpy.array(matB)matW=numpy.linalg.solve(matA,matB)
# numpy.linalg中的函数solve可以求解形如 Ax = b 的线性方程组，其中 A 为矩阵，b 为一维或二维的数组，x 是未知变量# 画出拟合后的曲线
# print(matW)
x_ = numpy.arange(-1,1.06,0.01)
y_ =[]
for i in range(0,len(xxa)):yy=0.0for j in range(0,order+1):dy = (x_[i]**j)*matW[j]# dy*=matW[j]yy += dyy_.append(yy)
ax.plot(x_,y_,color='g',linestyle='-',marker='')ax.legend()
plt.show()

相关参考：
https://blog.csdn.net/qq_27586341/article/details/90170839
https://blog.csdn.net/MoreAction_/article/details/106443383
https://www.cnblogs.com/nhyq-wyj/p/14898517.html
https://sikasjc.github.io/2018/10/24/curvefitting/

https://blog.csdn.net/piaoxuezhong/article/details/54973750
https://blog.csdn.net/tutu1583/article/details/82111060