【C++】线性拟合、指数曲线拟合的最小二乘法作业

曲线拟合的最小二乘法

引入

• 对于未知的函数$y=f(x)$，已知一些节点数据$A_m(x_i,y_i)$，希望构建函数$y=\phi (x)$逼近$y=f(x)$。

• 在插值中，一般严格要求在每个插值节点处都没有偏差，都满足${\delta }_i=\phi (x_i)-f(x_i)=0$；而有时候拟合时允许在节点处有偏差，但总的偏差应该尽可能小。

• 可以使偏差的绝对值之和尽可能小，但为了方便计算，采用偏差的平方和，即$\sum _{i=1}^m[\phi (x_i)-f(x_i){]}^2$尽可能小。

• 按此原则构建拟合曲线$y=\phi (x)$的方法，称为曲线拟合的最小二乘法。

丰富

考虑计算方便和实际意义，进行两处的丰富：

1. $\phi (x)$的函数类型需要提前选取好。
   常常在简单函数中选取，如代数多项式或三角多项式。
   可以选好$n(n<m)$个函数${\phi }_i(x)$，构成$\phi (x)=a_0{\phi }_0(x)+a_1{\phi }_1(x)+...+a_n{\phi }_n(x)$。$a_i$为参数，后续需要求解。

2. 节点数据$A_m$的重要性有时不同。
   有必要添加权系数${\omega }_i>0$，使得$\sum _{i=1}^m{\omega }_i[\phi (x_i)-f(x_i){]}^2$最小。

综上，根据最小二乘问题：$\sum _{i=1}^m{\omega }_i[\phi (x_i)-f(x_i){]}^2$最小，求出的$\phi (x)=a_0{\phi }_0(x)+a_1{\phi }_1(x)+...+a_n{\phi }_n(x)$则是该最小二乘问题的最小二乘解。

求解

${\phi }_i(x)$人为选取后，需要求解$\phi (x)=a_0{\phi }_0(x)+a_1{\phi }_1(x)+...+a_n{\phi }_n(x)$中的参数$a_i$来确定$\phi (x)$的形式，即是求多元函数$S(a_0,a_1,...a_n)=\sum _{i=1}^m{\omega }_i[\sum _{k=0}^n{a}_k{\phi }_k(x_i)-f(x_i){]}^2$的极小值点。

由多元函数求极小值的必要条件，有$\frac{\partial S}{\partial a_j}=0,j=0,1,...,n$

求偏导后，即

移项，即

下面是推导过程，看懂即可，此处不再手打。关于下图中的离散点集的内积，类似函数的内积，也就是两个函数在每个离散点处乘积的累加。由于离散，所以只在节点处累加，如果是函数，则是在定义域积分，可看作无数离散点、无数次累加。
不过，下图中写成内积形式，是为了以简洁的形式写出(3.5.9)的矩阵；不以简洁形式则类似公式(3.5.11)的形式，不过(3.5.11)代入了选取的${\phi }_i(x)$。

选取${\phi }_i(x)=x^i,i=0,1,...,n.$

为了简单地编写下面的作业，只需要利用公式(3.5.11)中的矩阵。

例题3.7：直线拟合

题目介绍

编程之前

对于线性函数的拟合，n=1，只需要取法方程组的前两行两列来计算a0和a1。

然后可以再利用前面一篇中求解线性方程组中的函数来求出$(a_0,a_1)$。

代码

//数据准备：x、y为已知的节点数据，w为权重数据，m为节点数
void PreapreDatasetA1(vector<double>& x, vector<double>& y, 
	vector<double>& w, int& m) {
	m = 5;
	x = { 1,2,3,4,5 };
	y = { 1.2,2.7,5.3,7.1,8.5 };
	w = { 2,1,1,3,1 };
}

//线性拟合：x、y为已知的节点数据，w为权重数据，m为结点数
void LinearFitting(vector<double> x, vector<double> y, 
	vector<double> w, int m) {
	//MA=B
	vector<vector<double>> M;//图3.5.11中的第1个矩阵，从左往右
	//矩阵初始化
	vector<double> r1 = { 0,0 };
	vector<double> r2 = { 0,0 };
	M.push_back(r1);
	M.push_back(r2);

	vector<double> A = { 0, 0 };//图3.5.11中的第2个矩阵
	vector<double> B = { 0, 0 };//图3.5.11中的第3个矩阵

	//计算出矩阵M和B
	for (int i = 0; i < m;i++) {
		M[0][0] += w[i];
		M[0][1] += w[i] * x[i];
		M[1][0] += w[i] * x[i];
		M[1][1] += w[i] * x[i] * x[i];
		B[0] += w[i] * y[i];
		B[1] += w[i] * x[i] * y[i];
	}

	//高斯消元
	GaussElimination_Sequence(M, n, B);

	//回代上三角
	CalcUpperTraiangularMatrix(M, n, B, A);

	cout << "计算结果：y="<< A[0]  << "+" << A[1] << "x" << endl;

}

void LineFittingTest() {
	vector<double> x, y, w;
	int m;
	PreapreDatasetA1(x, y, w, m);
	LinearFitting(x, y, w, m);
	
	cout << "H";
}

例题3.8：指数曲线拟合

题目介绍

代码

方式1

或许可以直接偷懒改一下输入数据，为表3.3，那么也就和上一题一样了：

void PreapreDatasetA2(vector<double>& x, vector<double>& y, 
	vector<double>& w, int& m) {
	m = 9;
	x = { 1,2,3,4,5,6,7,8,9 };
	y = { 0.25,0.35,0.438,0.573,0.648,0.725,0.84,0.947,1.04 };
	w = { 1,1,1,1,1,1,1,1,1, };
}

输出时修改一下打印方式：

cout << "计算结果：u="<< A[0]  << "+" << A[1] << "x" << endl;
cout << "计算结果：y=10^(" << A[0] << "+" << A[1] << "x)" << endl;

方式2

也可以输入表3.2的数据，根据$(x_i,y_i)$，先取对数，求出$u_i$（计算出表3.3的数据，保留的位数不同，精度也不同），再进行线性拟合。

void PreapreDatasetB1(vector<double>& x, vector<double>& y, 
	vector<double>& w, int& m) {
	m = 9;
	x = { 1,2,3,4,5,6,7,8,9 };
	y = { 1.78,2.24,2.74,3.74,4.45,5.31,6.92,8.85,10.97 };
	w = { 1,1,1,1,1,1,1,1,1, };

	for (int i = 0; i < m; i++)
	{
		y[i] = log10(y[i]);
	}
}