最小二乘拟合多项式（利用构造正交多项式的方法）C++

最新推荐文章于 2023-03-17 08:15:00 发布

weixin_45950984

最新推荐文章于 2023-03-17 08:15:00 发布

阅读量1.3k

点赞数 1

文章标签：算法 python 机器学习

本文链接：https://blog.csdn.net/weixin_45950984/article/details/129335863

版权

1，最小二乘介绍

最小二乘法（又称最优平方估计）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。设 f(x)是定义在点集𝑋 = {𝑥1, 𝑥,2 … , 𝑥𝑛}上的列表函数或者 f(x)是定义在区间 [a,b]上表达式复杂的连续函数，构造广义多项式

𝑝(𝑥) = ∑𝑐𝑖𝜑𝑖 𝑛 𝑖=1

使得

||𝑝 − 𝑓||2 2

达到最小，其中𝑐𝑖为待定参数，𝜑𝑖为一组线性无关的函数（称为基函数，它的选取与具体的问题有关）。取 p(x)作为 f(x)的近似表达式就是最优平方逼近问题。

2，最小二乘原理

在求最小二乘拟合函数时，基函数的选取十分重要，同一组数据在不同的基函数得到的拟合函数是不同的，基函数的选取将直接影响拟合的效果。本次报告中的基函数为三项递推正交基函数。

伪代码：

（1）根据所要拟合的阶数构造基函数。

（2）计算相应的系数。

（3）代入最优逼近函数：

3，计算实例

程序包含一个cpp文件和一个Polynomial.hpp的多项式类头文件。主函数为cpp文件。

运行前要将拟合的原始数据点保存在对应文件夹下的data,txt文件里面，然后运行程序。同样，data文件储存路径可以修改（全局变量path）。拟合多项式的次数为全局变量N，可以自行修改。

本算例参考课本176页计算实习5.1，数据如表2.1所示。

表2-1 多项式拟合数据

Xi	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
Yi	5.1234	5.3057	5.5687	5.9375	6.437	7.0978	7.9493	9.0253	10.362

通过程序的计算拟合最终得到多项式：

0.990967x4+3.00333x3+2.01065x2+0.992689x+5.00097

拟合误差为：5.744292×10 -4

图 2‑1是拟合之后的效果图，可以看出效果比较理想。

4，代码

#include<iostream>
#include<Map>
#include<cmath>
#include"Polynomial.hpp"
#include<fstream>
#include<string>
#include<vector>

using namespace std;

string filename = "data.txt";
int N=4;//最高次数

//加载数据
void LoadData(vector<vector<double>> &data)
{
	ifstream ifs;
	ifs.open(filename, ios::in);
	if (ifs.fail())
	{
		cout << "打开文件失败！" << endl;
	}
	else
	{
		int i = 0, j = 0;
		double a;
		vector<double> d;
		while (!ifs.eof())
		{
			ifs >> a;
			d.push_back(a);
			i++;
			if (i == 2)
			{
				data.push_back(d);
				i = 0;
				d.erase(d.begin(), d.end());
			}
		}
	}
}

double InnerProduct(Poly& P, vector<vector<double>>& data)
{
	int n = data.size();
	double gg = 0, gf = 0,c;
	for (int i = 0; i < n; i++)
	{
		gg += P.Yvale(data[i][0]) * P.Yvale(data[i][0]);
	}
	for (int i = 0; i < n; i++)
	{
		gf += P.Yvale(data[i][0]) * data[i][1];
	}
	//cout << gg << "  " << gf << endl;
	c = gf / gg;
	return c;
}

double InnerProduct1(Poly& P1, Poly& P2, vector<vector<double>>& data)
{
	int n = data.size();
	double gg = 0;
	for (int i = 0; i < n; i++)
	{
		gg+= P1.Yvale(data[i][0]) * P2.Yvale(data[i][0]);
	}
	return gg;
}

Poly LeastSquaresFitting(int& num, vector<vector<double>>& data)
{
	Poly P;
	//计算前两项g0，g1

	double a=0, r=0;
	double c0[1][2] = { 0,1 };
	Poly g0(1, c0);
	int n = data.size();
	for (int i = 0; i < n; i++)
	{
		a += data[i][0];
	}
	double b = a / double(n);
	double c1[2][2] = { 1,1,0,-b };
	Poly g1(2, c1);
	//g0.PrintPloy();
	//g1.PrintPloy();

	//多项式前两项
	double C0 = 0, C1 = 0;
	C0 = InnerProduct(g0, data);
	C1 = InnerProduct(g1, data);
	//cout << c0 << "  " << c1 << endl;
	Poly P1;
	P = C0 * g0;
	//P.PrintPloy();
	P1= C1 * g1;
	P = P + P1;

	//第三项之后
	if (num > 1)
	{
		Poly gk = g1, gk0 = g0;
		double x1[1][2] = { 1, 1 };
		Poly Px(1, x1);
		for (int i = 1; i < num; i++)
		{
			Poly gk1;
			double beta, gama,gama0,b=0,c=0;
			Poly xgk;
			xgk = Px * gk;
			beta = InnerProduct1(xgk, gk, data);
			gama = InnerProduct1(gk, gk,data);
			gama0= InnerProduct1(gk0, gk0, data);
			b = beta / gama;
			c = gama / gama0;
			//构造正交多项式
			double pb[2][2] = { 1,1,0,-b };
			Poly x_b(2, pb);
			Poly cgk0, xbg;
			cgk0 = c * gk0;
			xbg = x_b * gk;
			gk1 = xbg - cgk0;
			//计算第n项系数
			double cn = InnerProduct(gk1, data);
			Poly Pn;
			Pn = cn * gk1;

			P = P + Pn;

			//更新
			gk0 = gk;
			gk = gk1;
		}
	}

	return P;

}

int main()
{
	//int n = 3;
	//double ploy[3][2] = {
	//	1,1,
	//2, 5,
	//3,4
	//};//指数，系数

	//Poly P(n, ploy);
	//Poly P1(n, ploy);

	//Poly P2;
	//P2= ploy[2][0] * P;
	//P1.PrintPloy();
	//P2.PrintPloy();

	//cout << P.Yvale(1.5) << endl;

	//加载数据
	vector<vector<double>> data;
	LoadData(data);
	//cout << data.size()<<endl;

	//最小二乘拟合
	Poly P;
	P = LeastSquaresFitting(N, data);
	P.PrintPloy();



	system("pause");
	return 0;

}

代码中的多项式类可以参考：

https://blog.csdn.net/weixin_45950984/article/details/129335421?spm=1001.2014.3001.5501