利用Eigen完成一元线性回归

        学完Fortran之后，有一段时间人很麻，因为这个东西，看起来好用，但是吧，第一我发现这压根就没有啥库，IML还是IMFL，这个东西intel官网可以下载，但是我C盘空间只剩4G了，就光这一个IML就得10G，想了想，感觉自己好像点错了科技树，再回过头去看自己大学学的东西，python就混了个脸熟，具体让我用它干啥我搞不出来，java，嘿嘿，不瞒你们，学完就忘了，lz又不喜欢matlab，想来想去，现在能上手的就是c++了，所以又重新捡起来c++，又发现Eigen库和GSL库，完全可以帮助我完成大规模的矩阵计算，再加上matplotlibcpp的加持，完全可以解决很多问题，所以再就没管Fortran的事情，之后完成老师给的小任务后，我就再不搞Fortran了，还是c++香甜可口。

        说了这么多，言归正传，研一嘛，还是要上课的，刚好数量统计上到一元线性回归了，好熟悉的感觉，高中的时候，那会我爸老怕我拿计算器作弊，还振振有词说要锻炼我的手算能力，就没给买，上数学课的时候老师讲一元线性回归的时候，算的lz手麻，作业最后还是抄了好朋友的，今天看到了，跟着书上的推了推公式，想着用Eigen来编程实现一下，感觉做出来效果还很不错，哈哈，废话不多说，上公式：

 现在要求的就是三个量、、，公式在手，code起来。

这里面的io方法我用的是(6条消息) 【C++】读取 .txt 文件中的一列或多列数据(非常实用)_一米九零小胖子的博客-CSDN博客_c++读取txt文件每一行

这个博文里的读取多列的方法，我觉得写的非常好，题目用的是西安交大数理统计书给出的例题。

这是我的代码（没怎么写注释，我觉得对照着变量名和公式应该都能看懂的）：

#include <iostream>
#include<Eigen>
#include<vector>
#include<fstream>
#include<sstream>
#include"matplotlibcpp.h"
using namespace std;
using namespace Eigen;
namespace plt = matplotlibcpp;
int main() {
	fstream ifs;
	string filepath = "C:/Users/L/Desktop/data.txt";
	vector<double> x, y;
	vector<string> item;			
	string temp;						
	ifs.open(filepath,ios::in);
	while (getline(ifs, temp))      
	{
		item.push_back(temp);
	}
	for (auto it = item.begin(); it != item.end(); it++)
	{
		istringstream istr(*it); 
		string str;
		int count = 0;	
		double r;
		while (istr >> str)       
		{
			if (count == 0)
			{
				 r = atof(str.c_str());      
				x.push_back(r);
			}
			//获取第2列数据
			else if (count == 1)
			{
				r = atof(str.c_str());      
				y.push_back(r);
			}
			count++;
		}	
		
	}
	ifs.close();
	VectorXd xba(x.size()), yba(x.size()),sei1(x.size()), sei2(x.size());
	//初始化
	for (int i = 0; i < x.size(); i++) {
		xba(i) = x[i];
		yba(i) = y[i];
	}
	double alpha, beta,xmean,ymean,sigma, xmeanba,seita, seita1,seita2;
	xmean = xba.array().mean();
	ymean= yba.array().mean();
	for (int i = 0; i < x.size(); i++) {
		xba(i) = xba(i) - xmean;
		yba(i) = yba(i) - ymean;
	}
	sigma = xba.transpose() * yba;
	xmeanba = xba.transpose() * xba;
	beta = sigma / xmeanba;
	alpha = ymean - xmean * beta;
	for (int i = 0; i < x.size(); i++) {
		sei1(i) = x[i] * x[i];
		sei2(i) = y[i] * y[i];
	}
	seita1 = (sei1.array().sum() / x.size() - xmean * xmean);
	seita2 = (sei2.array().sum() / x.size() - ymean * ymean);
	seita = seita2 - (beta * beta) * seita1;
	vector<double> yimitate;
	for (int i = 0; i < x.size(); i++) {
		yimitate.push_back(alpha + beta * x[i]);
	}
	fstream ofs;
	string outfile = "G:/today/Univariate linear regression.txt";
	ofs.open(outfile,ios::out);
	if (ofs.is_open()) {
            ofs << "X" << '\t' <<"Y" << '\t' << "Yimitate" << '\t' << "abs" << endl;
		for (int i = 0; i < x.size(); i++) {
			ofs << x[i] << '\t' << y[i] << '\t' << yimitate[i] <<'\t'<<abs(y[i]-yimitate[i]) << endl;
		}
	}
	ofs.close();
	cout << " y= " << alpha << " + " << beta << "x" << endl;
	cout << "sigma*sigma=" << seita<< endl;
	plt::named_plot("y=a+bx", x, yimitate,"r-");
	plt::named_plot("original data", x, y, "bo");
	plt::xlabel("x");
	plt::ylabel("y");
	plt::title("Univariate linear regression");
	plt::legend();
	plt::show();
	return 0;
}

这是原始数据：

最后看看效果：

贴个大图：

还有写入的数据：

 

 完美，至于Eigen和matplotlibcpp的安装嘛，嘿嘿，不讲，Csdn上多的是。

行了，今天就到这里吧，great！！！