最小二乘法及算法实现

最小二乘法

最小二乘法是一种优化方法。通过最小化误差的平方和来寻找数据的最佳函数进行匹配。

线性函数模型：

线性函数模型：
$$Y={\hat{B}}_0+{\hat{B}}_{1}X$$

残差形式写为：
$$Y_i={\hat{B}}_0+{\hat{B}}_1{X}_1+e_i$$

可将$e_i$写为
$$e_i=Y_i-{\hat{B}}_0-{\hat{B}}_1{X}_1$$

e_i为样本$(x_i,y_i)$的误差。
所以，平方损失函数可以表示为
$$Q=\sum _{i=1}^n{e}_i^2=\sum _{i=1}^n(Y_i-{\hat{Y}}_i{)}^2=\sum _{i=1}^n(Y_i-{\hat{B}}_0-{\hat{B}}_1{X}_i)$$
即，使Q最小确定直线，Q可看作是以${\hat{B}}_0$,${\hat{B}}_1$为变量的Q的函数。
问题转换成一个极值问题：
再对Q求偏导。
$$\{\begin{array}{l}\frac{\partial Q}{\partial {\hat{B}}_0}=2{\sum }_{i=1}^n(Y_i-{\hat{B}}_0-{\hat{B}}_1{X}_i)\ast (-1)=0\\ \frac{\partial Q}{\partial {\hat{B}}_1}=2{\sum }_{i=1}^n(Y_i-{\hat{B}}_0-{\hat{B}}_1{X}_i)\ast (-X_i)=0\end{array}$$
解得
$$\{\begin{array}{l}{\hat{B}}_0=\frac{\sum X_i^2\sum Y_i-\sum X_i\sum X_i{Y}_i}{n\sum X_i^2-(\sum X_i{)}^2}\\ {\hat{B}}_1=\frac{n\sum X_i{Y}_i-\sum X_i\sum Y_i}{n\sum X_i^2-(\sum X_i{)}^2}\end{array}$$

矩阵表达形式

$$Y={\hat{B}}_0+{\hat{B}}_{1}X$$

推广到一般情况下，假如有更多的模型变量$x^1,x^2,...,x^m$(指样本里的模型相关的变量),可以用线性函数表示如下：
$$y(x^1,...,x^m;B_0,B_1,...,B_m)=B_0+B_1{x}^1+...+B_m{x}^m$$

对于n个样本来说，可以用如下线性方程组表示：

$$B_0+B_1{x}_1^1+...+B_m{x}_1^m=y_1$$

$$B_0+B_1{x}_2^1+...+B_m{x}_2^m=y_2$$

$$...$$

$$B_0+B_1{x}_i^1+...+B_m{x}_i^m=y_i$$

$$...$$

$$B_0+B_1{x}_n^1+...+B_m{x}_n^m=y_n$$

将上式记为矩阵形式为：
$$\left[\begin{array}{cccc}1& x_1^1& ...& x_1^m\\ 1& x_2^1& ...& x_2^m\\ ...& ...& ...& ...\\ 1& x_n^1& ...& x_n^m\end{array}\right]\left[\begin{array}{c}{B}_0\\ B_1\\ ...\\ B_m\end{array}\right]=\left[\begin{array}{c}y\\ y\\ ...\\ y_n\end{array}\right]$$
最终形式为$$A\cdot B=Y$$

最小二乘形式，可以表示为：$$min||AB-Y|{|}_2$$

最优解为：$$\hat{B}=(A^{T}A{)}^{-1}{A}^{T}Y$$

代码

/*
最小二乘法的实现
C++版
命令行输入数据文件
最后输入x得到预测的y值
*/
#include<iostream>
#include<fstream>
#include<vector>
using namespace std;

class LeastSquare {
	double b0, b1;
public:
	LeastSquare(const vector<double>& x, const vector<double>& y)
	{
		double t1 = 0, t2 = 0, t3 = 0, t4 = 0;
		for (int i = 0; i<x.size(); ++i)
		{
			t1 += x[i] * x[i];
			t2 += x[i];
			t3 += x[i] * y[i];
			t4 += y[i];
		}
		
		b0 = (t1*t4 - t2*t3) / (t1*x.size() - t2*t2);        // 求得 B0
		b1 = (t3*x.size() - t2*t4) / (t1*x.size() - t2*t2);  // 求得 B1 
	}

	double getY(const double x) const
	{
		return b0+b1*x;
	}

	void print() const
	{
		if (b1>=0)
			cout << "y = " << b0 << "+" << b1 << 'x' << "\n";
		else
			cout << "y = " << b0 << "" << b1 << 'x' << "\n";
	}

};

int main(int argc, char *argv[])
{
	if (argc != 2)
	{
		cout << " data.txt don't exit " << endl;
		return -1;
	}
	else
	{
		vector<double> x;
		vector<double> y;
		int count = 1;
		ifstream in(argv[1]);
		for (double d; in >> d; count++)
			if (count % 2 == 1)
				x.push_back(d);
			else
				y.push_back(d);
		LeastSquare ls(x, y);
		ls.print();

		cout << "Input x:\n";
		double x0;
		while (cin >> x0)
		{
			cout << "y = " << ls.getY(x0) << endl;
			cout << "Input x:\n";
		}
	}
	int endline;
	cin >> endline;
}

int main(int argc,char* argv[])
argc是命令行总的参数个数，argv[]是argc个参数，其中第0个参数是程序的全名，以后的参数命令行后面跟的用户输入的参数 比如：
int main(int argc, char* argv[])
{
}

两种方法：
第一种，无需调试的情况：
直接用dos命令进入到.exe目录下然后输入：*.exe pra1 pra2
第二种，需要调试的情况：

1. 先选择项目－〉右键－〉属性

2. 调试 －〉命令行参数

   在命令行参数里面输入命令行参数即可。
   需要注意的是，不需要像第一种那样样输入***.exe了。只需要输入 pra1 pra2 ，中间用空格隔开。

   例如：以上实现代码，需要输入一个data.txt，输入格式是(x,y)的点值。
   1 0
   2 1
   3 2
   0 1
   1 2
   2 3