学校网站建设的意义和应用,服务平台管理系统,奥林匹克做校服的网站,网站建设需要提供那些资料最近在做自车轨迹预测的工作#xff0c;遇到 曲线拟合、多项式拟合、最小二乘法这些概念有点不清晰#xff0c; 做一些概念区别的总结#xff1a;
曲线拟合用于查找一系列数据点的“最佳拟合”线或曲线。 大多数情况下#xff0c;曲线拟合将产生一个函数#xff0c;可用于…最近在做自车轨迹预测的工作遇到 曲线拟合、多项式拟合、最小二乘法这些概念有点不清晰 做一些概念区别的总结
曲线拟合用于查找一系列数据点的“最佳拟合”线或曲线。 大多数情况下曲线拟合将产生一个函数可用于在曲线的任何位置找到点。 在某些情况下也可能不关心找到函数而是只想使用曲线拟合来平滑数据并改善绘图的外观。
简而言之曲线拟合就是在受到一定约束条件的情况下构建可以表示一系列数据点的曲线或数学函数的过程。更加数学化一点的话对于一个给定的数据集曲线拟合是以方程形式引入因变量和自变量之间的数学关系的过程。
多项式拟合是假设函数模型是多项式是曲线拟合中的夹着数据符合的函数模型 还有其他曲线模型拟合如线性回归指数型函数其他非线性曲线拟合
曲线拟合结果可以有很多也就是说解有很多拟合结果的好坏如何评价如果找到最优解 最优解的评判标准是什么。那么“最佳”的准则是什么可以是所有观测点到直线的距离和最小也可以是所有观测点到直线的误差真实值-理论值绝对值和最小也可以是其它。 早在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};// 求解线性方程组Axbfor (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};// 求解线性方程组Axbfor (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 vectordouble doubleVector;
vectorpoint getFile(char *File); //获取文件数据
doubleVector getCoeff(vectorpoint sample, int n); //矩阵方程void main()
{int i, n;char *File XY.txt;vectorpoint 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(vectorpoint sample, int n)
{vectordoubleVector matFunX; //公式3左矩阵vectordoubleVector 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矩阵做怎样的变换呢//AXY在将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;
}//获取文件数据
vectorpoint getFile(char *File)
{int i 1;vectorpoint 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;//vectorpoint 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矩阵做怎样的变换呢//AXY在将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;//vectorpoint 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], ij); //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,colorr,linestyle-,marker)
# ,label(a*a-1)*(a*a-1)*(a*a-1)0.5
plt.plot(x,y,cred)# 2. 生成的曲线上的各个点偏移一下并放入到x_data,y_data中去
i0
x_data[]
y_data[]
for xx in x:yyy[i]dfloat(random.randint(60,140))/100#ax.plot([xx*d],[yy*d],colorm,linestyle,marker.)i1x_data.append(xx*d)y_data.append(yy*d)ax.plot(x_data,y_data,colorm,linestyle,marker.)# 3. 计算Axb中矩阵A、b
# 存储从0次到m次的所有冥方和
bigMat[]
for j in range(0, 2*order1):sum 0for i in range(0,len(xa)):sum (xa[i]**j)bigMat.append(sum)# 计算线性方程组系数矩阵A
matA[]
for rowNum in range(0,order1):rowbigMat[rowNum:rowNumorder1]matA.append(row)matAnumpy.array(matA)matB[]
for i in range(0,order1):ty0.0for k in range(0,len(xa)):tyya[k]*(xa[k]**i)matB.append(ty)matBnumpy.array(matB)matWnumpy.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)):yy0.0for j in range(0,order1):dy (x_[i]**j)*matW[j]# dy*matW[j]yy dyy_.append(yy)
ax.plot(x_,y_,colorg,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
