N元线性函数拟合的C++实现

一元线性方程可以看做是多元函数的特例，现在用矩阵形式表述多元函数情况下，最小二乘的一般形式：

 设拟合多项式为：

各店到这条曲线的距离之和，即偏差平方和如下：

对等式右边求ai的偏导数，得到：

                       ......

把这些等式表示成矩阵的形式，就可以得到下面的矩阵：

（3）

进行化简计算：

上面公式（3）可以写为：

代码如下：

#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'时，是左右两边同乘ratio
	double 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);
		}
		else
			temp.x = num;
		i++;
	}
	fclose(fp);
	return dst;
}

 XY.txt内容：

0    1.0
0.25 1.28
0.5  1.65
0.75 2.12
1    2.72

 

注：该博文是参考网友的博文https://blog.csdn.net/piaoxuezhong/article/details/54973750，知识点全是他总结的，代码修改了部分，增加了打印变换前后的矩阵。