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

2024-02-01 19:04

本文主要是介绍曲线拟合、多项式拟合、最小二乘法,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

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

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

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

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

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

请添加图片描述

早在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

这篇关于曲线拟合、多项式拟合、最小二乘法的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/668176

相关文章

poj 1258 Agri-Net(最小生成树模板代码)

感觉用这题来当模板更适合。 题意就是给你邻接矩阵求最小生成树啦。~ prim代码:效率很高。172k...0ms。 #include<stdio.h>#include<algorithm>using namespace std;const int MaxN = 101;const int INF = 0x3f3f3f3f;int g[MaxN][MaxN];int n

poj 1287 Networking(prim or kruscal最小生成树)

题意给你点与点间距离,求最小生成树。 注意点是,两点之间可能有不同的路,输入的时候选择最小的,和之前有道最短路WA的题目类似。 prim代码: #include<stdio.h>const int MaxN = 51;const int INF = 0x3f3f3f3f;int g[MaxN][MaxN];int P;int prim(){bool vis[MaxN];

poj 2349 Arctic Network uva 10369(prim or kruscal最小生成树)

题目很麻烦,因为不熟悉最小生成树的算法调试了好久。 感觉网上的题目解释都没说得很清楚,不适合新手。自己写一个。 题意:给你点的坐标,然后两点间可以有两种方式来通信:第一种是卫星通信,第二种是无线电通信。 卫星通信:任何两个有卫星频道的点间都可以直接建立连接,与点间的距离无关; 无线电通信:两个点之间的距离不能超过D,无线电收发器的功率越大,D越大,越昂贵。 计算无线电收发器D

poj 1734 (floyd求最小环并打印路径)

题意: 求图中的一个最小环,并打印路径。 解析: ans 保存最小环长度。 一直wa,最后终于找到原因,inf开太大爆掉了。。。 虽然0x3f3f3f3f用memset好用,但是还是有局限性。 代码: #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#incl

hdu 1102 uva 10397(最小生成树prim)

hdu 1102: 题意: 给一个邻接矩阵,给一些村庄间已经修的路,问最小生成树。 解析: 把已经修的路的权值改为0,套个prim()。 注意prim 最外层循坏为n-1。 代码: #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstri

poj 2175 最小费用最大流TLE

题意: 一条街上有n个大楼,坐标为xi,yi,bi个人在里面工作。 然后防空洞的坐标为pj,qj,可以容纳cj个人。 从大楼i中的人到防空洞j去避难所需的时间为 abs(xi - pi) + (yi - qi) + 1。 现在设计了一个避难计划,指定从大楼i到防空洞j避难的人数 eij。 判断如果按照原计划进行,所有人避难所用的时间总和是不是最小的。 若是,输出“OPETIMAL",若

poj 2135 有流量限制的最小费用最大流

题意: 农场里有n块地,其中约翰的家在1号地,二n号地有个很大的仓库。 农场有M条道路(双向),道路i连接着ai号地和bi号地,长度为ci。 约翰希望按照从家里出发,经过若干块地后到达仓库,然后再返回家中的顺序带朋友参观。 如果要求往返不能经过同一条路两次,求参观路线总长度的最小值。 解析: 如果只考虑去或者回的情况,问题只不过是无向图中两点之间的最短路问题。 但是现在要去要回

poj 3422 有流量限制的最小费用流 反用求最大 + 拆点

题意: 给一个n*n(50 * 50) 的数字迷宫,从左上点开始走,走到右下点。 每次只能往右移一格,或者往下移一格。 每个格子,第一次到达时可以获得格子对应的数字作为奖励,再次到达则没有奖励。 问走k次这个迷宫,最大能获得多少奖励。 解析: 拆点,拿样例来说明: 3 2 1 2 3 0 2 1 1 4 2 3*3的数字迷宫,走两次最大能获得多少奖励。 将每个点拆成两个

poj 2195 bfs+有流量限制的最小费用流

题意: 给一张n * m(100 * 100)的图,图中” . " 代表空地, “ M ” 代表人, “ H ” 代表家。 现在,要你安排每个人从他所在的地方移动到家里,每移动一格的消耗是1,求最小的消耗。 人可以移动到家的那一格但是不进去。 解析: 先用bfs搞出每个M与每个H的距离。 然后就是网络流的建图过程了,先抽象出源点s和汇点t。 令源点与每个人相连,容量为1,费用为

poj 3068 有流量限制的最小费用网络流

题意: m条有向边连接了n个仓库,每条边都有一定费用。 将两种危险品从0运到n-1,除了起点和终点外,危险品不能放在一起,也不能走相同的路径。 求最小的费用是多少。 解析: 抽象出一个源点s一个汇点t,源点与0相连,费用为0,容量为2。 汇点与n - 1相连,费用为0,容量为2。 每条边之间也相连,费用为每条边的费用,容量为1。 建图完毕之后,求一条流量为2的最小费用流就行了