参考文章：马同学马同学提供线性代数,微积分,概率论与数理统计,机器学习等知识讲解https://www.matongxue.com/madocs/818/
C++最小二乘法拟合-（线性拟合和多项式拟合）_尘中远的博客-CSDN博客_namespace gsl
最小二乘法—多项式拟合非线性函数 - 简书
最小二乘法，即：最小化误差的平方和：
其中，为样本点的数量，
为预测值，
为当前值。
-
直线拟合：线性拟合（一元函数）
如果想要根据输入数据拟合一条直线，那么可以表示为：
则为：
当a、b偏导数为0时，最小：
展开累加符号得：
根据克拉默（克莱姆）法则求解线性方程组（省略累加符号的上下限）：
C++代码
template<typename T>
bool linearFit(const T* x, const T* y,size_t length)
{
factor.resize(2,0);
typename T t1=0, t2=0, t3=0, t4=0;
for(int i=0; i<length; ++i)
{
t1 += x[i]*x[i];
t2 += x[i];
t3 += x[i]*y[i];
t4 += y[i];
}
factor[1] = (t3*length - t2*t4) / (t1*length - t2*t2);
factor[0] = (t1*t4 - t2*t3) / (t1*length - t2*t2);
return true;
}-
曲线拟合：多项式拟合（多元函数）
如果想要根据输入数据拟合一条曲线，那么可以表示为：
则为：
当、
、
、...
偏导数为0时，
最小：
展开累加符号得（省略累加符号的上下限）：
写成矩阵相乘的形式：
其中
A、b可根据样本数据集计算得到，然后使用高斯列主元消元法解线性方程组，得到x
计算A、b的C++代码
void polyfit(const T* x,const T* y,size_t length,int poly_n)
{
//将最优化偏导数的方程组写为矩阵乘法形式：Ax=b，其中A、b可根据待拟合的点集坐标计算得到，x为最终待求的多项式系数矩阵
//利用高斯列主元消元法来求解Ax=b
//3阶多项式f(x)=ax^3 + bx^2 + cx^1 + dx^0,可见n阶多项式，具有n+1个待求系数，需要n+1个方程组才能解出
std::vector<double> factor;
factor.resize(poly_n+1,0);
int i,j;
//求x加和的相关项（系数矩阵中有重复项）
std::vector<double> tempx(length, 1.0);
std::vector<double> sumxx(poly_n * 2 + 1);
for (i=0;i<2*poly_n+1;i++){
for (sumxx[i]=0,j=0;j<length;j++)
{
sumxx[i]+=tempx[j];
tempx[j]*=x[j];
}
}
//根据重复项规律求matrix_A
std::vector<double> matrix_A((poly_n + 1)*(poly_n + 1));
for (i=0;i<poly_n+1;i++)
for (j=0;j<poly_n+1;j++)
matrix_A[i*(poly_n+1)+j]=sumxx[i+j];
//求y加和的相关项matrix_b
std::vector<double> tempy(y, y + length);//通过y中0到length-1个元素初始化向量
std::vector<double> matrix_b(poly_n + 1);
for (i = 0; i < poly_n + 1; i++) {
for (matrix_b[i] = 0, j = 0; j < length; j++)
{
matrix_b[i] += tempy[j];
tempy[j] *= x[j];
}
}
//求多项式系数x，用变量factor代替
gauss_solve(poly_n+1, matrix_A, factor, matrix_b);
}高斯列主元消元法计算x的c++代码
template<typename T>
void gauss_solve(int n
, std::vector<typename T>& A
, std::vector<typename T>& x
, std::vector<typename T>& b)
{
GaussElimination(n, &A[0], &x[0], &b[0]);
}
template<typename T>
void GaussElimination(int n, T* A, T* x, T* b)
{
int i, j, k; //索引
int r; //主元所在行数
double maxValue; //当前列中绝对值最大的值，即列主元
double factor; //消元系数
double temp; //临时变量
//逐列进行消元
for (k = 0; k < n - 1; k++)
{
//获取列主元所在行数r
r = k;
maxValue = fabs(A[k*n + k]);
for (i = k + 1; i < n; i++)
{
if (maxValue < fabs(A[i*n + k]))
{
maxValue = fabs(A[i*n + k]);
r = i;
}
}
//交换行元素，当r == k时，不需要交换
if (r != k)
{
//交换等号左侧系数矩阵的行元素
for (i = 0; i < n; i++)
{
temp = A[k*n + i];
A[k*n + i] = A[r*n + i];
A[r*n + i] = temp;
}
//交换等号右侧向量的行元素
temp = b[k];
b[k] = b[r];
b[r] = temp;
}
//cout << "第" << k + 1 << "次选取主元并交换行的顺序" << endl;
//DispalyMatrix(n, (double*)A, b);
//按列消元过程
for (i = k + 1; i < n; i++)
{
factor = A[i*n + k] / A[k*n + k];
for (j = k+1; j < n; j++) //j = k时更便于观察结果，但是出现冗余计算
{
A[i*n + j] = A[i*n + j] - factor * A[k*n + j];
}
b[i] = b[i] - factor * b[k];
}
//cout << "第" << k + 1 << "次消元" << endl;
//DispalyMatrix(n, (double*)A, b);
}
//判断系数矩阵的奇异性
if (A[k*n + k] < 0.0001)
{
//cout << "系数矩阵是非奇异矩阵,方程无唯一解!" << "\n";
}
//回代法，将消元后的上三角矩阵变为单位矩阵，b即为系数
for (i = n - 1; i >= 0; i--)
{
temp = 0;
for (j = i + 1; j < n; j++)
{
temp = temp + A[i*n + j] * b[j];
}
b[i] = (b[i] - temp) / A[i*n + i];
x[i] = b[i];
}
//cout << "回代计算后的矩阵:" << "\n";
//for (i = 0; i < n; i++)
//{
// for (j = 0; j < n; j++)
// {
// if (i == j)
// A[i*n + j] = 1;
// else
// A[i*n + j] = 0;
// }
//}
//DispalyMatrix(n, (double*)A, b);
//cout << "最后得方程的根为:" << "\n";
//for (i = 0; i < n; i++)
// cout << "x" << i + 1 << "=" << b[i] << "\n";
}
//void DispalyMatrix(int n, const double *A, const double *b)
//{
// int i, j;
// for (i = 0; i < n; i++)
// {
// for (j = 0; j < n; j++)
// cout << setw(4) << setprecision(4) << *(A + i * n + j) << " ";
// cout << b[i] << endl;
// }
// cout << endl;
//}
